Add support for gzipped fasta files (#89)

* Add support for gzipped fasta files

* Add pytest with gzipped fastas

* Adjust samtools faidx rule in snakemake

* Fixes for pylint

* Account for differing assembly names in gzip pytest

Author: lcoombe

README.md

@@ -175,6 +175,7 @@ More information can be found on our [wiki page](https://github.com/BirolLab/ntS
 
 ### Tips / Visualization <a name=tips></a>
 - To lower the peak memory usage, increase the number of hash functions (--hashes) and/or increase the false positive rate (--fpr) for the constructed Bloom filter
+- Your input fasta files can be gzipped
 - Customize parameters such as --merge, --indel, --block_size and --w_rounds for your particular input data and research questions
 - For visualizing the multi-genome output synteny blocks, check out our streamlined visualization tool [ntSynt-viz](https://github.com/BirolLab/ntSynt-viz). Examples of the scripts used to generate the figures in the [ntSynt manuscript](https://link.springer.com/article/10.1186/s12915-025-02455-w) can be found in the sub-directory [visualization_scripts](https://github.com/BirolLab/ntSynt/tree/main/visualization_scripts)
 - If you do not know the approximate sequence divergence between the input assemblies, we recommend using [Mash](https://github.com/marbl/Mash) to estimate the divergences

bin/ntsynt_run_pipeline.smk

@@ -51,7 +51,11 @@ rule faidx:
     output: "{file}.fai"
     params: benchmarking=expand("{benchmark_path} -o {{file}}.faidx.time", benchmark_path=benchmark_path) if benchmark else []
     threads: 1
-    shell: "{params.benchmarking} samtools faidx -o {output} {input.fa}"
+    run:
+        if input.fa.endswith(".gz"):
+            shell("{params.benchmarking} gunzip -c {input.fa} | samtools faidx -o {output} -")
+        else:
+            shell("{params.benchmarking} samtools faidx -o {output} {input.fa}")
 
 rule make_common_bf:
     input: refs=references

bin/ntsynt_synteny.py

@@ -12,6 +12,9 @@
 import shlex
 import subprocess
 import sys
+import gzip
+import shutil
+import tempfile
 import intervaltree
 import pybedtools
 import btllib
@@ -132,6 +135,17 @@ def get_synteny_bed_lists(paths):
 
         return synteny_beds
 
+    def prepare_fasta_for_bedtools(self, fa_filename_full):
+        """If fasta is gzipped, decompress to a temp file and return its path.
+        Returns (path, is_temp) — caller should delete path if is_temp is True."""
+        if fa_filename_full.endswith(".gz"):
+            tmp = tempfile.NamedTemporaryFile(suffix=".fa", delete=False, dir=os.getcwd()) # pylint: disable=consider-using-with
+            with gzip.open(fa_filename_full, "rb") as f_in:
+                shutil.copyfileobj(f_in, tmp)
+            tmp.close()
+            return tmp.name, True
+        return fa_filename_full, False
+
     def mask_assemblies_with_synteny_extents(self, synteny_beds, w):
         "Mask each reference assembly with determined synteny blocks"
         mx_to_fa_dict = {}
@@ -142,10 +156,15 @@ def mask_assemblies_with_synteny_extents(self, synteny_beds, w):
             bed_str = "\n".join(bed_str)
             fa_filename = self.find_fa_name(assembly)
             fa_filename_full = assembly_to_fastas[fa_filename]
-            synteny_bed = pybedtools.BedTool(bed_str, from_string=True).slop(g=f"{fa_filename}.fai",
-                                                                             l=-1*(w+self.args.k),
-                                                                             r=-1*(w+self.args.k)).sort()
-            synteny_bed.mask_fasta(fi=fa_filename_full, fo=f"{fa_filename}_masked.fa.tmp")
+            fa_tmp, is_temp = self.prepare_fasta_for_bedtools(fa_filename_full)
+            try:
+                synteny_bed = pybedtools.BedTool(bed_str, from_string=True).slop(g=f"{fa_filename}.fai",
+                                                                                l=-1*(w+self.args.k),
+                                                                                r=-1*(w+self.args.k)).sort()
+                synteny_bed.mask_fasta(fi=fa_tmp, fo=f"{fa_filename}_masked.fa.tmp")
+            finally:
+                if is_temp:
+                    os.remove(fa_tmp)
 
             # pybedtools may output multi-line fasta which breaks btllib reading, need seqtk to make single line
             with open(f"{fa_filename}_masked.fa", "w", encoding="utf-8") as fout:

tests/ntsynt_tests.py

@@ -22,12 +22,18 @@ def launch_ntSynt_fof(genome_file_list, prefix="ntSynt", k=24, w=1000, **kwargs)
     ret_code = subprocess.call(cmd)
     assert ret_code == 0
 
+def strip_gz_from_line(line):
+    "Strip .gz from the assembly name (second column) of the synteny TSV line"
+    fields = line.split("\t")
+    fields[1] = fields[1].replace(".gz", "")
+    return "\t".join(fields)
+
 def are_expected_blocks(block_file1, block_file2):
     "Check that the given synteny blocks are as expected"
     with open(block_file1, 'r', encoding="utf-8") as fin1:
         with open(block_file2, 'r', encoding="utf-8") as fin2:
             while (line1 := fin1.readline()) and (line2 := fin2.readline()):
-                assert line1 == line2
+                assert strip_gz_from_line(line1) == strip_gz_from_line(line2)
 
 def test_prep_files():
     "Prep the test files"
@@ -68,3 +74,10 @@ def test_cleanup():
         cmd = shlex.split(f"pigz {infile}")
         ret_code = subprocess.call(cmd)
         assert ret_code == 0
+
+def test_3_genomes_gz():
+    "Testing ntSynt with three input gzipped genomes"
+    genome1, genome2, genome3 = "celegans-chrII-III.fa.gz", "celegans-chrII-III.A.fa.gz", "celegans-chrII-III.B.fa.gz"
+    launch_ntSynt(genome1, genome2, genome3, k=20, prefix="celegans-A-B-ntSynt_gz", indel=500, merge=3000)
+    are_expected_blocks("celegans-A-B-ntSynt_gz.synteny_blocks.tsv",
+                        "expected_result/celegans-A-B-ntSynt.synteny_blocks.tsv")