Document CRAM input and test CRAM/BAM coverage equivalence

cnvlib/commands.py

@@ -273,7 +273,7 @@ def first_paths(fnames: list[str]) -> dict[str, str]:
 
 
 def _cmd_batch(args: argparse.Namespace) -> None:
-    """Run the complete CNVkit pipeline on one or more BAM or bedGraph files."""
+    """Run the complete CNVkit pipeline on one or more BAM, CRAM, or bedGraph files."""
     logging.info("CNVkit %s", __version__)
     # Validate/restrict options, beyond what argparse mutual exclusion can do
     bad_args_msg = ""
@@ -426,7 +426,7 @@ def _cmd_batch(args: argparse.Namespace) -> None:
 P_batch.add_argument(
     "sample_fnames",
     nargs="*",
-    help="""Mapped sequence reads (.bam) or pre-computed per-base depth
+    help="""Mapped sequence reads (.bam or .cram) or pre-computed per-base depth
             (bedGraph .bed.gz with tabix index .tbi or .csi)""",
 )
 P_batch.add_argument(
@@ -491,10 +491,10 @@ def _cmd_batch(args: argparse.Namespace) -> None:
     "--normal",
     nargs="*",
     metavar="FILES",
-    help="""Normal samples (.bam or bedGraph .bed.gz) used to construct the pooled, paired,
-            or flat reference. If this option is used but no filenames are given, a "flat"
-            reference will be built. Otherwise, all filenames following this option will
-            be used.""",
+    help="""Normal samples (.bam, .cram, or bedGraph .bed.gz) used to construct the
+            pooled, paired, or flat reference. If this option is used but no filenames
+            are given, a "flat" reference will be built. Otherwise, all filenames
+            following this option will be used.""",
 )
 add_fasta(P_batch_newref)
 P_batch_newref.add_argument(
@@ -749,7 +749,7 @@ def _cmd_antitarget(args: argparse.Namespace) -> None:
 
 
 def _cmd_autobin(args: argparse.Namespace) -> None:
-    """Quickly calculate reasonable bin sizes from BAM read counts."""
+    """Quickly calculate reasonable bin sizes from BAM or CRAM read counts."""
     if args.method in ("hybrid", "amplicon") and not args.targets:
         raise RuntimeError(f"Sequencing method {args.method!r} requires targets (-t)")
     if args.method == "wgs":
@@ -830,7 +830,8 @@ def read_regions(bed_fname):
 P_autobin.add_argument(
     "bams",
     nargs="+",
-    help="""Sample BAM file(s) or bedGraph file(s) (.bed.gz) to test for target coverage""",
+    help="""Sample BAM or CRAM file(s), or bedGraph file(s) (.bed.gz), to test for
+            target coverage""",
 )
 add_fasta(P_autobin)
 P_autobin.add_argument(
@@ -925,7 +926,7 @@ def read_regions(bed_fname):
 
 
 def _cmd_coverage(args: argparse.Namespace) -> None:
-    """Calculate coverage in the given regions from BAM read depths or bedGraph."""
+    """Calculate coverage in the given regions from BAM/CRAM read depths or bedGraph."""
     pset = coverage.do_coverage(
         args.interval,
         args.bam_or_bg_file,
@@ -953,7 +954,7 @@ def _cmd_coverage(args: argparse.Namespace) -> None:
 P_coverage.add_argument(
     "bam_or_bg_file",
     metavar="input_file",
-    help="""Mapped sequence reads (.bam) or pre-computed per-base depth
+    help="""Mapped sequence reads (.bam or .cram) or pre-computed per-base depth
             (bedGraph .bed.gz with tabix index .tbi or .csi)""",
 )
 P_coverage.add_argument("interval", help="Intervals (.bed or .list)")

doc/pipeline.rst

@@ -422,7 +422,7 @@ command.
 ``coverage``
 ------------
 
-Calculate coverage in the given regions from BAM read depths.
+Calculate coverage in the given regions from BAM or CRAM read depths.
 
 By default, coverage is calculated via mean read depth from a pileup.
 Alternatively, the ``--count`` option scans each read individually with pysam
@@ -474,6 +474,12 @@ these reads when calculating read depth.
     physically removed. (Exception: for amplicon sequencing, do *not* mark
     duplicates at all -- see :doc:`nonhybrid`.)
 
+CNVkit does not read base quality scores when calculating coverage, so base
+quality score recalibration does not affect it. GATK ApplyBQSR rewrites the
+quality scores and, apart from GATK's default filter for malformed reads, passes
+every read through with its alignment and flags unchanged, so the MarkDuplicates
+output and the recalibrated output yield the same coverage.
+
 You will probably want to index the finished BAM file using `samtools
 <https://samtools.sourceforge.net/>`_ or SAMBAMBA.  But if you haven't done this
 beforehand, CNVkit will automatically do it for you.
@@ -511,6 +517,77 @@ beforehand, CNVkit will automatically do it for you.
     or manual indexing are required (equivalent to running ``samtools index
     -c`` yourself).
 
+CRAM input
+``````````
+
+CNVkit accepts CRAM files wherever it accepts BAM files, including the
+:ref:`batch`, :ref:`coverage`, and :ref:`autobin` commands. Most CRAM files store
+each read as its differences from an external reference genome, which is then
+needed to read the file. Supply the FASTA file that the CRAM files were written
+against with the ``-f``/``--fasta`` option::
+
+    cnvkit.py coverage Sample.cram targets.bed -f hg38.fa -o Sample.targetcoverage.cnn
+
+    cnvkit.py batch *Tumor.cram --normal *Normal.cram --targets my_targets.bed \
+        --fasta hg38.fa --output-reference my_reference.cnn --output-dir results/
+
+The :ref:`batch` command also accepts ``-f`` together with ``-r``/``--reference``,
+for applying an existing copy number reference to new CRAM files. As with BAM
+files, CNVkit builds a missing or outdated CRAM index (``.crai``) automatically.
+
+.. note::
+    **Without** ``-f``, **reading a CRAM file can take many hours or fail.**
+    CNVkit reads alignments through pysam and HTSlib. HTSlib looks for each
+    reference sequence, in order: by name in the FASTA file given with ``-f``;
+    by the MD5 checksum in the ``M5`` tag of the CRAM header's ``@SQ`` line,
+    first in the ``REF_CACHE`` directory and then in each location listed in
+    ``REF_PATH``; and by name in the local FASTA file named in the ``@SQ``
+    line's ``UR`` tag. If ``REF_PATH`` is set and none of these supplies the
+    sequence, HTSlib reports ``Unable to fetch reference`` and ``Failure to
+    decode slice``, and CNVkit stops with an error. If ``REF_PATH`` is not set,
+    the outcome depends on the HTSlib version:
+
+    - HTSlib 1.21 and earlier (bundled with pysam 0.23 and earlier) use the
+      European Bioinformatics Institute's reference server as the default
+      ``REF_PATH``, downloading each sequence and storing it in a local cache,
+      ``~/.cache/hts-ref`` unless ``REF_CACHE`` is set. Fetching every chromosome
+      of a large genome this way can take many hours, and is a common reason for
+      CRAM processing to run far slower than BAM processing of the same data.
+      HTSlib logs a warning beginning ``[W::cram_populate_ref] Creating
+      reference cache directory`` only when it creates the default cache
+      directory; later downloads into that directory are silent.
+    - Later HTSlib versions (bundled with pysam 0.24 and later) have no default
+      server, and fail with the errors above.
+
+    Passing the matching FASTA file with ``-f`` avoids both outcomes. Where that
+    is not possible, the ``seq_cache_populate.pl`` script distributed with
+    samtools builds a local directory of sequences named by their MD5
+    checksums, and prints the path template to assign to ``REF_CACHE``, such
+    as ``/data/ref_cache/%2s/%2s/%s``. Assign the same template to ``REF_PATH``,
+    so that HTSlib 1.21 and earlier do not fall back to the EBI server for a
+    sequence missing from the cache. HTSlib's `reference sequence
+    documentation <https://www.htslib.org/doc/reference_seqs.html>`_ describes
+    these variables in full.
+
+.. note::
+    **The FASTA file must contain the same sequences the CRAM files were written
+    against.** Where it differs, HTSlib reports an ``MD5 checksum reference
+    mismatch``, and CNVkit stops with an error rather than computing coverage
+    from misread alignments. A warning such as ``[W::sanitise_SQ_lines] Header
+    @SQ length mismatch for ref chr1, 249250621 vs 248956422`` means that the
+    FASTA file's sequence lengths disagree with the CRAM header, which typically
+    indicates a different genome assembly (here, a CRAM header from hg19 and a
+    FASTA file from hg38). Decoding can still succeed where the two sequences
+    happen to agree over the region being read, so treat the warning as a sign
+    that the wrong FASTA file was given, even when no error follows.
+
+A CRAM file and a BAM file holding the same alignments yield identical coverage,
+in both the default and ``--count`` calculations, and CNVkit's test suite checks
+this on a fixture that includes duplicate-marked reads. Duplicate flags are
+honored in CRAM input exactly as in BAM input (see the note on marking
+duplicates above), and since coverage does not depend on base quality scores,
+CRAM's lossy quality-score compression options do not affect it either.
+
 
 .. _reference:
 

test/test_preprocessing.py

@@ -231,43 +231,71 @@ def test_bedcov_depth_independent_of_process_count(self):
     def test_coverage_skips_marked_duplicates(self):
         """Coverage ignores reads flagged as duplicates in both the depth
         (bedcov) and by-count paths, so marking duplicates (e.g. Picard
-        MarkDuplicates) is equivalent to physically removing them (#689)."""
+        MarkDuplicates) is equivalent to physically removing them (#689).
+
+        The same holds for CRAM input (#1105). Base qualities play no part
+        either: half of the surviving reads carry Phred 2 at every base, so
+        a recalibration that rewrites only base qualities, such as GATK
+        ApplyBQSR, cannot change the result.
+        """
         seqlen = 50
-        n_total, n_dup = 10, 6  # 4 reads survive the duplicate filter
+        n_total, n_dup = 10, 6
+        n_kept = n_total - n_dup
         header = {
             "HD": {"VN": "1.6", "SO": "coordinate"},
             "SQ": [{"SN": "chr1", "LN": 1000}],
         }
         tmpdir = tempfile.mkdtemp()
-        bam = os.path.join(tmpdir, "dup.bam")
-        with pysam.AlignmentFile(bam, "wb", header=header) as out:
-            for i in range(n_total):
-                a = pysam.AlignedSegment()
-                a.query_name = f"r{i}"
-                a.query_sequence = "A" * seqlen
-                a.flag = 0
-                a.reference_id = 0
-                a.reference_start = 100
-                a.mapping_quality = 60
-                a.cigarstring = f"{seqlen}M"
-                a.query_qualities = pysam.qualitystring_to_array("I" * seqlen)
-                a.is_duplicate = i >= (n_total - n_dup)
-                out.write(a)
-        pysam.index(bam)
-        bed = os.path.join(tmpdir, "region.bed")
-        with open(bed, "w") as fh:
-            fh.write("chr1\t90\t160\tregionA\n")
-        n_kept = n_total - n_dup
         try:
-            # Depth path (default): basecount counts only non-duplicate reads
-            table = coverage.bedcov(bed, bam, 0)
-            self.assertEqual(float(table["basecount"].iloc[0]), n_kept * seqlen)
-            # By-count path: same surviving read count
-            bamfile = pysam.AlignmentFile(bam, "rb")
-            count, _row = coverage.region_depth_count(
-                bamfile, "chr1", 90, 160, "regionA", 0
-            )
-            self.assertEqual(count, n_kept)
+            # The CRAM is written against a FASTA that is then moved, so its
+            # @SQ UR tag names a missing file. With REF_PATH and REF_CACHE
+            # pointing at a nonexistent directory, HTSlib has no other source (and
+            # does not fall back to the EBI server), so decoding depends on the
+            # FASTA passed to CNVkit.
+            written_fasta = os.path.join(tmpdir, "written.fa")
+            with open(written_fasta, "w") as fh:
+                fh.write(">chr1\n" + "A" * 1000 + "\n")
+            alns = {}
+            for ext, mode in (("bam", "wb"), ("cram", "wc")):
+                aln = alns[ext] = os.path.join(tmpdir, f"dup.{ext}")
+                with pysam.AlignmentFile(
+                    aln, mode, header=header, reference_filename=written_fasta
+                ) as out:
+                    for i in range(n_total):
+                        a = pysam.AlignedSegment()
+                        a.query_name = f"r{i}"
+                        a.query_sequence = "A" * seqlen
+                        a.flag = 0
+                        a.reference_id = 0
+                        a.reference_start = 100
+                        a.mapping_quality = 60
+                        a.cigarstring = f"{seqlen}M"
+                        a.query_qualities = pysam.qualitystring_to_array(
+                            ("I" if i % 2 else "#") * seqlen
+                        )
+                        a.is_duplicate = i >= n_kept
+                        out.write(a)
+                pysam.index(aln)
+            fasta = os.path.join(tmpdir, "ref.fa")
+            os.rename(written_fasta, fasta)
+            bed = os.path.join(tmpdir, "region.bed")
+            with open(bed, "w") as fh:
+                fh.write("chr1\t90\t160\tregionA\n")
+            no_cache = os.path.join(tmpdir, "no-cache", "%s")
+            ref_env = {"REF_PATH": no_cache, "REF_CACHE": no_cache}
+            for ext, aln in alns.items():
+                with self.subTest(fmt=ext), mock.patch.dict(os.environ, ref_env):
+                    # Depth path (default): only non-duplicate reads contribute
+                    table = coverage.bedcov(bed, aln, 0, fasta)
+                    self.assertEqual(float(table["basecount"].iloc[0]), n_kept * seqlen)
+                    # By-count path: same surviving read count
+                    with pysam.AlignmentFile(
+                        aln, "rb", reference_filename=fasta
+                    ) as alnfile:
+                        count, _row = coverage.region_depth_count(
+                            alnfile, "chr1", 90, 160, "regionA", 0
+                        )
+                    self.assertEqual(count, n_kept)
         finally:
             shutil.rmtree(tmpdir)
 
@@ -554,18 +582,31 @@ def test_coverage(self):
             self.assertGreater(cna.log2.nunique(), 1)
 
     def test_coverage_cram(self):
-        """The 'coverage' command with CRAM input."""
+        """CRAM input yields the same coverage as a BAM of the same alignments.
+
+        CRAM decoding reproduces the alignments exactly, so the two formats
+        must agree exactly in both the depth and by-count paths. The fixture
+        holds duplicate-marked reads within the chrM target, so the comparison
+        also requires both formats to apply the duplicate flag alike (#1105).
+
+        The fixture's @SQ UR tags name an absolute path on the machine that
+        wrote it, so where that path exists, this test cannot show that the
+        FASTA argument is used; test_coverage_skips_marked_duplicates does.
+        """
         bed = "formats/my-targets.bed"
         cram = "formats/na12878-chrM-Y-trunc.cram"
         fasta = "formats/chrM-Y-trunc.hg19.fa"
         bam = "formats/na12878-chrM-Y-trunc.bam"
-        # Verify CRAM coverage matches BAM coverage
-        cna_cram = commands.do_coverage(bed, cram, fasta=fasta)
-        cna_bam = commands.do_coverage(bed, bam)
-        self.assertEqual(len(cna_cram), len(cna_bam))
-        np.testing.assert_array_almost_equal(
-            cna_cram["depth"].values, cna_bam["depth"].values, decimal=1
-        )
+        for by_count in (False, True):
+            with self.subTest(by_count=by_count):
+                cna_cram = commands.do_coverage(
+                    bed, cram, by_count=by_count, fasta=fasta
+                )
+                cna_bam = commands.do_coverage(bed, bam, by_count=by_count)
+                self.assertEqual(len(cna_cram), 4)  # One row per BED region
+                np.testing.assert_array_equal(
+                    cna_cram["depth"].values, cna_bam["depth"].values
+                )
 
     def test_coverage_bedgraph(self):
         """The 'coverage' command with bedGraph input."""