Merge branch 'rz_rnaseqmetrics_parameterize_endbias_length' of https://github.com/watchmaker-genomics/picard into rz_rnaseqmetrics_parameterize_endbias_length

Author: ramsey-watchmaker

Dockerfile

@@ -1,4 +1,4 @@
-FROM openjdk:8
+FROM adoptopenjdk:8
 MAINTAINER Broad Institute DSDE <[email protected]>
 
 ARG build_command=shadowJar

src/main/java/picard/fingerprint/ConvertHaplotypeDatabaseToVcf.java

@@ -0,0 +1,54 @@
+package picard.fingerprint;
+
+import htsjdk.samtools.util.IOUtil;
+import org.broadinstitute.barclay.argparser.Argument;
+import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
+import picard.PicardException;
+import picard.cmdline.CommandLineProgram;
+import picard.cmdline.StandardOptionDefinitions;
+import picard.cmdline.programgroups.DiagnosticsAndQCProgramGroup;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+
+@CommandLineProgramProperties(
+        summary = "Convert Haplotype database file to vcf" +
+                "<h3>Examples</h3>" +
+                "<pre>" +
+                "java -jar picard.jar ConvertHaplotypeDatabaseToVcf \\\n" +
+                "-I haplotype_database.txt \\\n" +
+                "-O haplotype_database.vcf.gz \\\n" +
+                "-R reference.fasta" +
+                "</pre>",
+        oneLineSummary = "Convert Haplotype database file to vcf",
+        programGroup = DiagnosticsAndQCProgramGroup.class
+)
+public class ConvertHaplotypeDatabaseToVcf extends CommandLineProgram {
+
+    @Argument(doc = "Haplotype database to be converted to VCF.", shortName = StandardOptionDefinitions.INPUT_SHORT_NAME)
+    public File INPUT;
+
+    @Argument(doc = "Where to write converted haplotype database VCF.", shortName = StandardOptionDefinitions.OUTPUT_SHORT_NAME)
+    public File OUTPUT;
+
+    @Override
+    protected boolean requiresReference() {
+        return true;
+    }
+
+    @Override
+    protected int doWork() {
+        IOUtil.assertFileIsReadable(INPUT);
+        IOUtil.assertFileIsWritable(OUTPUT);
+
+        final HaplotypeMap haplotypeMap = new HaplotypeMap(INPUT);
+
+        try {
+            haplotypeMap.writeAsVcf(OUTPUT, REFERENCE_SEQUENCE);
+        } catch (final FileNotFoundException ex) {
+            throw new PicardException("Problem writing haplotype map to " + OUTPUT, ex);
+        }
+
+        return 0;
+    }
+}

src/main/java/picard/fingerprint/CrosscheckFingerprints.java

@@ -289,6 +289,16 @@ public class CrosscheckFingerprints extends CommandLineProgram {
             "Should only be used with SECOND_INPUT. ", optional = true)
     public File SECOND_INPUT_SAMPLE_MAP;
 
+    @Argument(doc = "A tsv with two columns representing the individual with which each sample is associated.  The first column is the sample id, and the second " +
+            "column is the associated individual id.  Values in the first column must be unique.  If INPUT_SAMPLE_MAP or SECOND_INPUT_SAMPLE_MAP is also specified, " +
+            "then the values in the first column of this file should be the sample aliases specified in the second columns of INPUT_SAMPLE_MAP and SECOND_INPUT_SAMPLE_MAP, " +
+            "respectively.  When this input is specified, expectations for matches will be based on the equality or inequality of the individual ids associated with two " +
+            "samples, as opposed to the sample ids themselves.  Samples which are not listed in this file will have their sample id used as their individual id, for the " +
+            "purposes of match expectations.  This means that one sample id could be used as the individual id for another sample, but not included in the map itself, and " +
+            "these two samples would be considered to have come from the same individual.  Note that use of this parameter only affects labelling of matches and mismatches as " +
+            "EXPECTED or UNEXPECTED.  It has no affect on how data is grouped for crosschecking.", optional = true)
+    public File SAMPLE_INDIVIDUAL_MAP;
+
     @Argument(doc = "An argument that controls how crosschecking with both INPUT and SECOND_INPUT should occur. ")
     public Fingerprint.CrosscheckMode CROSSCHECK_MODE = CHECK_SAME_SAMPLE;
 
@@ -368,6 +378,7 @@ public class CrosscheckFingerprints extends CommandLineProgram {
     private double[][] crosscheckMatrix = null;
     private final List<String> lhsMatrixKeys = new ArrayList<>();
     private final List<String> rhsMatrixKeys = new ArrayList<>();
+    private Map<String, String> sampleIndividualMap;
 
     @Override
     protected String[] customCommandLineValidation() {
@@ -432,6 +443,9 @@ protected int doWork() {
         if (SECOND_INPUT_SAMPLE_MAP != null) {
             IOUtil.assertFileIsReadable(SECOND_INPUT_SAMPLE_MAP);
         }
+        if (SAMPLE_INDIVIDUAL_MAP != null) {
+            IOUtil.assertFileIsReadable(SAMPLE_INDIVIDUAL_MAP);
+        }
 
         final HaplotypeMap map = new HaplotypeMap(HAPLOTYPE_MAP);
         final FingerprintChecker checker = new FingerprintChecker(map);
@@ -782,12 +796,22 @@ private int crossCheckFingerprints(final Map<FingerprintIdDetails, Fingerprint>
         // use 1L to promote size() to a long and avoid possible overflow
         final long totalChecks = lhsFingerprintIdDetails.size() * ((long) rhsFingerprintIdDetails.size());
 
+        if (SAMPLE_INDIVIDUAL_MAP != null) {
+            final List<FingerprintIdDetails> allFingerprintIdDetails = new ArrayList<>(lhsFingerprintIdDetails);
+            allFingerprintIdDetails.addAll(rhsFingerprintIdDetails);
+            final Set<String> inputSamples = allFingerprintIdDetails.stream().map(id -> id.sample).collect(Collectors.toSet());
+
+            sampleIndividualMap = buildSampleIndividualsMap(SAMPLE_INDIVIDUAL_MAP, inputSamples);
+        }
+
         for (int row = 0; row < lhsFingerprintIdDetails.size(); row++) {
             final FingerprintIdDetails lhsId = lhsFingerprintIdDetails.get(row);
 
             for (int col = 0; col < rhsFingerprintIdDetails.size(); col++) {
                 final FingerprintIdDetails rhsId = rhsFingerprintIdDetails.get(col);
-                final boolean expectedToMatch = EXPECT_ALL_GROUPS_TO_MATCH || lhsId.sample.equals(rhsId.sample);
+                final String lhsMatchId = resolveIndividualIfPossible(lhsId.sample);
+                final String rhsMatchId = resolveIndividualIfPossible(rhsId.sample);
+                final boolean expectedToMatch = EXPECT_ALL_GROUPS_TO_MATCH || lhsMatchId.equals(rhsMatchId);
 
                 final MatchResults results = FingerprintChecker.calculateMatchResults(lhsFingerprints.get(lhsId), rhsFingerprints.get(rhsId),
                         GENOTYPING_ERROR_RATE, LOSS_OF_HET_RATE, false, CALCULATE_TUMOR_AWARE_RESULTS);
@@ -926,4 +950,44 @@ private FingerprintResult getMatchResults(final boolean expectedToMatch, final M
             }
         }
     }
+
+    private String resolveIndividualIfPossible(final String sampleID) {
+        if (sampleIndividualMap != null) {
+            return sampleIndividualMap.getOrDefault(sampleID, sampleID);
+        }
+
+        return sampleID;
+    }
+
+    private Map<String, String> buildSampleIndividualsMap(final File sampleIndividualMapFile, final Set<String> sampleIDsInput) {
+        if (SAMPLE_INDIVIDUAL_MAP != null) {
+            final Map<String, String> individualMap = getStringStringMap(sampleIndividualMapFile, "SAMPLE_INDIVIDUAL_MAP");
+
+            //check for allowed but unusual situations that may indicate a mistake
+            //are there any sample ids used as individual ids?
+            for (final String individualID : individualMap.values()) {
+                if (sampleIDsInput.contains(individualID)) {
+                    log.warn("Sample " + individualID + " is used as an individual ID in " + sampleIndividualMapFile);
+                }
+            }
+
+            //are there any sample ids in input not found in the map?
+            for (final String sampleId : sampleIDsInput) {
+                if (!individualMap.keySet().contains(sampleId)) {
+                    log.warn("Sample " + sampleId + " in input data not found in " + sampleIndividualMapFile);
+                }
+            }
+
+            //are there any sample ids listed in the map that are not found in input data?
+            for (final String sampleId : individualMap.keySet()) {
+                if (!sampleIDsInput.contains(sampleId)) {
+                    log.warn("Sample " + sampleId + " in " + sampleIndividualMapFile + " not found in input data");
+                }
+            }
+
+            return individualMap;
+        }
+
+        return null;
+    }
 }

src/test/java/picard/fingerprint/ConvertHaplotypeDatabaseToVcfTest.java

@@ -0,0 +1,89 @@
+package picard.fingerprint;
+
+import org.testng.Assert;
+import org.testng.annotations.Test;
+import picard.cmdline.CommandLineProgramTest;
+
+import java.io.File;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+
+public class ConvertHaplotypeDatabaseToVcfTest extends CommandLineProgramTest {
+    private final File TEST_DATA_DIR = new File("testdata/picard/fingerprint/");
+    private final File TEST_MAP =
+            new File(TEST_DATA_DIR, "haplotypeMap_small.txt");
+
+    private final File TEST_FASTA =
+            new File(TEST_DATA_DIR, "reference.fasta");
+
+    @Override
+    public String getCommandLineProgramName() {
+        return ConvertHaplotypeDatabaseToVcf.class.getSimpleName();
+    }
+
+    @Test
+    public void testConvertHaplotypeDatabaseToVcf() throws IOException {
+        final File outfile = getTempOutputFile("testConvertHaplotypeDatabase", ".vcf");
+        final String[] args = new String[]{
+                "I=" + TEST_MAP.getAbsolutePath(),
+                "R=" + TEST_FASTA.getAbsolutePath(),
+                "O=" + outfile.getAbsolutePath()
+        };
+
+        runPicardCommandLine(args);
+
+        final HaplotypeMap hapDBFromTxt = new HaplotypeMap(TEST_MAP);
+        final HaplotypeMap hapDBFromVcf = new HaplotypeMap(outfile);
+
+        assertHaplotypeMapsAreEquivalent(hapDBFromTxt, hapDBFromVcf);
+    }
+
+    private void assertHaplotypeMapsAreEquivalent(final HaplotypeMap hapMap1, final HaplotypeMap hapMap2) {
+        final List<HaplotypeBlock> hapBlocks1 = new ArrayList<>(hapMap1.getHaplotypes());
+        final List<HaplotypeBlock> hapBlocks2 = new ArrayList<>(hapMap2.getHaplotypes());
+
+        Assert.assertEquals(hapBlocks1.size(), hapBlocks2.size());
+
+        Collections.sort(hapBlocks1);
+        Collections.sort(hapBlocks2);
+
+        for (int i = 0; i < hapBlocks1.size(); i++) {
+            assertHaplotypeBlocksAreEquivalent(hapBlocks1.get(i), hapBlocks2.get(i));
+        }
+    }
+
+    private void assertHaplotypeBlocksAreEquivalent(final HaplotypeBlock hapBlock1, final HaplotypeBlock hapBlock2) {
+        //equals method only checks for equality of domain
+        Assert.assertEquals(hapBlock1, hapBlock2);
+
+        final List<Snp> snps1 = new ArrayList<>(hapBlock1.getSnps());
+        final List<Snp> snps2 = new ArrayList<>(hapBlock2.getSnps());
+
+        Assert.assertEquals(snps1.size(), snps2.size());
+
+        Collections.sort(snps1);
+        Collections.sort(snps2);
+
+        for (int i = 0; i < snps1.size(); i++) {
+            assertSnpsAreEquivalent(snps1.get(i), snps2.get(i));
+        }
+    }
+
+    private void assertSnpsAreEquivalent(final Snp snp1, final Snp snp2) {
+        //equals method only checks for equality of domain
+        Assert.assertEquals(snp1, snp2);
+
+        Assert.assertEquals(snp1.getName(), snp2.getName());
+
+        //allele order may have been swapped
+        if (snp1.getAllele1() == snp2.getAllele1() && snp1.getAllele2() == snp2.getAllele2()) {
+            Assert.assertEquals(snp1.getMaf(), snp2.getMaf(), 0.0005); //equal rounded to third decimal
+        } else if (snp1.getAllele1() == snp2.getAllele2() && snp1.getAllele2() == snp2.getAllele1()) {
+            Assert.assertEquals(snp1.getMaf(), 1 - snp2.getMaf(), 0.0005); //equal rounded to third decimal
+        } else {
+            throw new AssertionError("Snps " + snp1.getName() + " have different alleles");
+        }
+    }
+}

src/test/java/picard/fingerprint/CrosscheckFingerprintsTest.java

@@ -61,6 +61,9 @@ public class CrosscheckFingerprintsTest extends CommandLineProgramTest {
     private File NA12891_r1_shifted_bam, NA12891_r2_shifted_bam, NA12892_r1_shifted_bam, NA12892_r2_shifted_bam;
     private final File referenceForCrams = new File(TEST_DATA_DIR, "reference.shifted.for.crams.fasta");
 
+    private final File NA12891_NA12892_diff_individuals = new File(TEST_DATA_DIR, "NA12891_NA12892_different_individuals.txt");
+    private final File NA12891_NA12892_same_individual = new File(TEST_DATA_DIR, "NA12891_NA12892_same_individual.txt");
+
     private final int NA12891_r1_RGs = 27;
     private final int NA12891_r2_RGs = 26;
     private final int NA12892_r1_RGs = 25;
@@ -432,6 +435,40 @@ public void testCrossCheckSMs(final File file1, final File file2, final int expe
         Assert.assertEquals(matrixParser.columnLabelsList().size(), numberOfSamples + 1 );
     }
 
+    @DataProvider(name = "bamFilesIndividuals")
+    public Object[][] bamFilesIndividuals() {
+            return new Object[][] {
+                    {NA12891_r1, NA12892_r1, NA12891_NA12892_diff_individuals, 0, 2},
+                    {NA12891_r1, NA12892_r1, NA12891_NA12892_same_individual, 1, 2}, // same individual, so expects match
+                    {NA12891_r2, NA12891_named_NA12892_r1, NA12891_NA12892_diff_individuals, 1, 2}, // diff individuals, so do not expect match
+                    {NA12891_r2, NA12891_named_NA12892_r1, NA12891_NA12892_same_individual, 0, 2} //same individual, so expect match
+            };
+    }
+
+    @Test(dataProvider = "bamFilesIndividuals")
+    public void testCrossCheckIndividuals(final File file1, final File file2, final File individualsMap, final int expectedRetVal, final int numberOfSamples) throws IOException {
+        File metrics = File.createTempFile("Fingerprinting", "NA1291.SM.crosscheck_metrics");
+        metrics.deleteOnExit();
+
+        File matrix = File.createTempFile("Fingerprinting", "NA1291.SM.matrix");
+        matrix.deleteOnExit();
+
+        final String[] args = new String[]{
+                "INPUT=" + file1.getAbsolutePath(),
+                "INPUT=" + file2.getAbsolutePath(),
+                "OUTPUT=" + metrics.getAbsolutePath(),
+                "MATRIX_OUTPUT=" + matrix.getAbsolutePath(),
+                "HAPLOTYPE_MAP=" + HAPLOTYPE_MAP,
+                "LOD_THRESHOLD=" + -1.0,
+                "CROSSCHECK_BY=SAMPLE",
+                "SAMPLE_INDIVIDUAL_MAP=" + individualsMap.getAbsolutePath()
+        };
+        doTest(args, metrics, expectedRetVal, numberOfSamples * numberOfSamples  , CrosscheckMetric.DataType.SAMPLE);
+
+        TabbedTextFileWithHeaderParser matrixParser = new TabbedTextFileWithHeaderParser(matrix);
+        Assert.assertEquals(matrixParser.columnLabelsList().size(), numberOfSamples + 1 );
+    }
+
     @DataProvider(name = "checkSamplesData")
     public Iterator<Object[]> checkSamplesData() {
         List<Object[]> tests = new ArrayList<>();

testdata/picard/fingerprint/NA12891_NA12892_different_individuals.txt

@@ -0,0 +1,2 @@
+NA12891	Individual1
+NA12892	Individual2

testdata/picard/fingerprint/NA12891_NA12892_same_individual.txt

@@ -0,0 +1,2 @@
+NA12891	Individual1
+NA12892	Individual1

testdata/picard/fingerprint/haplotypeMap_small.txt

@@ -0,0 +1,22 @@
+@HD	VN:1.6	SO:unsorted
+@SQ	SN:chr1	LN:101	M5:bd01f7e11515bb6beda8f7257902aa67
+@SQ	SN:chr2	LN:101	M5:31c33e2155b3de5e2554b693c475b310
+@SQ	SN:chr3	LN:101	M5:631593c6dd2048ae88dcce2bd505d295
+@SQ	SN:chr4	LN:101	M5:c60cb92f1ee5b78053c92bdbfa19abf1
+@SQ	SN:chr5	LN:101	M5:07ebc213c7611db0eacbb1590c3e9bda
+@SQ	SN:chr6	LN:101	M5:7be2f5e7ee39e60a6c3b5b6a41178c6d
+@SQ	SN:chr7	LN:202	M5:93763aaf6a455871c7d7a7718bff9ccf
+@SQ	SN:chr8	LN:202	M5:d339678efce576d5546e88b49a487b63
+#CHROMOSOME	POSITION	NAME	MAJOR_ALLELE	MINOR_ALLELE	MAF	ANCHOR_SNP	PANELS
+chr1	75	rs1	C	A	0.159601		foo
+chr1	100	rs2	A	G	0.223794		foo
+chr2	75	rs3	C	G	0.437188		foo
+chr3	25	rs9	C	T	0.722731		foo
+chr3	50	rs7	T	C	0.740000	rs6	foo
+chr3	75	rs6	A	T	0.623026		foo
+chr3	100	rs8	A	T	0.260000	rs6	foo
+chr4	75	rs10	A	G	0.205686		foo
+chr5	75	rs11	C	G	0.558528		foo
+chr6	75	rs15	G	C	0.452579		foo
+chr7	75	rs16	A	T	0.662917		foo
+chr8	75	rs17	G	A	0.392887		foo