Skip to content
Snippets Groups Projects
Commit 8793fa4a authored by mariabernard's avatar mariabernard
Browse files

1000RNASeq_chicken ASE : perform rmdup for single end read

parent 274dafc9
No related branches found
No related tags found
No related merge requests found
......@@ -130,7 +130,7 @@ for sample in table["sample_name"].unique():
final_outputs.append("Results/STAR_Aln_2/" + sample + "_rg_genomic_pp_rmdup_uniq_split.bam")
# single end, no filter on properly paired
elif len(sample_table["reverse_read"].dropna() ) == 0:
final_outputs.append("Results/STAR_Aln_2/" + sample + "_rg_genomic_uniq_split.bam")
final_outputs.append("Results/STAR_Aln_2/" + sample + "_rg_genomic_rmdup_uniq_split.bam")
# # counter
for sample in table["sample_name"].unique():
......@@ -146,4 +146,4 @@ for sample in table["sample_name"].unique():
rule all:
input:
final_outputs
\ No newline at end of file
final_outputs
......@@ -18,7 +18,7 @@ def get_splitNCigar_bam(wildcards):
return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_pp_rmdup_uniq_split.bam"
# single end, no filter
elif len(sample_table["reverse_read"].dropna() ) == 0:
return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_uniq_split.bam"
return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_rmdup_uniq_split.bam"
else :
raise Exception(wildcards.sample + " is sequenced in pair end & single end mode\n This workflow do not accept this case!\n")
......
......@@ -8,7 +8,7 @@ Align RNASeq reads on masked genome
'''
wildcard_constraints:
properlyPaired_rmDup = '_pp_rmdup_|_'
properlyPaired_rmDup = '_pp_|_'
STAR_INDEX_FILE=['genomeParameters.txt', 'chrName.txt', 'chrLength.txt', 'chrStart.txt', 'chrNameLength.txt', 'exonGeTrInfo.tab', 'geneInfo.tab', 'transcriptInfo.tab', 'exonInfo.tab', 'sjdbList.fromGTF.out.tab', 'sjdbInfo.txt', 'sjdbList.out.tab', 'Genome', 'SA', 'SAindex']
......@@ -161,8 +161,8 @@ rule rmdup:
input:
bam = get_properly_paired_bam
output:
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired_rmDup}.bam",
metrics = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired_rmDup}.metrics.txt",
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired}.bam",
metrics = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired}.metrics.txt",
params:
mem=str(int(config["rmdup"]["mem"].replace("G","")) -4 )+"G" if int(config["rmdup"]["mem"].replace("G","")) -4 > 1 else config["rmdup"]["mem"] ,
jar=find_jar("picard.jar")
......@@ -173,9 +173,9 @@ rule rmdup:
rule remove_multimap:
input:
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired_rmDup}.bam"
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired}rmdup.bam"
output:
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired_rmDup}uniq.bam"
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired}rmdup_uniq.bam"
shell:
"""
samtools view -h -q 255 {input.bam} -bo {output.bam}
......@@ -184,9 +184,9 @@ rule remove_multimap:
rule SplitNCigarReads:
input:
fasta = 'Results/genomeMasked/' + os.path.splitext(os.path.basename(config['fasta_ref']))[0] + '_masked.fa',
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired_rmDup}uniq.bam"
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired}rmdup_uniq.bam"
output:
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPaired_rmDup}uniq_split.bam"
bam = "Results/STAR_Aln_2/{sample}_rg_genomic{properlyPairedrmdup_uniq_split.bam"
params:
mem=config["SplitNCigarReads"]["mem"],
jar=find_jar("GenomeAnalysisTK.jar")
......@@ -194,4 +194,4 @@ rule SplitNCigarReads:
"""
java -Xmx{params.mem} -jar {params.jar} -T SplitNCigarReads -R {input.fasta} -I {input.bam} -o {output.bam} -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS
"""
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment