added rules to generate normalised coverage track #45

Snakefile

@@ -103,7 +103,19 @@ rule finish:
             ]
         ),
         salmon_merge_genes = expand(os.path.join(config["output_dir"], "summary_salmon", "quantmerge", "genes_{salmon_merge_on}.tsv"), salmon_merge_on=["tpm", "numreads"]),
-        salmon_merge_transcripts = expand(os.path.join(config["output_dir"], "summary_salmon", "quantmerge", "transcripts_{salmon_merge_on}.tsv"), salmon_merge_on=["tpm", "numreads"])
+        salmon_merge_transcripts = expand(os.path.join(config["output_dir"], "summary_salmon", "quantmerge", "transcripts_{salmon_merge_on}.tsv"), salmon_merge_on=["tpm", "numreads"]),
+        star_rpm = expand(os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str1.out.bg"),
+            zip,
+            sample=[i for i in list(samples_table.index.values)],
+            seqmode=[
+                samples_table.loc[i, 'seqmode']
+                for i in list(samples_table.index.values)
+            ])
 
 rule create_index_star:
     """Create index for STAR alignments"""

tests/test_integration_workflow/expected_output.md5

@@ -78,3 +78,11 @@ e72f5d798c99272f8c0166dc77247db1  results/single_end/synthetic_10_reads_mate_1_s
 92bcd0592d22a6a58d0360fc76103e56  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/salmon_quant/aux_info/observed_bias
 92bcd0592d22a6a58d0360fc76103e56  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/salmon_quant/aux_info/observed_bias_3p
 d41d8cd98f00b204e9800998ecf8427e  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/salmon_quant/aux_info/unmapped_names.txt
+0139e75ddbfe6eb081c2c2d9b9108ab4  results/paired_end/synthetic_10_reads_paired_synthetic_10_reads_paired/STAR_coverage/synthetic_10_reads_paired_synthetic_10_reads_paired_Signal.UniqueMultiple.str1.out.bg
+c266d31e0a2ad84975cb9de335891e64  results/paired_end/synthetic_10_reads_paired_synthetic_10_reads_paired/STAR_coverage/synthetic_10_reads_paired_synthetic_10_reads_paired_Signal.UniqueMultiple.str2.out.bg
+0139e75ddbfe6eb081c2c2d9b9108ab4  results/paired_end/synthetic_10_reads_paired_synthetic_10_reads_paired/STAR_coverage/synthetic_10_reads_paired_synthetic_10_reads_paired_Signal.Unique.str1.out.bg
+c266d31e0a2ad84975cb9de335891e64  results/paired_end/synthetic_10_reads_paired_synthetic_10_reads_paired/STAR_coverage/synthetic_10_reads_paired_synthetic_10_reads_paired_Signal.Unique.str2.out.bg
+ea91b4f85622561158bff2f7c9c312b3  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/STAR_coverage/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1_Signal.UniqueMultiple.str1.out.bg
+ede14ac41c10067838f375106fce4852  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/STAR_coverage/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1_Signal.UniqueMultiple.str2.out.bg
+ea91b4f85622561158bff2f7c9c312b3  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/STAR_coverage/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1_Signal.Unique.str1.out.bg
+ede14ac41c10067838f375106fce4852  results/single_end/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1/STAR_coverage/synthetic_10_reads_mate_1_synthetic_10_reads_mate_1_Signal.Unique.str2.out.bg

workflow/rules/paired_end.snakefile.smk

@@ -338,3 +338,59 @@ rule pe_genome_quantification_kallisto:
         {params.directionality} \
         {input.reads1} {input.reads2} > {output.pseudoalignment}) &> {log}"
     
+rule star_rpm_paired_end:
+    ''' Create stranded bedgraph coverage with STARs RPM normalisation '''
+    input: 
+        bam = os.path.join(
+            config["output_dir"],
+            "paired_end",
+            "{sample}",
+            "map_genome",
+            "{sample}_Aligned.sortedByCoord.out.bam")
+    output:
+        str1 = (os.path.join(
+            config["output_dir"],
+            "paired_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.Unique.str1.out.bg"),
+            os.path.join(
+            config["output_dir"],
+            "paired_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str1.out.bg")),
+        str2 = (os.path.join(
+            config["output_dir"],
+            "paired_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.Unique.str2.out.bg"), 
+            os.path.join(
+            config["output_dir"],
+            "paired_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str2.out.bg"))
+    params:
+        out_dir = lambda wildcards, output: os.path.dirname(output.str1[0]),
+        prefix = lambda wildcards, output: os.path.join(os.path.dirname(output.str1[0]),
+            str(wildcards.sample) + "_"),
+        stranded = "Stranded"
+    singularity:
+        "docker://zavolab/star:2.6.0a"
+    log: os.path.join(config["log_dir"], "paired_end", "{sample}", "star_rpm_paired_end.log")
+    threads: 4
+    shell:
+        """
+        (mkdir -p {params.out_dir}; \
+        chmod -R 777 {params.out_dir}; \
+        STAR \
+        --runMode inputAlignmentsFromBAM \
+        --runThreadN {threads} \
+        --inputBAMfile {input.bam} \
+        --outWigType "bedGraph" \
+        --outWigStrand "{params.stranded}" \
+        --outWigNorm "RPM" \
+        --outFileNamePrefix {params.prefix}) &> {log}
+        """

workflow/rules/single_end.snakefile.smk

@@ -264,3 +264,59 @@ rule genome_quantification_kallisto:
         {params.directionality} \
         {input.reads} > {output.pseudoalignment}) &> {log}"
 
+rule star_rpm_single_end:
+    ''' Create stranded bedgraph coverage with STARs RPM normalisation '''
+    input: 
+        bam = os.path.join(
+            config["output_dir"],
+            "single_end",
+            "{sample}",
+            "map_genome",
+            "{sample}_Aligned.sortedByCoord.out.bam")
+    output:
+        str1 = (os.path.join(
+            config["output_dir"],
+            "single_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.Unique.str1.out.bg"),
+            os.path.join(
+            config["output_dir"],
+            "single_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str1.out.bg")),
+        str2 = (os.path.join(
+            config["output_dir"],
+            "single_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.Unique.str2.out.bg"),
+            os.path.join(
+            config["output_dir"],
+            "single_end",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str2.out.bg"))
+    params:
+        out_dir = lambda wildcards, output: os.path.dirname(output.str1[0]),
+        prefix = lambda wildcards, output: os.path.join(os.path.dirname(output.str1[0]),
+            str(wildcards.sample) + "_"),
+        stranded = "Stranded"
+    singularity:
+        "docker://zavolab/star:2.6.0a"
+    log: os.path.join(config["log_dir"], "single_end", "{sample}", "star_rpm_single_end.log")
+    threads: 4
+    shell:
+        """
+        (mkdir -p {params.out_dir}; \
+        chmod -R 777 {params.out_dir}; \
+        STAR \
+        --runMode inputAlignmentsFromBAM \
+        --runThreadN {threads} \
+        --inputBAMfile {input.bam} \
+        --outWigType "bedGraph" \
+        --outWigStrand {params.stranded} \
+        --outWigNorm "RPM" \
+        --outFileNamePrefix {params.prefix}) &> {log}
+        """