Java Code Examples for htsjdk.samtools.SAMReadGroupRecord#setPlatform()

/** Creates a simple header with the values provided on the command line. */
public SAMFileHeader createSamFileHeader() {
    final SAMReadGroupRecord rgroup = new SAMReadGroupRecord(this.READ_GROUP_NAME);
    rgroup.setSample(this.SAMPLE_NAME);
    if (this.LIBRARY_NAME != null) rgroup.setLibrary(this.LIBRARY_NAME);
    if (this.PLATFORM != null) rgroup.setPlatform(this.PLATFORM);
    if (this.PLATFORM_UNIT != null) rgroup.setPlatformUnit(this.PLATFORM_UNIT);
    if (this.SEQUENCING_CENTER != null) rgroup.setSequencingCenter(SEQUENCING_CENTER);
    if (this.PREDICTED_INSERT_SIZE != null) rgroup.setPredictedMedianInsertSize(PREDICTED_INSERT_SIZE);
    if (this.DESCRIPTION != null) rgroup.setDescription(this.DESCRIPTION);
    if (this.RUN_DATE != null) rgroup.setRunDate(this.RUN_DATE);
    if (this.PLATFORM_MODEL != null) rgroup.setPlatformModel(this.PLATFORM_MODEL);
    if (this.PROGRAM_GROUP != null) rgroup.setProgramGroup(this.PROGRAM_GROUP);

    final SAMFileHeader header = new SAMFileHeader();
    header.addReadGroup(rgroup);

    for (final String comment : COMMENT) {
        header.addComment(comment);
    }

    header.setSortOrder(this.SORT_ORDER);
    return header ;
}
 

/**
 * Section of code shared between the two recalibration walkers which uses the command line arguments to adjust attributes of the read such as quals or platform string
 *
 * @param read The read to adjust
 * @param RAC  The list of shared command line arguments
 */
public static void parsePlatformForRead(final GATKRead read, final SAMFileHeader header, final RecalibrationArgumentCollection RAC) {
    final SAMReadGroupRecord readGroup = ReadUtils.getSAMReadGroupRecord(read, header);

    if (RAC.FORCE_PLATFORM != null && (readGroup.getPlatform() == null || !readGroup.getPlatform().equals(RAC.FORCE_PLATFORM))) {
        readGroup.setPlatform(RAC.FORCE_PLATFORM);
    }

    if (readGroup.getPlatform() == null) {
        if (RAC.DEFAULT_PLATFORM != null) {
            if (!warnUserNullPlatform) {
                Utils.warnUser("The input .bam file contains reads with no platform information. " +
                        "Defaulting to platform = " + RAC.DEFAULT_PLATFORM + ". " +
                        "First observed at read with name = " + read.getName());
                warnUserNullPlatform = true;
            }
            readGroup.setPlatform(RAC.DEFAULT_PLATFORM);
        }
        else {
            throw new UserException.MalformedRead(read, "The input .bam file contains reads with no platform information. First observed at read with name = " + read.getName());
        }
    }
}
 

protected SAMReadGroupRecord createReadGroupRecord(
        String RGID, String RGLB, String RGPL, 
        String RGPU, String RGSM, String RGCN, 
        String RGDS, Iso8601Date RGDT, Integer RGPI) {
    SAMReadGroupRecord rg = new SAMReadGroupRecord(RGID);
    rg.setLibrary(RGLB);
    rg.setPlatform(RGPL);
    rg.setSample(RGSM);
    rg.setPlatformUnit(RGPU);
    if(RGCN != null)
        rg.setSequencingCenter(RGCN);
    if(RGDS != null)
        rg.setDescription(RGDS);
    if(RGDT != null)
        rg.setRunDate(RGDT);
    if(RGPI != null)
        rg.setPredictedMedianInsertSize(RGPI);
    return rg;
}
 

public List<SAMReadGroupRecord> getReadGroups(final List<String> cellBarcodes) {
	List<SAMReadGroupRecord> g = new ArrayList<>(cellBarcodes.size());
	for (String id: cellBarcodes) {
		SAMReadGroupRecord rg = new SAMReadGroupRecord(id);
		rg.setLibrary(id);
	    rg.setPlatform(id);
	    rg.setSample(id);
	    rg.setPlatformUnit(id);
		g.add(rg);
	}
	return (g);
}
 

public void setupTest1(final int ID, final String readGroupId, final SAMReadGroupRecord readGroupRecord, final String sample,
                  final String library, final SAMFileHeader header, final SAMRecordSetBuilder setBuilder)
        throws IOException {

    final String separator = ":";
    final int contig1 = 0;
    final int contig2 = 1;
    readGroupRecord.setSample(sample);
    readGroupRecord.setPlatform(platform);
    readGroupRecord.setLibrary(library);
    readGroupRecord.setPlatformUnit(readGroupId);
    header.addReadGroup(readGroupRecord);
    setBuilder.setReadGroup(readGroupRecord);
    setBuilder.setUseNmFlag(true);

    setBuilder.setHeader(header);

    final int max = 800;
    final int min = 1;
    final Random rg = new Random(5);

    //add records that align to chrM and O but not N
    for (int i = 0; i < NUM_READS; i++) {
        final int start = rg.nextInt(max) + min;
        final String newReadName = READ_NAME + separator + ID + separator + i;

        if (i != NUM_READS - 1) {
            setBuilder.addPair(newReadName, contig1, start + ID, start + ID + LENGTH);
        } else {
            setBuilder.addPair(newReadName, contig2, start + ID, start + ID + LENGTH);
        }
    }
}
 

public void setupTest2(final int ID, final String readGroupId, final SAMReadGroupRecord readGroupRecord, final String sample,
                       final String library, final SAMFileHeader header, final SAMRecordSetBuilder setBuilder)
        throws IOException {

    final String separator = ":";
    final int contig1 = 0;
    final int contig2 = 1;
    final int contig3 = 2;
    readGroupRecord.setSample(sample);
    readGroupRecord.setPlatform(platform);
    readGroupRecord.setLibrary(library);
    readGroupRecord.setPlatformUnit(readGroupId);
    setBuilder.setReadGroup(readGroupRecord);
    setBuilder.setUseNmFlag(true);

    setBuilder.setHeader(header);

    final int max = 800;
    final int min = 1;
    final Random rg = new Random(5);

    //add records that align to all 3 chr in reference file
    for (int i = 0; i < NUM_READS; i++) {
        final int start = rg.nextInt(max) + min;
        final String newReadName = READ_NAME + separator + ID + separator + i;

        if (i<=NUM_READS/3) {
            setBuilder.addPair(newReadName, contig1, start + ID, start + ID + LENGTH);
        } else if (i< (NUM_READS - (NUM_READS/3))) {
            setBuilder.addPair(newReadName, contig2, start + ID, start + ID + LENGTH);
        } else {
            setBuilder.addPair(newReadName, contig3, start + ID, start + ID + LENGTH);
        }
    }
}
 

void setup(final int numReads,
           final String readName,
           final int ID,
           final String readGroupId,
           final SAMReadGroupRecord readGroupRecord,
           final String sample,
           final String library,
           final SAMFileHeader header,
           final SAMRecordSetBuilder setBuilder) throws IOException {
    final String separator = ":";
    readGroupRecord.setSample(sample);
    readGroupRecord.setPlatform(platform);
    readGroupRecord.setLibrary(library);
    readGroupRecord.setPlatformUnit(readGroupId);
    header.addReadGroup(readGroupRecord);
    setBuilder.setReadGroup(readGroupRecord);
    setBuilder.setUseNmFlag(true);
    setBuilder.setHeader(header);

    final int max = 15000;
    final int min = 1;
    final Random rg = new Random(5);

    for (int i = 0; i < numReads; i++) {
        final int start = rg.nextInt(max) + min;
        final String newReadName = readName + separator + ID + separator + i;
        setBuilder.addPair(newReadName, 0, start + ID, start + ID + 99);
    }
}
 

@BeforeClass
public void init() {
    RAC = new RecalibrationArgumentCollection();
    covariate = new CycleCovariate(RAC);
    illuminaReadGroup = new SAMReadGroupRecord("MY.ID");
    illuminaReadGroup.setPlatform("illumina");
}
 

@Test(enabled = true, dataProvider = "LIBS_NotHoldingTooManyReads")
public void testLIBS_NotHoldingTooManyReads(final int nReadsPerLocus, final int downsampleTo, final int payloadInBytes) {
    logger.warn(String.format("testLIBS_NotHoldingTooManyReads %d %d %d", nReadsPerLocus, downsampleTo, payloadInBytes));
    final int readLength = 10;

    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeader(1, 1, 100000);
    final int nSamples = 1;
    final List<String> samples = new ArrayList<>(nSamples);
    for ( int i = 0; i < nSamples; i++ ) {
        final SAMReadGroupRecord rg = new SAMReadGroupRecord("rg" + i);
        final String sample = "sample" + i;
        samples.add(sample);
        rg.setSample(sample);
        rg.setPlatform(NGSPlatform.ILLUMINA.getDefaultPlatform());
        header.addReadGroup(rg);
    }

    final boolean downsample = downsampleTo != -1;
    final DownsamplingMethod downsampler = downsample
            ? new DownsamplingMethod(DownsampleType.BY_SAMPLE, downsampleTo, null)
            : new DownsamplingMethod(DownsampleType.NONE, null, null);

    final WeakReadTrackingIterator iterator = new WeakReadTrackingIterator(nReadsPerLocus, readLength, payloadInBytes, header);

    final LocusIteratorByState li;
    li = new LocusIteratorByState(
            iterator,
            downsampler,
            false,
            samples,
            header,
            true
    );

    while ( li.hasNext() ) {
        final AlignmentContext next = li.next();
        Assert.assertTrue(next.getBasePileup().size() <= downsampleTo, "Too many elements in pileup " + next);
        // TODO -- assert that there are <= X reads in memory after GC for some X
    }
}
 

/**
 * A unit test that creates an artificial read for testing some code that uses reads
 */
@Test(dataProvider = "TestMappings")
public void testPLFromReadWithRG(final String plField, final NGSPlatform expected) {
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeader(sequenceDictionary);
    final String rgID = "ID";
    final SAMReadGroupRecord rg = new SAMReadGroupRecord(rgID);
    if ( plField != null )
        rg.setPlatform(plField);
    header.addReadGroup(rg);
    final GATKRead read = ArtificialReadUtils.createArtificialRead(header, "myRead", 0, 1, 10);
    read.setAttribute("RG", rgID);
    Assert.assertEquals(NGSPlatform.fromRead(read, header), expected);
}
 

@Override
protected int mainExec(OutputStream out, LogStream initLog) throws IOException {
  final PrintStream outStream = new PrintStream(out);
  try {
    final File output = (File) mFlags.getValue(CommonFlags.OUTPUT_FLAG);
    try (PrintStream summaryStream = new PrintStream(new FileOutputStream(new File(output, CommonFlags.SUMMARY_FILE)))) {
      final Collection<File> inputFiles = CommonFlags.getFileList(mFlags, CommonFlags.INPUT_LIST_FLAG, null, false);
      final Integer maximumNs = (Integer) mFlags.getValue(MAXIMUM_NS);

      final boolean useQuality = !mFlags.isSet(NO_QUALITY);
      final SAMReadGroupRecord samReadGroupRecord;
      if (mFlags.isSet(SamCommandHelper.SAM_RG)) {
        samReadGroupRecord = SamCommandHelper.validateAndCreateSamRG((String) mFlags.getValue(SamCommandHelper.SAM_RG), SamCommandHelper.ReadGroupStrictness.REQUIRED);
        final String platform = samReadGroupRecord.getPlatform();
        final MachineType mt = MachineType.COMPLETE_GENOMICS_2; // Default to new read structure
        if (!mt.compatiblePlatform(platform)) {
          if (platform == null || platform.length() == 0) {
            Diagnostic.warning("Read group platform not set, defaulting to \"" + mt.platform() + "\"");
            samReadGroupRecord.setPlatform(mt.platform());
          }
        }
      } else {
        samReadGroupRecord = null;
      }
      final ByteArrayOutputStream bos = new ByteArrayOutputStream();
      try (final PrintStream ps = new PrintStream(bos)) {
        performPreread(inputFiles, output, maximumNs, useQuality, ps, (Boolean) mFlags.getValue(COMPRESS_FLAG), mFlags.isSet(KEEP_NAMES), samReadGroupRecord);
        ps.flush();
        outStream.print(bos.toString());
        summaryStream.print(bos.toString());
      } catch (final IOException e) {
        if (output.getUsableSpace() == 0) {
          throw new NoTalkbackSlimException(e, ErrorType.DISK_SPACE, output.getPath());
        } else {
          throw e;
        }
      }
      return 0;
    }
  } finally {
    outStream.flush();
  }
}
 

protected int doWork() {
    IOUtil.assertInputIsValid(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);

    final SamReader in = SamReaderFactory.makeDefault()
        .referenceSequence(REFERENCE_SEQUENCE)
        .open(SamInputResource.of(INPUT));

    // create the read-group we'll be using
    final SAMReadGroupRecord rg = new SAMReadGroupRecord(RGID);
    rg.setLibrary(RGLB);
    rg.setPlatform(RGPL);
    rg.setSample(RGSM);
    rg.setPlatformUnit(RGPU);
    if (RGCN != null) rg.setSequencingCenter(RGCN);
    if (RGDS != null) rg.setDescription(RGDS);
    if (RGDT != null) rg.setRunDate(RGDT);
    if (RGPI != null) rg.setPredictedMedianInsertSize(RGPI);
    if (RGPG != null) rg.setProgramGroup(RGPG);
    if (RGPM != null) rg.setPlatformModel(RGPM);
    if (RGKS != null) rg.setKeySequence(RGKS);
    if (RGFO != null) rg.setFlowOrder(RGFO);

    log.info(String.format("Created read-group ID=%s PL=%s LB=%s SM=%s%n", rg.getId(), rg.getPlatform(), rg.getLibrary(), rg.getSample()));

    // create the new header and output file
    final SAMFileHeader inHeader = in.getFileHeader();
    final SAMFileHeader outHeader = inHeader.clone();
    outHeader.setReadGroups(Collections.singletonList(rg));
    if (SORT_ORDER != null) outHeader.setSortOrder(SORT_ORDER);

    final SAMFileWriter outWriter = new SAMFileWriterFactory().makeSAMOrBAMWriter(outHeader,
            outHeader.getSortOrder() == inHeader.getSortOrder(),
            OUTPUT);

    final ProgressLogger progress = new ProgressLogger(log);
    for (final SAMRecord read : in) {
        read.setAttribute(SAMTag.RG.name(), RGID);
        outWriter.addAlignment(read);
        progress.record(read);
    }

    // cleanup
    CloserUtil.close(in);
    outWriter.close();
    return 0;
}
 

@Test
public void testOutput() {
    final int length = 100;

    List<Byte> quals = new ArrayList<>(QualityUtils.MAX_SAM_QUAL_SCORE + 1);
    List<Long> counts = new ArrayList<>(QualityUtils.MAX_SAM_QUAL_SCORE + 1);

    for (int i = 0;  i<= QualityUtils.MAX_SAM_QUAL_SCORE; i++) {
        quals.add((byte) i);
        counts.add(1L);
    }

    final QuantizationInfo quantizationInfo = new QuantizationInfo(quals, counts);
    final RecalibrationArgumentCollection RAC = new RecalibrationArgumentCollection();

    quantizationInfo.noQuantization();
    final String readGroupID = "id";
    final StandardCovariateList covariateList = new StandardCovariateList(RAC, Collections.singletonList(readGroupID));

    final SAMReadGroupRecord rg = new SAMReadGroupRecord(readGroupID);
    rg.setPlatform("illumina");
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeaderWithReadGroup(rg);

    final GATKRead read = ArtificialReadUtils.createRandomRead(header, length, false);
    read.setReadGroup(rg.getReadGroupId());

    final byte [] readQuals = new byte[length];
    for (int i = 0; i < length; i++)
        readQuals[i] = 20;
    read.setBaseQualities(readQuals);

    final int expectedKeys = expectedNumberOfKeys(length, RAC.INDELS_CONTEXT_SIZE, RAC.MISMATCHES_CONTEXT_SIZE);
    int nKeys = 0;  // keep track of how many keys were produced
    final ReadCovariates rc = RecalUtils.computeCovariates(read, header, covariateList, true, new CovariateKeyCache());

    final RecalibrationTables recalibrationTables = new RecalibrationTables(covariateList);
    final NestedIntegerArray<RecalDatum> rgTable = recalibrationTables.getReadGroupTable();
    final NestedIntegerArray<RecalDatum> qualTable = recalibrationTables.getQualityScoreTable();

    for (int offset = 0; offset < length; offset++) {

        for (EventType errorMode : EventType.values()) {

            final int[] covariates = rc.getKeySet(offset, errorMode);
            final int randomMax = errorMode == EventType.BASE_SUBSTITUTION ? 10000 : 100000;

            rgTable.put(createRandomRecalDatum(randomMax, 10), covariates[0], errorMode.ordinal());
            qualTable.put(createRandomRecalDatum(randomMax, 10), covariates[0], covariates[1], errorMode.ordinal());
            nKeys += 2;
            for (NestedIntegerArray<RecalDatum> covTable : recalibrationTables.getAdditionalTables()){
                Covariate cov = recalibrationTables.getCovariateForTable(covTable);
                final int covValue = covariates[covariateList.indexByClass(cov.getClass())];
                if ( covValue >= 0 ) {
                    covTable.put(createRandomRecalDatum(randomMax, 10), covariates[0], covariates[1], covValue, errorMode.ordinal());
                    nKeys++;
                }
            }
        }
    }
    Assert.assertEquals(nKeys, expectedKeys);
}
 

@Test
public void testCovariateGeneration() {
    final RecalibrationArgumentCollection RAC = new RecalibrationArgumentCollection();

    final String[] readGroups = {"RG1", "RG2", "RGbla"};
    ReadGroupCovariate rgCov = new ReadGroupCovariate(RAC, Arrays.asList(readGroups));
    QualityScoreCovariate qsCov = new QualityScoreCovariate(RAC);
    ContextCovariate coCov = new ContextCovariate(RAC);
    CycleCovariate cyCov = new CycleCovariate(RAC);

    StandardCovariateList covariates = new StandardCovariateList(RAC, Arrays.asList(readGroups));

    final int NUM_READS = 100;
    final Random rnd = Utils.getRandomGenerator();
    final CovariateKeyCache keyCache= new CovariateKeyCache();

    for (int idx = 0; idx < NUM_READS; idx++) {
        for (final String readGroupID : readGroups) {
            final SAMReadGroupRecord readGroupRecord = new SAMReadGroupRecord(readGroupID);
            readGroupRecord.setPlatform("illumina");
            final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeaderWithReadGroup(readGroupRecord);

            final int length = 10 + rnd.nextInt(100); // random read length, at least 10 bp long
            final GATKRead read = ArtificialReadUtils.createRandomRead(header, length, false);
            read.setIsReverseStrand(rnd.nextBoolean());
            read.setReadGroup(readGroupID);

            final byte[] mQuals = read.getBaseQualities();
            final byte[] iQuals = ReadUtils.getBaseInsertionQualities(read);
            final byte[] dQuals = ReadUtils.getBaseDeletionQualities(read);
            ReadCovariates rc = RecalUtils.computeCovariates(read, header, covariates, true, keyCache);

            // check that the length is correct
            Assert.assertEquals(rc.getMismatchesKeySet().length, length);
            Assert.assertEquals(rc.getInsertionsKeySet().length, length);
            Assert.assertEquals(rc.getDeletionsKeySet().length, length);

            for (int i = 0; i < length; i++) {
                // check that read group is always the same
                Assert.assertEquals(rgCov.formatKey(rc.getMismatchesKeySet(i)[0]), readGroupID);
                Assert.assertEquals(rgCov.formatKey(rc.getInsertionsKeySet(i)[0]), readGroupID);
                Assert.assertEquals(rgCov.formatKey(rc.getDeletionsKeySet(i)[0]), readGroupID);

                // check quality score
                Assert.assertEquals(qsCov.formatKey(rc.getMismatchesKeySet(i)[1]), String.valueOf(mQuals[i]));
                Assert.assertEquals(qsCov.formatKey(rc.getInsertionsKeySet(i)[1]), String.valueOf(iQuals[i]));
                Assert.assertEquals(qsCov.formatKey(rc.getDeletionsKeySet(i)[1]),  String.valueOf(dQuals[i]));

                // check context
                Assert.assertEquals(coCov.formatKey(rc.getMismatchesKeySet(i)[2]), ContextCovariateUnitTest.expectedContext(read, i, RAC.MISMATCHES_CONTEXT_SIZE, RAC.LOW_QUAL_TAIL), "read: " +idx  + " readGroup:" + readGroupID + " context mismatch key at position:" + i);
                Assert.assertEquals(coCov.formatKey(rc.getInsertionsKeySet(i)[2]), ContextCovariateUnitTest.expectedContext(read, i, RAC.INDELS_CONTEXT_SIZE, RAC.LOW_QUAL_TAIL), "read: " +idx  + " readGroup:" + readGroupID + " context insertion key at position:" + i);
                Assert.assertEquals(coCov.formatKey(rc.getDeletionsKeySet(i)[2]),  ContextCovariateUnitTest.expectedContext(read, i, RAC.INDELS_CONTEXT_SIZE, RAC.LOW_QUAL_TAIL), "read: " +idx  + " readGroup:" + readGroupID + " context deletion key at position:" + i);

                // check cycle
                final int expectedCycleMismatch = CycleCovariateUnitTest.expectedCycle(read, i, false, RAC.MAXIMUM_CYCLE_VALUE);
                final int expectedCycleIndel = CycleCovariateUnitTest.expectedCycle(read, i, true, RAC.MAXIMUM_CYCLE_VALUE);
                Assert.assertEquals(cyCov.formatKey(rc.getMismatchesKeySet(i)[3]), String.valueOf(expectedCycleMismatch), "read: " + idx + " cycle mismatch key at position:" + i);
                Assert.assertEquals(cyCov.formatKey(rc.getInsertionsKeySet(i)[3]), String.valueOf(expectedCycleIndel),  "read: " +idx + " cycle insertion key at position:" + i);
                Assert.assertEquals(cyCov.formatKey(rc.getDeletionsKeySet(i)[3]), String.valueOf(expectedCycleIndel),  "read: " +idx + " cycle deletion key at position:" + i);
            }

        }

    }

}