1000RNASEQ chicken ASE: change bam index output name, use .bam.bai for phASER...

1000RNASEQ chicken ASE: change bam index output name, use .bam.bai for phASER and .bai for ASEReadCounter, change phASER output csv file

Snakemake/1000RNASeq_chicken/ASE/Snakefile

@@ -138,14 +138,14 @@ for sample in table["sample_name"].unique():
 # # counter
 for sample in table["sample_name"].unique():
 	final_outputs.append('Results/Counter/' + sample + '_ASEReadCounter.csv')
-	final_outputs.append('Results/Counter/' + sample + '_phASER.csv')
+	final_outputs.append('Results/Counter/' + sample + '_phASER.allelic_counts.txt')
 
 # print(final_outputs)
 # ['Results/genomeMasked/reference_masked.fa', 
 #    'Results/STAR_Aln_2/samplePE_rg_genomic_pp_uniq_rmdup_split.bam', 
 #    'Results/STAR_Aln_2/sampleSE_rg_genomic_uniq_rmdup_split.bam', 
 #    'Results/Counter/samplePE_ASEReadCounter.csv', 'Results/Counter/sampleSE_ASEReadCounter.csv', 
-#    'Results/Counter/samplePE_phASER.csv', 'Results/Counter/sampleSE_phASER.csv']
+#    'Results/Counter/samplePE_phASER.allelic_counts.txt', 'Results/Counter/sampleSE_phASER.allelic_counts.txt']
 
 rule all:
     input:

Snakemake/1000RNASeq_chicken/ASE/rules/ASE_counter.smk

@@ -28,10 +28,10 @@ def get_splitNCigar_bai(wildcards):
     
     # paired end, filter on properly paired and rmdup
     if len(sample_table["forward_read"]) == len(sample_table["reverse_read"].dropna()) :
-        return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_pp_rmdup_uniq_split.bai"
+        return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_pp_rmdup_uniq_split.bam.bai"
     # single end, no filter
     elif len(sample_table["reverse_read"].dropna() ) == 0:
-        return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_rmdup_uniq_split.bai"
+        return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_rmdup_uniq_split.bam.bai"
     else :
     	raise Exception(wildcards.sample + " is sequenced in pair end & single end mode\n 	This workflow do not accept this case!\n")
 
@@ -39,8 +39,8 @@ rule ASEReadCounter:
 	input:
 		ref = config["fasta_ref"],
 		vcf = 'Results/hard_filtering/' + os.path.basename(config['vcf']).replace('.vcf.gz','') + '_SNP_filtered.vcf.gz',
-		bam = get_splitNCigar_bam ,
-		bai = get_splitNCigar_bai 
+		bam = get_splitNCigar_bam,
+		bai = get_splitNCigar_bai
 	output:
 		cvs = 'Results/Counter/{sample}_ASEReadCounter.csv'
 	log:
@@ -52,6 +52,7 @@ rule ASEReadCounter:
    		baseQuality = config['baseQuality']
 	shell:
 		"""
+        ln -s {input.bam}.bai `sed 's/.bam.bai/.bai/'`
 		java -Xmx{params.mem} -jar {params.jar} -T ASEReadCounter -R {input.ref} -o {output} -I {input.bam} -sites {input.vcf}  -U ALLOW_N_CIGAR_READS --minMappingQuality {params.mappingQuality} --minBaseQuality {params.baseQuality} 2> {log}
 		"""
 
@@ -61,7 +62,7 @@ rule phASER:
 		bam = get_splitNCigar_bam,
 		bai = get_splitNCigar_bai
 	output:
-		cvs = 'Results/Counter/{sample}_phASER.csv'
+		csv = 'Results/Counter/{sample}_phASER.allelic_counts.txt'
 	log:
 	threads:  config["phASER"]["cpu"]
 	params:

Snakemake/1000RNASeq_chicken/ASE/rules/basics.smk

@@ -35,7 +35,7 @@ rule bam_index:
         "{file}.bam"
         
     output:
-        temp("{file}.bai")
+        temp("{file}.bam.bai")
         
     shell:
         "samtools index {input} "