CRAM 3.1 read support

.github/workflows/gatk-tests.yml

@@ -118,6 +118,9 @@ jobs:
         if: needs.check-secrets.outputs.google-credentials == 'true'
         uses: google-github-actions/setup-gcloud@v2
 
+      - name: 'Install Samtools'
+        run: scripts/install-samtools.sh
+
       - name: pull lfs files
         run: git lfs pull
 
@@ -190,6 +193,9 @@ jobs:
         if: needs.check-secrets.outputs.google-credentials == 'true'
         uses: google-github-actions/setup-gcloud@v2
 
+      - name: 'Install Samtools'
+        run: scripts/install-samtools.sh
+
       - name: build test jars
         run: ./gradlew clean shadowTestClassJar shadowTestJar
 

scripts/install-samtools.sh

@@ -0,0 +1,12 @@
+#!/bin/sh
+set -ex
+#ubuntu specific
+sudo apt-get update
+sudo apt-get upgrade
+sudo apt-get install -y libncurses-dev libbz2-dev liblzma-dev
+
+#install from the github tar
+export SAMTOOLS_VERSION=1.21
+wget https://github.com/samtools/samtools/releases/download/${SAMTOOLS_VERSION}/samtools-${SAMTOOLS_VERSION}.tar.bz2
+tar -xjvf samtools-${SAMTOOLS_VERSION}.tar.bz2
+cd samtools-${SAMTOOLS_VERSION} && ./configure --prefix=/usr && make && sudo make install

src/test/java/org/broadinstitute/hellbender/GATKBaseTest.java

@@ -75,6 +75,9 @@ public abstract class GATKBaseTest extends BaseTest {
     public static final String NA12878_20_21_WGS_bam = largeFileTestDir + "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.bam";
     public static final String NA12878_20_21_WGS_mmp2_bam = largeFileTestDir + "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.mmp2.bam";
     public static final String NA12878_20_21_WGS_cram = largeFileTestDir + "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.0.samtools.cram";
+    // created using samtools archive mode ("--output-fmt cram,version=3.1,archive"); contains slices that exercise all
+    // of the cram 3.1 codecs (rANSnx16, adaptive arithmetic, fqzcomp, and name tokenization)
+    public static final String NA12878_20_21_WGS_cram_31 = largeFileTestDir + "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.1.samtools.archive.cram";
 
     public static final String NA12878_20_21_covered_regions = publicTestDir + "wgs_calling_regions.v1.chr20_chr21.interval_list";
 

src/test/java/org/broadinstitute/hellbender/PrintFileDiagnosticsIntegrationTest.java

@@ -1,6 +1,9 @@
 package org.broadinstitute.hellbender;
 
+import htsjdk.beta.plugin.IOUtils;
+import htsjdk.io.IOPath;
 import org.apache.commons.lang3.tuple.Pair;
+import org.broadinstitute.hellbender.engine.GATKPath;
 import org.broadinstitute.hellbender.testutils.ArgumentsBuilder;
 import org.broadinstitute.hellbender.testutils.IntegrationTestSpec;
 import org.testng.annotations.DataProvider;
@@ -14,11 +17,12 @@ public class PrintFileDiagnosticsIntegrationTest extends CommandLineProgramTest
 
     @DataProvider(name = "fileDiagnosticsTestCases")
     public Object[][] getFileDiagnosticsTestCases() {
+        // the pathnames used by the diagnostics tool wind up embedded in the diagnostics output file, so for these
+        // tests use just a relative pathname as input (instead of the named constants, i.e., NA12878_20_21_WGS_cram,
+        // which are full path names) in order to avoid test failures caused by the full pathname varying in
+        // different environments, i.e. in CI
         return new Object[][]{
                 {
-                        //this pathname is embedded in the diagnostics output file, so we use a relative pathname
-                        // instead of the named constant NA12878_20_21_WGS_cram in order to avoid test failures
-                        // caused by the full pathname varying in different environments
                         "src/test/resources/large/CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.0.samtools.cram",
                         List.of(Pair.of("count-limit", "10")),
                         "src/test/resources/filediagnostics/CEUTrio.HiSeq.WGS.b37.NA12878.20.21.txt"
@@ -33,6 +37,12 @@ public Object[][] getFileDiagnosticsTestCases() {
                         null,
                         "src/test/resources/filediagnostics/cram_with_crai_index.cram.crai.txt"
                 },
+                {
+                    // cram file that uses all the new 3.1 codecs (fqzcomp, name tok, ransNx16, and adaptive arithmetic)
+                        "src/test/resources/large/CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.1.samtools.archive.cram",
+                        List.of(Pair.of("count-limit", "20")),
+                        "src/test/resources/filediagnostics/CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.1.samtools.archive.cram.txt"
+                }
         };
     }
 
@@ -41,15 +51,22 @@ public void testFileDiagnostics(
             final String inputPath,
             final List<Pair<String, String>> extraArgs,
             final String expectedOutputPath) throws IOException {
-        final File outFile = createTempFile("testFileDiagnostics", ".txt");
-        ArgumentsBuilder argBuilder = new ArgumentsBuilder();
-        argBuilder.addInput(inputPath);
-        argBuilder.addOutput(outFile);
+        final IOPath outFile = IOUtils.createTempPath("testFileDiagnostics", ".txt");
+        runFileDiagnosticsTool(new GATKPath(inputPath), extraArgs, outFile);
+        IntegrationTestSpec.assertEqualTextFiles(outFile.toPath().toFile(), new File(expectedOutputPath));
+    }
+
+    private void runFileDiagnosticsTool(
+            final IOPath inputPath,
+            final List<Pair<String, String>> extraArgs,
+            final IOPath outputPath) {
+        final ArgumentsBuilder argBuilder = new ArgumentsBuilder();
+        argBuilder.addInput(inputPath.getRawInputString());
+        argBuilder.addOutput(outputPath.getRawInputString());
         if (extraArgs != null) {
             extraArgs.forEach(argPair -> argBuilder.add(argPair.getKey(), argPair.getValue()));
         }
         runCommandLine(argBuilder.getArgsList());
-
-        IntegrationTestSpec.assertEqualTextFiles(outFile, new File(expectedOutputPath));
     }
+
 }

src/test/java/org/broadinstitute/hellbender/engine/CRAMSupportIntegrationTest.java

@@ -43,13 +43,13 @@ public Object[][] getRoundTripCRAMTests() {
                 // read equality; at least some of which are because they are unmapped/unplaced, but have cigar
                 // strings that both samtools and htsjdk drop when roundtripping
                 {NA12878_20_21_WGS_bam, b37_reference_20_21, true, false},
+                // roundtrip a v3.0 file
                 // this cram is the result of converting the above bam to cram using samtools; once the file is
                 // converted, we can use full read equality when roundtripping through cram, so we don't need to
                 // be lenient
                 {NA12878_20_21_WGS_cram, b37_reference_20_21, false, false},
-                // roundtrip a v2.1 file
-                { largeFileTestDir + "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.0.samtools.cram",
-                        b37_reference_20_21, false, false },
+                // roundtrip a v3.1 file
+                {NA12878_20_21_WGS_cram_31, b37_reference_20_21, false, false },
         };
     }
 

src/test/java/org/broadinstitute/hellbender/testutils/SamtoolsTestUtilsTest.java

@@ -0,0 +1,80 @@
+package org.broadinstitute.hellbender.testutils;
+
+import htsjdk.io.HtsPath;
+import htsjdk.io.IOPath;
+import htsjdk.samtools.SAMRecord;
+import htsjdk.samtools.SamReader;
+import htsjdk.samtools.SamReaderFactory;
+import htsjdk.samtools.ValidationStringency;
+import htsjdk.samtools.util.CloseableIterator;
+import htsjdk.samtools.util.ProcessExecutor;
+import org.broadinstitute.hellbender.GATKBaseTest;
+import org.testng.Assert;
+import org.testng.SkipException;
+import org.testng.annotations.Test;
+
+import java.io.File;
+import java.io.IOException;
+
+public class SamtoolsTestUtilsTest extends GATKBaseTest {
+    private static final File TEST_DATA_DIR = new File("src/test/resources/org/broadinstitute/hellbender/tools/");
+
+    @Test
+    public void testSamtoolsIsAvailable() {
+        Assert.assertTrue(SamtoolsTestUtils.isSamtoolsAvailable());
+    }
+
+    @Test
+    public void testSamtoolsVersion() {
+        if (isGATKDockerContainer()) {
+            // since this test confirms that we're running a recent (specific) version of samtools for cram 3.1
+            // validation, skip it when running on the GATK Docker container, since the docker doesn't have the same
+            // recent samtools version that we use to verify cram 3.1 functionality
+            throw new SkipException("Samtools not available in GATK Docker container");
+        }
+        if (!SamtoolsTestUtils.isSamtoolsAvailable()) {
+            throw new SkipException("Samtools not available on local device");
+        }
+        // If this test runs, but fails because version validation fails, then the local samtools version is
+        // not the one expected by the htsjdk tests
+        final ProcessExecutor.ExitStatusAndOutput processStatus = SamtoolsTestUtils.executeSamToolsCommand("--version");
+        Assert.assertTrue(processStatus.stdout.contains(SamtoolsTestUtils.expectedSamtoolsVersion));
+    }
+
+    @Test(expectedExceptions = RuntimeException.class)
+    public void testSamtoolsPresentButCommandFails() {
+        if (!SamtoolsTestUtils.isSamtoolsAvailable()) {
+            throw new SkipException("Samtools not available on local device");
+        }
+        SamtoolsTestUtils.executeSamToolsCommand("--notASamtoolsCommand");
+    }
+
+    @Test
+    public void testCRAMConversion()throws IOException {
+        if (!SamtoolsTestUtils.isSamtoolsAvailable()) {
+            throw new SkipException("Samtools not available on local device");
+        }
+
+        // Validates CRAM 3.1 conversion.
+        final File sourceFile = new File(TEST_DATA_DIR, "print_reads.cram");
+        final File cramReference = new File(TEST_DATA_DIR, "print_reads.fasta");
+        // This also validates that any extra command line arguments are passed through to samtools by requesting
+        // that NM/MD values are synthesized in the output file (which is required for the output records to match).
+        final IOPath tempSamtoolsPath = SamtoolsTestUtils.convertToCRAM(
+                new HtsPath(sourceFile.getAbsolutePath()),
+                new HtsPath(cramReference.getAbsolutePath()),
+                "--output-fmt cram,version=3.0,fast");
+        final SamReaderFactory factory = SamReaderFactory.makeDefault()
+                .validationStringency(ValidationStringency.LENIENT)
+                .referenceSequence(cramReference);
+        try (final SamReader originalReader = factory.open(sourceFile);
+             final SamReader samtoolsCopyReader = factory.open(tempSamtoolsPath.toPath());
+             final CloseableIterator<SAMRecord> originalIt = originalReader.iterator();
+             final CloseableIterator<SAMRecord> samtoolsIt = samtoolsCopyReader.iterator()) {
+            while (originalIt.hasNext() && samtoolsIt.hasNext()) {
+                Assert.assertEquals(originalIt.next(), samtoolsIt.next());
+            }
+            Assert.assertEquals(samtoolsIt.hasNext(), originalIt.hasNext());
+        }
+    }
+}