Added pipeline to call viral consensus.

Author: gkarthik

pipeline_consensus/Snakefile

@@ -0,0 +1,84 @@
+import os
+from datetime import datetime
+import time
+import re
+from shutil import copyfile
+
+import pandas as pd
+import re
+
+configfile: "config.json"
+
+reference = config["reference"]
+bed_file = config["bed_file"]
+out_dir = config["out_dir"]
+samples_path = config["samples_path"]
+
+df = pd.read_table(samples_path, sep="\t")
+_ = df.groupby("sample").sum()
+lib_delim = config["library_delimiter"]
+
+_["sample_library"] = _["sample_library"].apply(lambda x: x.split(lib_delim)[0] +"_" + "_".join(re.findall("L[0-9]+", x)))
+
+rule all:
+    input:
+        expand("{out_dir}/consensus_sequences/{sample}.fa", out_dir = out_dir, sample = _["sample_library"])
+
+rule call_consensus:
+    input:
+        "{out_dir}/trimmed_bams/{sample}.trimmed.sorted.bam"
+    output:
+        "{out_dir}/consensus_sequences/{sample}.fa"
+    shell:
+        """
+        samtools mpileup -A -Q 0 -d 300000 {input} | ivar consensus -p {output} -n N -m 10
+        """
+
+rule trim_reads:
+    input:
+        "{out_dir}/merged_aligned_bams/{sample}.sorted.bam"
+    output:
+        "{out_dir}/trimmed_bams/{sample}.trimmed.sorted.bam"
+    params:
+        bed="{bed}".format(bed = bed_file),
+        tmp="{out_dir}/trimmed_bams/{sample}.trimmed.bam"
+    shell:
+        """
+        ivar trim -i {input} -b {params.bed} -p {params.tmp}
+        samtools sort -T {wildcards.sample}.trim -o {output} {params.tmp}
+        rm {params.tmp}
+        """
+
+rule merge_multiple_libraries:
+    input:
+        bams=lambda wildcards: df[df["sample"] == _[_["sample_library"] == wildcards.sample].index.values[0]]["sample_library"].apply(lambda x: os.path.join(out_dir, "aligned_bams", x +".sorted.bam")).tolist(),
+        forward=lambda wildcards: df[df["sample"] == _[_["sample_library"] == wildcards.sample].index.values[0]]["forward"].sort_values().tolist(),
+        reverse=lambda wildcards: df[df["sample"] == _[_["sample_library"] == wildcards.sample].index.values[0]]["reverse"].sort_values().tolist()
+    output:
+        bam="{out_dir}/merged_aligned_bams/{sample}.sorted.bam",
+        fastq=expand("{{out_dir}}/merged_fastq/{{sample}}_R{readno}.fastq.gz", readno=[1,2])
+    params:
+        tmp="{out_dir}/merged_aligned_bams/{sample}.bam"
+    shell:
+        """
+        samtools merge {params.tmp} {input.bams}
+        samtools sort -T {wildcards.sample}.merge -o {output.bam} {params.tmp}
+        rm {params.tmp}
+        cat {input.forward} > {output.fastq[0]}
+        cat {input.reverse} > {output.fastq[1]}
+        """
+
+rule align_reads:
+    input:
+        lambda wildcards: df[df["sample_library"]==wildcards.sample][["forward", "reverse"]].values[0].tolist()
+    output:
+        temp("{out_dir}/aligned_bams/{sample}.sorted.bam")
+    params:
+        ref= "{ref}".format(ref = reference),
+        tmp="{out_dir}/aligned_bams/{sample}.sorted.tmp.bam"
+    shell:
+        """
+        bwa mem {params.ref} {input[0]} {input[1]} | samtools view -F 4 -Sb | samtools sort -T {wildcards.sample}.align -o {params.tmp}
+        samtools addreplacerg -r "ID:{wildcards.sample}" -o {output} {params.tmp}
+        rm {params.tmp}
+        """

pipeline_consensus/config.json

@@ -0,0 +1,7 @@
+{
+  "reference": "/path/to/bwa/ref",
+  "bed_file": "/path/to/primer/bed/file",
+  "out_dir": "/path/to/outdir",
+  "samples_path": "example_samples.tsv",
+  "library_delimiter": "_"
+}

pipeline_consensus/example_samples.tsv

@@ -0,0 +1,10 @@
+forward	reverse	sample	sample_library
+Z025_L1_R1.fastq.gz	Z025_L1_R2.fastq.gz	Z025	Z025_L1
+Z019_L1_R1.fastq.gz	Z019_L1_R2.fastq.gz	Z019	Z019_L1
+Z019_L2_R1.fastq.gz	Z019_L2_R2.fastq.gz	Z019	Z019_L2
+Z019_L3_R1.fastq.gz	Z019_L3_R2.fastq.gz	Z019	Z019_L3
+Z023_L1_R1.fastq.gz	Z023_L1_R2.fastq.gz	Z023	Z023_L1
+Z023_L2_R1.fastq.gz	Z023_L2_R2.fastq.gz	Z023	Z023_L2
+Z023_L3_R1.fastq.gz	Z023_L3_R2.fastq.gz	Z023	Z023_L3
+Z829_L1_R1.fastq.gz	Z829_L1_R2.fastq.gz	Z829	Z829_L1
+Z829_L2_R1.fastq.gz	Z829_L2_R2.fastq.gz	Z829	Z829_L2