Snakemake RNAseqVariantCalling : report 1000RRNASeq_chicken/calling update:...

Snakemake RNAseqVariantCalling : report 1000RRNASeq_chicken/calling update: quality trimming and trimgalore_summary

Snakemake/RNAseqVariantCalling/Snakefile

@@ -168,12 +168,13 @@ rule trim_se:
     output:
         out=temp("Results/TrimGalore/{prefix}_trim.fastq.gz")
     params:
-        min_trim=lambda wildcards : str(int([ int(table[table["forward_read"] == f]["read_length"].tolist()[0]) for f in table["forward_read"] if f.startswith(wildcards.prefix)][0] / 3 ))
+        min_trim=lambda wildcards : str(int([ int(table[table["forward_read"] == f]["read_length"].tolist()[0]) for f in table["forward_read"] if f.startswith(wildcards.prefix)][0] / 3 )),
+        qual=config["trimming_quality"]
     log:
         "Results/TrimGalore/{prefix}_trimSe.log"
     shell:
         """
-        trim_galore --no_report_file --length {params.min_trim} --quality 0 -o `dirname {output.out}` {input} 2> {log}
+        trim_galore --no_report_file --length {params.min_trim} --quality {params.qual} -o `dirname {output.out}` {input} 2> {log}
         mv `dirname {output.out}`/{wildcards.prefix}*trimmed.fq.gz {output.out}
         """
 
@@ -186,10 +187,11 @@ rule trim_pe:
     log:
         "Results/TrimGalore/{prefix}_trimPe.log"
     params:
-        min_trim=lambda wildcards : str(int([ int(table[table["forward_read"] == f]["read_length"].tolist()[0]) for f in table["forward_read"] if f.startswith(wildcards.prefix)][0] / 3 ))
+        min_trim=lambda wildcards : str(int([ int(table[table["forward_read"] == f]["read_length"].tolist()[0]) for f in table["forward_read"] if f.startswith(wildcards.prefix)][0] / 3 )),
+        qual=config["trimming_quality"]
     shell:
         """
-        trim_galore --paired --no_report_file --length {params.min_trim} --quality 0 -o `dirname {output.out1}` {input} 2> {log}
+        trim_galore --paired --no_report_file --length {params.min_trim} --quality {params.qual} -o `dirname {output.out1}` {input} 2> {log}
         mv `dirname {output.out1}`/{wildcards.prefix}*_val_1.fq.gz {output.out1}
         mv `dirname {output.out1}`/{wildcards.prefix}*_val_2.fq.gz {output.out2}
         """

Snakemake/RNAseqVariantCalling/example/config_calling.yaml

@@ -20,6 +20,9 @@ gtf_ref: <ABS_PATH>/reference.gtf
 # known_vcf file is set of known variants used to recalibrate bases quality in GATK preprocessing steps RealignerTargetCreator and BaseRecalibrator
 known_vcf: <ABS_PATH>/refence_known_var.vcf
 
+# quality trimming threshold used in trimgalore to remove low quality bases.
+trimming_quality : 15
+
 # computing ressources, also give to --cluster-config snakemake option if executed on a cluster
 # this yaml file defined default resources (mem and cpu at least) in a __default__ section, and specific resources either all or one resource for particular rule if different from the default
 resources: <ABS_PATH>/resources.yaml

Snakemake/RNAseqVariantCalling/script/trimgalore_summary.py

@@ -122,11 +122,15 @@ else:
 sample_names = sample_paired + sample_single
 res_lines = OrderedDict()
 
-for head in res_paired_lines:
-    if head in res_single_lines:
-        res_lines[head] = res_paired_lines[head] + res_single_lines[head]
-    else:
-        res_lines[head] = res_paired_lines[head] + ['0']*len(sample_single)
+if len(res_paired_lines) > 0:
+    for head in res_paired_lines:
+        if head in res_single_lines:
+            res_lines[head] = res_paired_lines[head] + res_single_lines[head]
+        else:
+            res_lines[head] = res_paired_lines[head] + ['0']*len(sample_single)
+else:
+    for head in res_single_lines:
+        res_lines[head] = res_single_lines[head]
 
 # write results
 with open(snakemake.output.out, "w") as summary_file: