Merge pull request #7 from avilab/recalib

Recalib
hivlab · Dec 21, 2020 · 284627e · 284627e
2 parents 9dae03e + a80fb34
commit 284627e
Showing 1 changed file with 74 additions and 34 deletions.
diff --git a/workflow/Snakefile b/workflow/Snakefile
@@ -76,7 +76,7 @@ rule reformat:
     input:
         unpack(get_fastq),
     output:
-        out=temp("results/{sample}/{run}/interleaved.fq"),
+        out="results/{sample}/{run}/interleaved.fq",
     log:
         "results/{sample}/{run}/log/reformat.log",
     params:
@@ -95,7 +95,7 @@ rule trim:
     input:
         input=rules.reformat.output.out,
     output:
-        out=temp("results/{sample}/{run}/trimmed.fq"),
+        out="results/{sample}/{run}/trimmed.fq",
     log:
         "results/{sample}/{run}/log/trim.log",
     params:
@@ -130,15 +130,15 @@ rule correct1:
     input:
         input=rules.filter.output.out,
     output:
-        out=temp("results/{sample}/{run}/ecco.fq"),
+        out="results/{sample}/{run}/ecco.fq",
     params:
         extra="ecco mix vstrict ordered",
     log:
         "results/{sample}/{run}/log/correct1.log",
     resources:
         runtime=120,
         mem_mb=4000,
-    threads: 8
+    threads: 4
     wrapper:
         f"{WRAPPER_PREFIX}/v0.6/bbtools/bbmerge"
 
@@ -155,6 +155,7 @@ rule correct2:
     resources:
         runtime=120,
         mem_mb=lambda wildcards, input: round(4000 + 6 * input.size_mb),
+    threads: 4
     wrapper:
         f"{WRAPPER_PREFIX}/v0.6/bbtools/tadpole"
 
@@ -214,7 +215,7 @@ rule samtools_sort:
     resources:
         mem_mb=4000,
         runtime=lambda wildcards, attempt: attempt * 240,
-    threads: 8
+    threads: 4
     wrapper:
         "0.68.0/bio/samtools/sort"
 
@@ -237,64 +238,103 @@ rule mark_duplicates:
         "0.68.0/bio/picard/markduplicates"
 
 
-rule realign:
+rule lofreq1:
+    """
+    Variant calling.
+    """
     input:
         ref=REF_GENOME,
         bam=rules.mark_duplicates.output.bam,
     output:
-        "results/{sample}/realign.bam",
+        "results/{sample}/lofreq1.vcf",
     log:
-        "results/{sample}/log/realign.log",
+        "results/{sample}/log/lofreq1.log",
     params:
-        extra="--verbose --keepflags",
+        extra="",
     resources:
         runtime=120,
         mem_mb=4000,
-    threads: 8
+    threads: 4
     wrapper:
-        f"{WRAPPER_PREFIX}/master/lofreq/viterbi"
+        f"{WRAPPER_PREFIX}/v0.6/lofreq/call"
 
 
-rule indelqual:
+rule indexfeaturefile:
     """
-    Indel recalibration.
+    Index vcf vile.
     """
     input:
-        ref=REF_GENOME,
-        bam=rules.realign.output[0],
+        "results/{sample}/lofreq1.vcf",
     output:
-        "results/{sample}/indelqual.bam",
+        "results/{sample}/lofreq1.vcf.idx",
     log:
-        "results/{sample}/log/indelqual.log",
+        "results/{sample}/log/indexfeaturefile.log",
     params:
-        extra="--verbose",
+        extra="",
     resources:
         runtime=120,
         mem_mb=4000,
-    threads: 8
+    threads: 1
     wrapper:
-        f"{WRAPPER_PREFIX}/v0.6/lofreq/indelqual"
+        f"{WRAPPER_PREFIX}/v0.6.1/gatk/indexfeaturefile"
 
 
-rule alnqual:
+rule gatk_baserecalibrator:
     input:
         ref=REF_GENOME,
-        bam=rules.indelqual.output[0],
+        bam=rules.mark_duplicates.output.bam,
+        dict=REF_GENOME_DICT,
+        known="results/{sample}/lofreq1.vcf",
+        feature_index=rules.indexfeaturefile.output[0],
     output:
-        "results/{sample}/alnqual.bam",
+        recal_table="results/{sample}/recal_table.grp",
     log:
-        "results/{sample}/log/alnqual.log",
-    params:
-        extra="-b",
+        "results/{sample}/log/baserecalibrator.log",
+    resources:
+        runtime=120,
+        mem_mb=4000,
+    wrapper:
+        "0.68.0/bio/gatk/baserecalibrator"
+
+
+rule applybqsr:
+    """
+    Inserts indel qualities into BAM.
+    """
+    input:
+        ref=REF_GENOME,
+        bam=rules.mark_duplicates.output.bam,
+        recal_table="results/{sample}/recal_table.grp",
+    output:
+        bam="results/{sample}/recalibrated.bam",
+    log:
+        "results/{sample}/log/applybqsr.log",
     resources:
         runtime=120,
         mem_mb=4000,
-    threads: 8
     wrapper:
-        f"{WRAPPER_PREFIX}/v0.7.0/lofreq/alnqual"
+        "0.68.0/bio/gatk/applybqsr"
 
 
-LOFREQ_OPTS = "--call-indels --min-cov 50 --max-depth 1000000  --min-bq 30 --min-alt-bq 30 --min-mq 20 --max-mq 255 --min-jq 0 --min-alt-jq 0 --def-alt-jq 0 --sig 0.01 --bonf dynamic --no-default-filter"
+rule indelqual:
+    """
+    Insert indel qualities.
+    """
+    input:
+        ref=REF_GENOME,
+        bam=rules.applybqsr.output.bam,
+    output:
+        "results/{sample}/indelqual.bam",
+    log:
+        "results/{sample}/log/indelqual.log",
+    params:
+        extra="--verbose",
+    resources:
+        runtime=120,
+        mem_mb=4000,
+    threads: 8
+    wrapper:
+        f"{WRAPPER_PREFIX}/v0.6/lofreq/indelqual"
 
 
 rule lofreq:
@@ -303,13 +343,13 @@ rule lofreq:
     """
     input:
         ref=REF_GENOME,
-        bam=rules.alnqual.output[0],
+        bam=rules.indelqual.output[0],
     output:
         "results/{sample}/lofreq.vcf",
     log:
         "results/{sample}/log/lofreq.log",
     params:
-        extra=LOFREQ_OPTS,
+        extra="--call-indels",
     resources:
         runtime=120,
         mem_mb=4000,
@@ -323,7 +363,7 @@ rule pileup:
     Calculate coverage.
     """
     input:
-        input=rules.realign.output[0],
+        input=rules.applybqsr.output.bam,
         ref=REF_GENOME,
     output:
         out="results/{sample}/covstats.txt",
@@ -350,7 +390,7 @@ rule vcffilter:
     log:
         "results/{sample}/log/vcffilter.log",
     params:
-        extra="-f 'AF > 0.5'",
+        extra="-f 'AF > 0.9'",
     resources:
         runtime=120,
         mem_mb=4000,
@@ -536,7 +576,7 @@ rule bamstats:
     Host genome mapping stats.
     """
     input:
-        rules.realign.output[0],
+        rules.applybqsr.output.bam,
     output:
         "results/{sample}/bamstats.txt",
     resources: