Fix bedtools genomecov deprecation (coverage calculation)

autometa/common/coverage.py

@@ -69,32 +69,6 @@ def from_spades_names(records):
     return coverages.to_frame()
 
 
-def make_length_table(fasta, out):
-    """Writes a tab-delimited length table to `out` given an input `fasta`.
-
-    Parameters
-    ----------
-    fasta : str
-        </path/to/assembly.fasta>
-    out : str
-        </path/to/lengths.tsv>
-
-    Returns
-    -------
-    str
-        </path/to/lengths.tsv>
-    """
-    seqs = {record.id: len(record) for record in SeqIO.parse(fasta, "fasta")}
-    lengths = pd.Series(seqs, name="length")
-    lengths.index.name = "contig"
-    lengths.to_csv(out, sep="\t", index=True, header=True)
-    return out
-
-
-def normalize(df):
-    raise NotImplementedError
-
-
 @utilities.timeit
 def get(
     fasta,
@@ -105,7 +79,6 @@ def get(
     se_reads=None,
     sam=None,
     bam=None,
-    lengths=None,
     bed=None,
     cpus=1,
 ):
@@ -153,8 +126,6 @@ def get(
         </path/to/alignments.sam>
     bam : str, optional
         </path/to/alignments.bam>
-    lengths : str, optional
-        </path/to/lengths.tsv>
     bed : str, optional
         </path/to/alignments.bed>
     cpus : int, optional
@@ -184,17 +155,14 @@ def get(
         tempdir = tempfile.mkdtemp(suffix=None, prefix="cov-alignments", dir=outdir)
         bed = bed if bed else os.path.join(tempdir, "alignment.bed")
         bam = bam if bam else os.path.join(tempdir, "alignment.bam")
-        lengths = lengths if lengths else os.path.join(tempdir, "lengths.tsv")
         sam = sam if sam else os.path.join(tempdir, "alignment.sam")
         db = os.path.join(tempdir, "alignment.db")
 
         def parse_bed(bed=bed, out=out):
             return bedtools.parse(bed, out)
 
-        def make_bed(lengths=lengths, fasta=fasta, bam=bam, bed=bed):
-            if not os.path.exists(lengths):
-                lengths = make_length_table(fasta, lengths)
-            bedtools.genomecov(bam, lengths, bed)
+        def make_bed(bam=bam, bed=bed):
+            bedtools.genomecov(bam, bed)
 
         def sort_samfile(sam=sam, bam=bam, cpus=cpus):
             samtools.sort(sam, bam, cpus=cpus)
@@ -273,10 +241,6 @@ def main():
     )
     parser.add_argument("--sam", help="</path/to/alignments.sam>")
     parser.add_argument("--bam", help="</path/to/alignments.bam>")
-    parser.add_argument(
-        "--lengths",
-        help="Path to tab-delimited lengths table with columns of contig & length.",
-    )
     parser.add_argument("--bed", help="</path/to/alignments.bed>")
     parser.add_argument(
         "--cpus",
@@ -317,7 +281,6 @@ def main():
         rev_reads=args.rev_reads,
         sam=args.sam,
         bam=args.bam,
-        lengths=args.lengths,
         bed=args.bed,
         cpus=args.cpus,
         out=args.out,

autometa/common/external/bedtools.py

@@ -36,16 +36,14 @@
 logger = logging.getLogger(__name__)
 
 
-def genomecov(ibam, lengths, out, force=False):
+def genomecov(ibam, out, force=False):
     """Run bedtools genomecov with input `ibam` and `lengths` to retrieve
     metagenome coverages.
 
     Parameters
     ----------
     ibam : str
         </path/to/indexed/BAM/file.ibam>. Note: BAM *must* be sorted by position.
-    lengths : str
-        </path/to/genome/lengths.tsv> tab-delimited cols=[contig,length]
     out : str
         </path/to/alignment.bed>
         The bedtools genomecov output is a tab-delimited file with the following columns:
@@ -71,7 +69,7 @@ def genomecov(ibam, lengths, out, force=False):
         Why the exception is raised.
 
     """
-    cmd = f"bedtools genomecov -ibam {ibam} -g {lengths}"
+    cmd = f"bedtools genomecov -ibam {ibam}"
     if os.path.exists(out) and not force:
         logger.debug(f"{out} already exists. skipping...")
         return out
@@ -126,7 +124,8 @@ def parse(bed, out=None, force=False):
     if out and (not os.path.exists(out) or (os.path.exists(out) and force)):
         dff.to_csv(out, sep="\t", index=True, header=True)
         logger.debug(f"{out} written")
-    logger.debug(f"{os.path.basename(out)} shape: {dff.shape}")
+    msg = f"{os.path.basename(out)} shape: {dff.shape}" if out else f"shape: {dff.shape}"
+    logger.debug(msg)
     return dff[["coverage"]]
 
 
@@ -144,12 +143,6 @@ def main():
     parser.add_argument(
         "--ibam", metavar="filepath", help="Path to sorted alignment.bam", required=True
     )
-    parser.add_argument(
-        "--lengths",
-        metavar="filepath",
-        help="Path to genome lengths.tsv; tab-delimited cols=[contig,length]",
-        required=True,
-    )
     parser.add_argument(
         "--bed",
         metavar="filepath",
@@ -167,15 +160,13 @@ def main():
     )
     parser.add_argument(
         "--force-cov",
-        help="force overwrite `--coverage`",
+        help="force overwrite `--output`",
         action="store_true",
         default=False,
     )
     args = parser.parse_args()
-    bed = genomecov(
-        ibam=args.ibam, lengths=args.lengths, out=args.bed, force=args.force_bed
-    )
-    parse(bed=bed, out=args.coverage, force=args.force_cov)
+    bed = genomecov(ibam=args.ibam, out=args.bed, force=args.force_bed)
+    parse(bed=bed, out=args.output, force=args.force_cov)
 
 
 if __name__ == "__main__":

modules/local/parse_bed.nf

@@ -28,7 +28,7 @@ process PARSE_BED {
     script:
         def software = getSoftwareName(task.process)
         """
-        autometa-parse-bed \\
+        autometa-bedtools-genomecov \\
             --ibam $bam \\
             --lengths $lengths \\
             --bed $bed_out \\

setup.py

@@ -30,7 +30,7 @@ def read(fname):
             "autometa-length-filter = autometa.common.metagenome:main",
             "autometa-markers = autometa.common.markers:main",
             "autometa-orfs = autometa.common.external.prodigal:main",
-            "autometa-parse-bed = autometa.common.external.bedtools:main",
+            "autometa-bedtools-genomecov = autometa.common.external.bedtools:main",
             "autometa-hmmsearch-filter = autometa.common.external.hmmsearch:main",
             "autometa-taxonomy = autometa.taxonomy.vote:main",
             "autometa-taxonomy-lca = autometa.taxonomy.lca:main",

tests/unit_tests/test_coverage.py

@@ -189,7 +189,6 @@ def __init__(self):
             self.rev_reads = rev_reads
             self.sam = sam_alignment
             self.bam = bam_alignment
-            self.lengths = "lengths.tsv"
             self.bed = bed_alignment
             self.cpus = 2
             self.out = out