Merge pull request #112 from martinghunt/varcall_no_trim_option

Varcall no trim option

Author: martinghunt

python/clockwork/tasks/variant_call_one_sample.py

@@ -22,4 +22,5 @@ def run(options):
         cortex_mem_height=options.mem_height,
         debug=options.debug,
         keep_bam=options.keep_bam,
+        trim_reads=not options.no_trim,
     )

python/clockwork/tests/var_call_one_sample_pipeline_test.py

@@ -35,35 +35,42 @@ def test_run(self):
             directory=os.path.join(root_outdir, "ref_dir")
         )
         ref_dir.make_index_files(ref_fa, False, True, cortex_mem_height=21)
-        var_call_out = os.path.join(root_outdir, "varcall")
-        var_call_one_sample_pipeline.run(
-            [reads1],
-            [reads2],
-            ref_dir.directory,
-            var_call_out,
-            sample_name="test_sample",
-            debug=False,
-            keep_bam=True,
-            cortex_mem_height=21,
-        )
+        for trim_reads in (True, False):
+            var_call_out = os.path.join(root_outdir, "varcall")
+            if trim_reads:
+                var_call_out += ".trim_reads"
+            else:
+                var_call_out += ".no_trim_reads"
+
+            var_call_one_sample_pipeline.run(
+                [reads1],
+                [reads2],
+                ref_dir.directory,
+                var_call_out,
+                sample_name="test_sample",
+                debug=False,
+                keep_bam=True,
+                cortex_mem_height=21,
+                trim_reads=trim_reads,
+            )
 
-        got_files = sorted(list(os.listdir(var_call_out)))
-        expect_files = [
-            "cortex.vcf",
-            "final.gvcf",
-            "final.gvcf.fasta",
-            "final.vcf",
-            "map.bam",
-            "map.bam.bai",
-            "samtools.vcf",
-        ]
-        self.assertEqual(got_files, expect_files)
+            got_files = sorted(list(os.listdir(var_call_out)))
+            expect_files = [
+                "cortex.vcf",
+                "final.gvcf",
+                "final.gvcf.fasta",
+                "final.vcf",
+                "map.bam",
+                "map.bam.bai",
+                "samtools.vcf",
+            ]
+            self.assertEqual(got_files, expect_files)
 
-        with open(os.path.join(var_call_out, "final.vcf")) as f:
-            calls = [x for x in f if not x.startswith("#")]
-        self.assertEqual(len(calls), 1)
-        fields = calls[0].split("\t")
-        self.assertEqual(fields[1], "500")
-        self.assertEqual(fields[3], "A")
-        self.assertEqual(fields[4], "T")
+            with open(os.path.join(var_call_out, "final.vcf")) as f:
+                calls = [x for x in f if not x.startswith("#")]
+            self.assertEqual(len(calls), 1)
+            fields = calls[0].split("\t")
+            self.assertEqual(fields[1], "500")
+            self.assertEqual(fields[3], "A")
+            self.assertEqual(fields[4], "T")
         utils.syscall(f"rm -r {root_outdir}")

python/clockwork/var_call_one_sample_pipeline.py

@@ -13,6 +13,7 @@ def run(
     cortex_mem_height=22,
     debug=False,
     keep_bam=False,
+    trim_reads=True,
 ):
     if len(reads1_list) != len(reads2_list):
         raise Exception(
@@ -24,6 +25,11 @@ def run(
     trimmed_reads_1 = []
     trimmed_reads_2 = []
     for i in range(len(reads1_list)):
+        if not trim_reads:
+            trimmed_reads_1.append(reads1_list[i])
+            trimmed_reads_2.append(reads2_list[i])
+            continue
+
         trimmed_reads_1.append(os.path.join(outdir, f"trimmed_reads.{i}.1.fq.gz"))
         trimmed_reads_2.append(os.path.join(outdir, f"trimmed_reads.{i}.2.fq.gz"))
         read_trim.run_trimmomatic(
@@ -44,7 +50,7 @@ def run(
         read_group=("1", sample_name),
     )
     utils.syscall(f"samtools index {rmdup_bam}")
-    if not debug:
+    if trim_reads and not debug:
         for filename in trimmed_reads_1 + trimmed_reads_2:
             os.unlink(filename)
 

python/scripts/clockwork

@@ -981,6 +981,11 @@ subparser_variant_call_one_sample.add_argument(
     action="store_true",
     help="Debug mode: do not clean up any files",
 )
+subparser_variant_call_one_sample.add_argument(
+    "--no_trim",
+    action="store_true",
+    help="Skip read trimming stage at start of pipeline",
+)
 subparser_variant_call_one_sample.add_argument(
     "ref_dir", help="Directory of reference files, made by clockwork reference_prepare"
 )