Added tigmint steps to makefile (#70)

* Added tigmint steps to makefile

* Updated targets in makefile

* Update makefile targets and README.MD

* Update test file

* Update test file names

* Correct output files names

* Update tigmint-make

Co-authored-by: Puneet Sidhu <[email protected]>

Author: PuneetSS

README.md

@@ -187,7 +187,9 @@ tigmint-make tigmint-long draft=myassembly reads=myreads span=auto G=gsize dist=
 - When aligning long reads with Minimap2, use the `-y` option to include the barcode in the BX tag of the alignments.
 - When using long reads, the minimum spanning molecule thresholds (`span`) should be no greater than 1/4 of the sequence coverage. Setting the parameter `span=auto` allows the appropriate parameter value to be selected automatically (this setting requires the parameter `G` as well).
 - When using long reads, the edit distance threshold (`nm`) is automatically set to the cut length (`cut`) to compensate for the higher error rate and length. This parameter should be kept relatively high to include as many alignments as possible.
-
+- Each Tigmint (and Tigmint-long) step can be run separately through specifying the target using `tigmint-make`. For example, the bwa index step for Tigmint with linked reads can be launched with `tigmint-make tigmint-index`
+  - Tigmint steps\targets (for linked reads): `tigmint-index`, `tigmint-align`, `tigmint-molecule`, `tigmint-cut`
+  - Tigmint-long steps\targets: `tigmint-long-estimate`, `tigmint-long-to-linked`, `tigmint-long-cut`
 # Using stLFR linked reads
 
 To use stLFR linked reads with Tigmint, you will need to re-format the reads to have the barcode in a `BX:Z:` tag in the read header.

bin/tigmint-make

@@ -146,41 +146,66 @@ draft_metrics: \
 
 # Correct misassemblies using Tigmint
 tigmint: \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.fa \
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.fa \
 	$(draft).tigmint.fa
 
+# Index the linked reads 
+tigmint-index: $(draft).fa.bwt
+
+# Align the linked reads (using bwa mem) 
+tigmint-align: $(draft).$(reads).sortbx.bam
+
+# Perform the tigmint-molecule step for linked reads 
+tigmint-molecule: $(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).bed
+
+# Perform the final tigmint-cut step for linked reads
+tigmint-cut: \
+        $(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.fa \
+        $(draft).tigmint.fa
+
 # Correct misassemblies for long reads using Tigmint
 tigmint-long: \
-	$(draft).$(reads).cut$(cut).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.fa \
+	$(draft).$(reads).cut$(cut).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.fa \
 	$(draft).cut$(cut).tigmint.fa
 
+# Estimate the tigmint parameters for long reads 
+tigmint-long-estimate: $(reads).tigmint-long.params.tsv
+
+# Covert the long reads to pseudo-linked reads, perform minimap2 alignment and the tigmint-molecule step
+tigmint-long-to-linked: $(draft).$(reads).cut$(cut).molecule.size$(minsize).dist$(dist).bed
+
+# Perform the final tigmint-cut step for long reads
+tigmint-long-cut: \
+        $(draft).$(reads).cut$(cut).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.fa \
+        $(draft).cut$(cut).tigmint.fa
+
 # Calculate assembly metrics after Tigmint
 tigmint_metrics: \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.abyss-fac.tsv \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.scaftigs.abyss-fac.tsv \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.scaftigs.$(ref).samtobreak.tsv
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.abyss-fac.tsv \
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.scaftigs.abyss-fac.tsv \
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.scaftigs.$(ref).samtobreak.tsv
 
 # Scaffold using ARCS
 arcs: \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.fa \
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.fa \
 	$(draft).tigmint.arcs.fa
 
 # Calculate assembly metrics after ARCS
 arcs_metrics: \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.abyss-fac.tsv \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.scaftigs.abyss-fac.tsv \
-	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.scaftigs.$(ref).samtobreak.tsv
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.abyss-fac.tsv \
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.scaftigs.abyss-fac.tsv \
+	$(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.scaftigs.$(ref).samtobreak.tsv
 
 # Symlink the Tigmint results
-$(draft).tigmint.fa: $(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.fa
+$(draft).tigmint.fa: $(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.fa
 	ln -sf $< $@
 
 # Symlink the long read Tigmint results
-$(draft).cut$(cut).tigmint.fa: $(draft).$(reads).cut$(cut).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.fa
+$(draft).cut$(cut).tigmint.fa: $(draft).$(reads).cut$(cut).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.fa
 	ln -sf $< $@
 
 # Symlink the ARCS results
-$(draft).tigmint.arcs.fa: $(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.fa
+$(draft).tigmint.arcs.fa: $(draft).$(reads).as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).trim$(trim).window$(window).span$(span).breaktigs.$(reads).c$c_e$e_r$r.arcs.a$a_l$l.links.fa
 	ln -sf $< $@
 
 # BWA
@@ -207,30 +232,14 @@ ifeq ($G, 0)
 endif
 endif
 
-# Segment long reads from gzipped fasta file, optionally calculating tigmint-long parameters.
-$(reads).cut$(cut).fq.gz: $(longreads) check_span_g
-ifeq ($(span), auto)
-ifeq ($(dist), auto)
-	$(gtime) $(gzip) -dc $(longreads) | $(bin)/long-to-linked -l$(cut) -m$(minsize) -g$G -s -d -o $(reads).tigmint-long.params.tsv | $(gzip) > $@
-else
-	$(gtime) $(gzip) -dc $(longreads) | $(bin)/long-to-linked -l$(cut) -m$(minsize) -g$G -s -o $(reads).tigmint-long.params.tsv | $(gzip) > $@
-endif
-else
-ifeq ($(dist), auto)
-	$(gtime) $(gzip) -dc $(longreads) | $(bin)/long-to-linked -l$(cut) -m$(minsize) -d -o $(reads).tigmint-long.params.tsv | $(gzip) > $@
-else
-	$(gtime) $(gzip) -dc $(longreads) | $(bin)/long-to-linked -l$(cut) -m$(minsize) | $(gzip) > $@
-endif
-endif
-
 # samtools
 
 # Index a FASTA file.
 %.fa.fai: %.fa
 	samtools faidx $<
 
 # Create molecule extents BED
-%.as$(as).nm$(nm).molecule.size$(minsize).bed: %.sortbx.bam
+%.as$(as).nm$(nm).molecule.size$(minsize).dist$(dist).bed: %.sortbx.bam
 	$(gtime) $(bin)/tigmint_molecule.py -a$(as) -n$(nm) -q$(mapq) -d$(dist) -s$(minsize) $< | sort -k1,1 -k2,2n -k3,3n >$@
 
 # Estimate dist parameter for tigmint-long
@@ -240,7 +249,13 @@ $(reads).tigmint-long.params.tsv: $(longreads)
 	$(bin)/tigmint_estimate_dist.py - -n $(dist_sample) -o $@'
 
 # Create molecule extents BED using cut long reads
-$(draft).$(reads).cut$(cut).molecule.size$(minsize).bed: $(longreads) $(draft).fa $(reads).tigmint-long.params.tsv check_span_g
+$(draft).$(reads).cut$(cut).molecule.size$(minsize).dist$(dist).bed: $(longreads) $(draft).fa $(reads).tigmint-long.params.tsv
+ifeq ($(span), auto)
+ifeq ($G, 0)
+	$(error Must set genome size parameter (G) to calculate span automatically)
+endif
+endif
+
 ifeq ($(dist), auto)
 	$(gtime) $(bin)/../src/long-to-linked-pe -l $(cut) -m$(minsize) -g$G -s -b $(reads).barcode-multiplicity.tsv --bx -t$t --fasta -f $(reads).tigmint-long.params.tsv $< | \
 	minimap2 -y -t$t -x map-$(longmap) --secondary=no $(draft).fa - | \

tests/test_installation/expected_outputs/test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.bed tests/test_installation/expected_outputs/test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.bed

@@ -1 +1 @@
-test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa
+test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa

tests/test_installation/expected_outputs/test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa tests/test_installation/expected_outputs/test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa

@@ -1 +1 @@
-test_contig_long.test_longreads.cut500.molecule.size2000.trim0.window1000.span10.breaktigs.fa
+test_contig_long.test_longreads.cut500.molecule.size2000.distauto.trim0.window1000.span10.breaktigs.fa

tests/test_installation/expected_outputs/test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa.bed tests/test_installation/expected_outputs/test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa.bed

@@ -66,12 +66,12 @@ def tigmint_pipeline():
     outfiles = ["test_contig.fa.amb", "test_contig.fa.ann", "test_contig.fa.bwt",
                 "test_contig.fa.fai", "test_contig.fa.pac", "test_contig.fa.sa",
                 "test_contig_long.cut500.tigmint.fa", "test_contig.tigmint.fa", "test_contig_long.fa.fai",
-                "test_contig_long.test_longreads.cut500.molecule.size2000.bed",
-                "test_contig_long.test_longreads.cut500.molecule.size2000.trim0.window1000.span10.breaktigs.fa",
-                "test_contig_long.test_longreads.cut500.molecule.size2000.trim0.window1000.span10.breaktigs.fa.bed",
-                "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.bed",
-                "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa",
-                "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa.bed",
+                "test_contig_long.test_longreads.cut500.molecule.size2000.distauto.bed",
+                "test_contig_long.test_longreads.cut500.molecule.size2000.distauto.trim0.window1000.span10.breaktigs.fa",
+                "test_contig_long.test_longreads.cut500.molecule.size2000.distauto.trim0.window1000.span10.breaktigs.fa.bed",
+                "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.bed",
+                "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa",
+                "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa.bed",
                 "test_contig.test_linkedreads.sortbx.bam",
                 "test_longreads.cut500.fq.gz", "test_longreads.tigmint-long.span.txt",
                 "test_longreads.tigmint-long.params.tsv"]
@@ -231,9 +231,9 @@ def test_pipeline(tigmint_pipeline):
         assert exp_alignment == obs_bam[i]
 
     # Other readable files
-    tigmint_outputs = ["test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.bed",
-        "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa",
-        "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.trim0.window1000.span20.breaktigs.fa.bed"]
+    tigmint_outputs = ["test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.bed",
+        "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa",
+        "test_contig.test_linkedreads.as0.65.nm5.molecule.size2000.dist50000.trim0.window1000.span20.breaktigs.fa.bed"]
     for output in tigmint_outputs:
         expected_output = "expected_outputs/" + output
         with open(expected_output) as exp:
@@ -251,9 +251,9 @@ def test_pipeline(tigmint_pipeline):
 
     # Other output files
     tigmint_long_outputs = ["test_longreads.tigmint-long.params.tsv",
-        "test_contig_long.test_longreads.cut500.molecule.size2000.bed",
-        "test_contig_long.test_longreads.cut500.molecule.size2000.trim0.window1000.span10.breaktigs.fa",
-        "test_contig_long.test_longreads.cut500.molecule.size2000.trim0.window1000.span10.breaktigs.fa.bed"]
+        "test_contig_long.test_longreads.cut500.molecule.size2000.distauto.bed",
+        "test_contig_long.test_longreads.cut500.molecule.size2000.distauto.trim0.window1000.span10.breaktigs.fa",
+        "test_contig_long.test_longreads.cut500.molecule.size2000.distauto.trim0.window1000.span10.breaktigs.fa.bed"]
     for output in tigmint_long_outputs:
         expected_output = "expected_outputs/" + output
         with open(expected_output) as exp: