Java Code Examples for htsjdk.samtools.reference.ReferenceSequence

The following examples show how to use htsjdk.samtools.reference.ReferenceSequence. These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source Project: Drop-seq   Source File: MaskReferenceSequence.java    License: MIT License 6 votes vote down vote up
private void processByPartialContig (final ReferenceSequenceFile ref, final FastaSequenceFileWriter writer, final File intervalListFile) {
	SAMSequenceDictionary sd = ref.getSequenceDictionary();
	// validate that the intervals and the reference have the same sequence dictionary.
	IntervalList iList = IntervalList.fromFile(intervalListFile);
	iList.getHeader().getSequenceDictionary().assertSameDictionary(sd);
	// map the intervals to a map to each contig.
	Map<String, List<Interval>> intervalsPerContig = getIntervalsForContig(iList);

	for (SAMSequenceRecord r: sd.getSequences()) {
		String contig = r.getSequenceName();
		log.info("Processing partial contig " + contig);
		// this list can be null.
		List<Interval> intervalsToMask = intervalsPerContig.get(contig);
		ReferenceSequence rs = ref.getSequence(contig);
		writeSequence(rs, intervalsToMask, writer);
	}
}
 
Example 2
Source Project: VarDictJava   Source File: ReferenceResource.java    License: MIT License 6 votes vote down vote up
/**
 * Method convert fasta data to array contains header and nucleotide bases.
 * @param fasta path to fasta file
 * @param chr chromosome name of region
 * @param start start position of region
 * @param end end position of region
 * @return array of nucleotide bases in the region of fasta
 */
public String[] retrieveSubSeq(String fasta, String chr, int start, int end) {
    try {
        IndexedFastaSequenceFile idx = fetchFasta(fasta);
        ReferenceSequence seq = idx.getSubsequenceAt(chr, start, end);
        byte[] bases = seq.getBases();
        return new String[]{">" + chr + ":" + start + "-" + end, bases != null ? new String(bases) : ""};
    } catch (SAMException e){
        if (UNABLE_FIND_CONTIG.matcher(e.getMessage()).find()){
            throw new WrongFastaOrBamException(chr, e);
        } else if (WRONG_START_OR_END.matcher(e.getMessage()).find()){
            throw new RegionBoundariesException(chr, start, end, e);
        } else {
            throw e;
        }
    }
}
 
Example 3
Source Project: picard   Source File: GcBiasUtils.java    License: MIT License 6 votes vote down vote up
public static int[] calculateRefWindowsByGc(final int windows, final File referenceSequence, final int windowSize) {
    final ReferenceSequenceFile refFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceSequence);
    ReferenceSequence ref;

    final int [] windowsByGc = new int [windows];

    while ((ref = refFile.nextSequence()) != null) {
        final byte[] refBases = ref.getBases();
        StringUtil.toUpperCase(refBases);
        final int refLength = refBases.length;
        final int lastWindowStart = refLength - windowSize;

        final CalculateGcState state = new GcBiasUtils().new CalculateGcState();

        for (int i = 1; i < lastWindowStart; ++i) {
            final int windowEnd = i + windowSize;
            final int gcBin = calculateGc(refBases, i, windowEnd, state);
            if (gcBin != -1) windowsByGc[gcBin]++;
        }
    }

    return windowsByGc;
}
 
Example 4
/**
 * Create one SAMSequenceRecord from a single fasta sequence
 */
private SAMSequenceRecord makeSequenceRecord(final ReferenceSequence refSeq) {
  final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(), refSeq.length());

  // Compute MD5 of upcased bases
  final byte[] bases = refSeq.getBases();
  for (int i = 0; i < bases.length; ++i) {
    bases[i] = StringUtil.toUpperCase(bases[i]);
  }

  ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
  if (GENOME_ASSEMBLY != null) {
    ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, GENOME_ASSEMBLY);
  }
  ret.setAttribute(SAMSequenceRecord.URI_TAG, URI);
  if (SPECIES != null) {
    ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, SPECIES);
  }
  return ret;
}
 
Example 5
Source Project: picard   Source File: FastWgsMetricsCollector.java    License: MIT License 6 votes vote down vote up
@Override
public void addInfo(final AbstractLocusInfo<EdgingRecordAndOffset> info, final ReferenceSequence ref, boolean referenceBaseN) {
    prepareCollector(info);
    for (final EdgingRecordAndOffset record : info.getRecordAndOffsets()) {
        final String readName = record.getReadName();
        Optional<Set<EdgingRecordAndOffset>> recordsAndOffsetsForName = Optional.ofNullable((Set<EdgingRecordAndOffset>)readsNames.get(readName));
        if (record.getType() == EdgingRecordAndOffset.Type.BEGIN) {
            processRecord(info.getPosition(), ref, record, recordsAndOffsetsForName.orElse(new HashSet<>()));
        } else {
            recordsAndOffsetsForName.ifPresent(
                    edgingRecordAndOffsets -> removeRecordFromMap(record,
                            edgingRecordAndOffsets));
        }
    }
    if (!referenceBaseN) {
        final int readNamesSize = pileupSize.get(info.getPosition());
        final int highQualityDepth = Math.min(readNamesSize, coverageCap);
        if (highQualityDepth < readNamesSize) {
            basesExcludedByCapping += readNamesSize - coverageCap;
        }
        highQualityDepthHistogramArray[highQualityDepth]++;
        unfilteredDepthHistogramArray[unfilteredDepthSize.get(info.getPosition())]++;
    }
}
 
Example 6
Source Project: picard   Source File: BaitDesigner.java    License: MIT License 6 votes vote down vote up
@Override
List<Bait> design(final BaitDesigner designer, final Interval target, final ReferenceSequence reference) {
    final List<Bait> baits = new LinkedList<Bait>();
    final int baitSize = designer.BAIT_SIZE;
    final int baitOffset = designer.BAIT_OFFSET;
    final int lastPossibleBaitStart = Math.min(target.getEnd(), reference.length() - baitSize);
    final int baitCount = 1 + (int) Math.floor((lastPossibleBaitStart - target.getStart()) / (double) baitOffset);

    int i = 0;
    for (int start = target.getStart(); start < lastPossibleBaitStart; start += baitOffset) {
        final Bait bait = new Bait(target.getContig(),
                start,
                CoordMath.getEnd(start, baitSize),
                target.isNegativeStrand(),
                designer.makeBaitName(target.getName(), ++i, baitCount));
        bait.addBases(reference, designer.DESIGN_ON_TARGET_STRAND);
        baits.add(bait);
    }
    return baits;
}
 
Example 7
Source Project: picard   Source File: SequenceDictionaryUtils.java    License: MIT License 6 votes vote down vote up
/**
 * Create one SAMSequenceRecord from a single fasta sequence
 */
private SAMSequenceRecord makeSequenceRecord(final ReferenceSequence refSeq) {
    final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(), refSeq.length());

    // Compute MD5 of upcased bases
    final byte[] bases = refSeq.getBases();
    for (int i = 0; i < bases.length; ++i) {
        bases[i] = StringUtil.toUpperCase(bases[i]);
    }

    ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
    if (genomeAssembly != null) {
        ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, genomeAssembly);
    }
    ret.setAttribute(SAMSequenceRecord.URI_TAG, uri);
    if (species != null) {
        ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, species);
    }
    return ret;
}
 
Example 8
Source Project: picard   Source File: AbstractWgsMetricsCollectorTest.java    License: MIT License 6 votes vote down vote up
@Test
public void testForCollectorWithoutData(){
    long[] templateQualHistogram =  new long[127];
    long[] templateHistogramArray = new long[11];
    CollectWgsMetrics collectWgsMetrics = new CollectWgsMetrics();
    AbstractWgsMetricsCollector collector = new AbstractWgsMetricsCollector(collectWgsMetrics,
            10, createIntervalList()) {
        @Override
        public void addInfo(AbstractLocusInfo info, ReferenceSequence ref, boolean referenceBaseN) {
        }
    };
    assertEquals(templateHistogramArray, collector.highQualityDepthHistogramArray);
    assertEquals(templateQualHistogram, collector.unfilteredBaseQHistogramArray);
    assertEquals(0, collector.basesExcludedByCapping);
    assertEquals(0, collector.basesExcludedByOverlap);
    assertEquals(0, collector.basesExcludedByBaseq);
}
 
Example 9
Source Project: cramtools   Source File: Utils.java    License: Apache License 2.0 5 votes vote down vote up
public static byte[] getBasesFromReferenceFile(String referenceFilePath, String seqName, int from, int length) {
	ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(new File(
			referenceFilePath));
	ReferenceSequence sequence = referenceSequenceFile.getSequence(seqName);
	byte[] bases = referenceSequenceFile.getSubsequenceAt(sequence.getName(), from, from + length).getBases();
	return bases;
}
 
Example 10
Source Project: Drop-seq   Source File: ValidateReference.java    License: MIT License 5 votes vote down vote up
private void validateReferenceBases(File referenceFile) {
    final ReferenceSequenceFile refSeqFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceFile, true);
    ReferenceSequence sequence;
    while ((sequence = refSeqFile.nextSequence()) != null) {
        for (final byte base: sequence.getBases()) {
            if (!IUPAC_TABLE[base]) {
                messages.baseErrors = String.format("WARNING: AT least one invalid base '%c' (decimal %d) in reference sequence named %s",
                        StringUtil.byteToChar(base), base, sequence.getName());
                break;
            }
        }
    }
}
 
Example 11
@Test(dataProvider = "provideDataForGetCodonChangeString")
void testGetCodonChangeString( final String refAllele,
                               final String altAllele,
                               final int alleleStart,
                               final int codingSequenceAlleleStart,
                               final int alignedCodingSequenceAlleleStart,
                               final String alignedCodingSeqRefAllele,
                               final String alignedCodingSeqAltAllele,
                               final String alignedAlternateAllele,
                               final int alignedRefAlleleStop,
                               final String contig,
                               final Strand strand,
                               final ReferenceSequence codingSequence,
                               final Locatable startCodon,
                               final String expected ) {

    final SequenceComparison seqComp = new SequenceComparison();
    seqComp.setReferenceAllele(refAllele);
    seqComp.setAlternateAllele(altAllele);
    seqComp.setAlleleStart(alleleStart);
    seqComp.setCodingSequenceAlleleStart(codingSequenceAlleleStart);
    seqComp.setAlignedCodingSequenceAlleleStart(alignedCodingSequenceAlleleStart);
    seqComp.setAlignedCodingSequenceReferenceAllele(alignedCodingSeqRefAllele);
    seqComp.setAlignedCodingSequenceAlternateAllele(alignedCodingSeqAltAllele);
    seqComp.setAlignedAlternateAllele(alignedAlternateAllele);
    seqComp.setAlignedReferenceAlleleStop(alignedRefAlleleStop);
    seqComp.setContig(contig);
    seqComp.setStrand(strand);
    seqComp.setTranscriptCodingSequence(codingSequence);

    Assert.assertEquals(FuncotatorUtils.getCodonChangeString(seqComp, startCodon), expected);
}
 
Example 12
Source Project: Drop-seq   Source File: GatherGeneGCLength.java    License: MIT License 5 votes vote down vote up
private List<GCResult> calculateGCContentGene (final Gene gene, final ReferenceSequence fastaRef, final SAMSequenceDictionary dict) {
	List<GCResult> result = new ArrayList<>();

	for (Transcript t : gene) {
		String seq=getTranscriptSequence(t, fastaRef, dict);
		GCResult gc = new GCResult(seq);
		gc.setTranscript(t);

		// check for GQuadruplexes.
		List<GQuadruplex> gq = GQuadruplex.find(t.name, seq);
		gc.incrementGQuadruplexCount(gq.size());
		result.add(gc);
	}
	return (result);
}
 
Example 13
Source Project: Drop-seq   Source File: GatherGeneGCLength.java    License: MIT License 5 votes vote down vote up
public void writeTranscriptSequence (final Gene gene, final ReferenceSequence fastaRef, final SAMSequenceDictionary dict, final FastaSequenceFileWriter outSequence ) {

		for (Transcript t : gene) {
			String sequence=getTranscriptSequence(t, fastaRef, dict);
			outSequence.	writeSequence(gene.getName(), sequence, t.name);
		}
	}
 
Example 14
Source Project: Drop-seq   Source File: MaskReferenceSequence.java    License: MIT License 5 votes vote down vote up
private void processByWholeContig (final ReferenceSequenceFile ref, final FastaSequenceFileWriter writer, final List<String> contigPatternToIgnore) {
	SAMSequenceDictionary sd = ref.getSequenceDictionary();
	Set<String> contigsToIgnore = selectContigsToIgnore(sd, contigPatternToIgnore);
	// write out each contig.

	for (SAMSequenceRecord r: sd.getSequences()) {
		String contig = r.getSequenceName();
		log.info("Processing complete contig " + contig);
		ReferenceSequence rs = ref.getSequence(contig);
		boolean setSequenceToN = contigsToIgnore.contains(contig);
		writeSequence(rs, setSequenceToN, writer);
	}

}
 
Example 15
Source Project: Drop-seq   Source File: MaskReferenceSequence.java    License: MIT License 5 votes vote down vote up
private void writeSequence (final ReferenceSequence rs, final List<Interval> intervalsToMask, final FastaSequenceFileWriter writer) {
	String sequence = rs.getBaseString();
	if (intervalsToMask!=null && intervalsToMask.size()>0) {
		char [] seqArray = sequence.toCharArray();
		for (Interval i: intervalsToMask)
			for (int pos=i.getStart()-1; pos<i.getEnd(); pos++)
				seqArray[pos]='N';
		sequence=new String (seqArray);
	}
	writer.writeSequence(rs.getName(), sequence);
}
 
Example 16
private static void addTrinucleotideContext(@NotNull final ImmutableEnrichedSomaticVariant.Builder builder,
        @NotNull final SomaticVariant variant, @NotNull IndexedFastaSequenceFile reference) {
    final int chromosomeLength = reference.getSequenceDictionary().getSequence(variant.chromosome()).getSequenceLength();
    if (variant.position() < chromosomeLength) {
        final ReferenceSequence sequence =
                reference.getSubsequenceAt(variant.chromosome(), Math.max(1, variant.position() - 1), variant.position() + 1);
        builder.trinucleotideContext(sequence.getBaseString());
    } else {
        LOGGER.warn("Requested ref sequence beyond contig length! variant = {}", variant);
    }
}
 
Example 17
@NotNull
private VariantHotspot correctRef(@NotNull final VariantHotspot hotspot) {
    final ReferenceSequence refSequence = hotspot.isSimpleInsert()
            ? sequenceFile.getSubsequenceAt(hotspot.chromosome(), hotspot.position(), hotspot.position())
            : sequenceFile.getSubsequenceAt(hotspot.chromosome(), hotspot.position(), hotspot.position() + hotspot.ref().length() - 1);

    return ImmutableVariantHotspotImpl.builder().from(hotspot).ref(refSequence.getBaseString()).build();
}
 
Example 18
Source Project: hmftools   Source File: RefSequence.java    License: GNU General Public License v3.0 5 votes vote down vote up
@VisibleForTesting
RefSequence(final ReferenceSequence sequence) {
    this.sequence = sequence;
    this.start = sequence.getContigIndex() + 1;
    this.end = start + sequence.getBases().length - 1;
    this.indexedBases = new IndexedBases(start, 0, sequence.getBases());
}
 
Example 19
private void testSequential(final CachingIndexedFastaSequenceFile caching, final Path unzipped, final int querySize) throws IOException {
    Assert.assertTrue(caching.isPreservingCase(), "testSequential only works for case preserving CachingIndexedFastaSequenceFile readers");

    //use an unzipped version of the fasta because it's slow to repeatedly query the zipped version without caching
    try( final ReferenceSequenceFile uncached = ReferenceSequenceFileFactory.getReferenceSequenceFile(unzipped)) {

        SAMSequenceRecord contig = uncached.getSequenceDictionary().getSequence(0);
        for (int i = 0; i < contig.getSequenceLength(); i += STEP_SIZE) {
            int start = i;
            int stop = start + querySize;
            if (stop <= contig.getSequenceLength()) {
                ReferenceSequence cachedVal = caching.getSubsequenceAt(contig.getSequenceName(), start, stop);
                ReferenceSequence uncachedVal = uncached.getSubsequenceAt(contig.getSequenceName(), start, stop);

                Assert.assertEquals(cachedVal.getName(), uncachedVal.getName());
                Assert.assertEquals(cachedVal.getContigIndex(), uncachedVal.getContigIndex());
                Assert.assertEquals(cachedVal.getBases(), uncachedVal.getBases());
            }
        }

        // asserts for efficiency.  We are going to make contig.length / STEP_SIZE queries
        // at each of range: start -> start + querySize against a cache with size of X.
        // we expect to hit the cache each time range falls within X.  We expect a hit
        // on the cache if range is within X.  Which should happen at least (X - query_size * 2) / STEP_SIZE
        // times.
        final int minExpectedHits = (int) Math.floor(
                (Math.min(caching.getCacheSize(), contig.getSequenceLength()) - querySize * 2.0) / STEP_SIZE);
        caching.printEfficiency(Level.DEBUG);
        Assert.assertTrue(caching.getCacheHits() >= minExpectedHits,
                          "Expected at least " + minExpectedHits + " cache hits but only got " + caching.getCacheHits());
    }

}
 
Example 20
@Override
public ReferenceBases getReferenceBases(final SimpleInterval interval) throws IOException {
    try ( ReferenceSequenceFile referenceSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(getReferencePath()) ) {
        ReferenceSequence sequence = referenceSequenceFile.getSubsequenceAt(interval.getContig(), interval.getStart(), interval.getEnd());
        return new ReferenceBases(sequence.getBases(), interval);
    }
}
 
Example 21
@Test(dataProvider = "ReferenceIntervalDataProvider")
public void testQueryAndPrefetch(final Path testReference, final SimpleInterval interval, final String expectedBases ) {
    try (ReferenceDataSource reference = new ReferenceFileSource(testReference))  {
        ReferenceSequence queryResult = reference.queryAndPrefetch(interval);

        Assert.assertEquals(new String(queryResult.getBases()), expectedBases,
                "Wrong bases returned from queryAndPrefetch() for interval " + interval);
    }
}
 
Example 22
@Override
protected InsertSizeMetricsCollectorArgs makeArg(SAMRecord samRecord, ReferenceSequence refSeq) {
    // inferred insert size is negative if the mate maps to lower position than the read, so use abs
    final int insertSize = Math.abs(samRecord.getInferredInsertSize());
    final SamPairUtil.PairOrientation orientation = SamPairUtil.getPairOrientation(samRecord);

    return new InsertSizeMetricsCollectorArgs(insertSize, orientation);
}
 
Example 23
Source Project: gatk   Source File: FuncotatorUtils.java    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
/**
 * Gets the next (+ strand) or previous (- strand) complete in-frame codon from the given {@link ReferenceSequence}
 * according to the current codon position and strand.
 * @param referenceSequence The {@link ReferenceSequence} containing the complete coding sequence for the transcript on which the current variant occurs.  Must not be {@code null}.
 * @param currentAlignedCodingSequenceAlleleStart The starting position (1-based, inclusive) of the current codon.  Must be > 0.
 * @param currentAlignedCodingSequenceAlleleStop The ending position (1-based, inclusive) of the current codon.  Must be > 0.
 * @param strand The {@link Strand} on which the current codon resides.  Must not be {@code null}.  Must not be {@link Strand#NONE}.
 * @return The next (+ strand) or previous (- strand) codon in frame with the current codon as specified by the given current codon positions.
 */
private static String getAdjacentReferenceCodon(final ReferenceSequence referenceSequence,
                                                final int currentAlignedCodingSequenceAlleleStart,
                                                final int currentAlignedCodingSequenceAlleleStop,
                                                final Strand strand) {

    Utils.nonNull( referenceSequence );
    ParamUtils.isPositive(currentAlignedCodingSequenceAlleleStart, "Genomic positions must be > 0.");
    ParamUtils.isPositive(currentAlignedCodingSequenceAlleleStop, "Genomic positions must be > 0.");
    assertValidStrand(strand);

    final String nextRefCodon;
    if ( strand == Strand.POSITIVE ) {

        // Add AminoAcid.CODON_LENGTH to get the "next" codon on the - strand:
        final int endex = currentAlignedCodingSequenceAlleleStop + AminoAcid.CODON_LENGTH;

        // Make sure we don't try to get bases after the end of our reference sequence:
        if ( endex >= referenceSequence.getBaseString().length() ) {
            nextRefCodon = "";
        }
        else {
            nextRefCodon = referenceSequence.getBaseString().substring(currentAlignedCodingSequenceAlleleStop, endex);
        }
    }
    else {
        // Make sure we don't try to get bases before the start of our reference sequence:
        if ( currentAlignedCodingSequenceAlleleStart == 1 ) {
            nextRefCodon = "";
        }
        else {
            // Subtract 1 because of 1-inclusive genomic positions
            // Subtract AminoAcid.CODON_LENGTH to get the "next" codon on the - strand:
            nextRefCodon = referenceSequence.getBaseString().substring(currentAlignedCodingSequenceAlleleStart - 1 - AminoAcid.CODON_LENGTH, currentAlignedCodingSequenceAlleleStart - 1);
        }
    }
    return nextRefCodon;
}
 
Example 24
Source Project: picard   Source File: MultiLevelCollector.java    License: MIT License 5 votes vote down vote up
/**
 * Construct a argument of ARGTYPE using the given SAMRecord and ReferenceSequence then pass
 * this value to all collectors that should include this record
 */
public void acceptRecord(final SAMRecord record, final ReferenceSequence refSeq) {
    final ARGTYPE arg = makeArg(record, refSeq);

    for(final Distributor collector : outputOrderedDistributors) {
        collector.acceptRecord(arg, record.getReadGroup());
    }
}
 
Example 25
Source Project: picard   Source File: CollectBaseDistributionByCycle.java    License: MIT License 5 votes vote down vote up
@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
    if ((PF_READS_ONLY) && (rec.getReadFailsVendorQualityCheckFlag())) {
        return;
    }
    if ((ALIGNED_READS_ONLY) && (rec.getReadUnmappedFlag())) {
        return;
    }
    if (rec.isSecondaryOrSupplementary()) {
        return;
    }
    hist.addRecord(rec);
}
 
Example 26
Source Project: picard   Source File: WgsMetricsProcessorImpl.java    License: MIT License 5 votes vote down vote up
/**
 * Method gets the data from iterator for each locus and processes it with the help of collector.
 */
@Override
public void processFile() {
    long counter = 0;

    while (iterator.hasNext()) {
        final AbstractLocusInfo<T> info = iterator.next();
        final ReferenceSequence ref = refWalker.get(info.getSequenceIndex());
        boolean referenceBaseN = collector.isReferenceBaseN(info.getPosition(), ref);
        collector.addInfo(info, ref, referenceBaseN);
        if (referenceBaseN) {
            continue;
        }

        progress.record(info.getSequenceName(), info.getPosition());
        if (collector.isTimeToStop(++counter)) {
            break;
        }
        collector.setCounter(counter);
    }
    // check that we added the same number of bases to the raw coverage histogram and the base quality histograms
    final long sumBaseQ = Arrays.stream(collector.unfilteredBaseQHistogramArray).sum();
    final long sumDepthHisto = LongStream.rangeClosed(0, collector.coverageCap).map(i -> (i * collector.unfilteredDepthHistogramArray[(int) i])).sum();
    if (sumBaseQ != sumDepthHisto) {
        log.error("Coverage and baseQ distributions contain different amount of bases!");
    }
}
 
Example 27
/**
 * Construct a argument of ARGTYPE using the given SAMRecord and ReferenceSequence then pass
 * this value to all collectors that should include this record
 */
public void acceptRecord(final SAMRecord record, final ReferenceSequence refSeq) {
    final ARGTYPE arg = makeArg(record, refSeq);

    for(final Distributor collector : outputOrderedDistributors) {
        collector.acceptRecord(arg, record.getReadGroup());
    }
}
 
Example 28
Source Project: picard   Source File: AlignmentSummaryMetricsCollector.java    License: MIT License 5 votes vote down vote up
public void addRecord(final SAMRecord record, final ReferenceSequence ref) {
    if (record.getNotPrimaryAlignmentFlag()) {
        // only want 1 count per read so skip non primary alignments
        return;
    }

    collectReadData(record);
    collectQualityData(record, ref);
}
 
Example 29
Source Project: picard   Source File: CollectSequencingArtifactMetrics.java    License: MIT License 5 votes vote down vote up
private String getRefContext(final ReferenceSequence ref, final int contextStartIndex, final int contextFullLength) {
    // cache the upper-cased string version of this reference so we don't need to create a string for every base in every read
    if (currentRefIndex != ref.getContigIndex()) {
        currentRefString = new String(ref.getBases()).toUpperCase();
        currentRefIndex = ref.getContigIndex();
    }
    return currentRefString.substring(contextStartIndex, contextStartIndex + contextFullLength);
}
 
Example 30
Source Project: picard   Source File: InsertSizeMetricsCollector.java    License: MIT License 5 votes vote down vote up
@Override
protected InsertSizeCollectorArgs makeArg(SAMRecord samRecord, ReferenceSequence refSeq) {
    final int insertSize = Math.abs(samRecord.getInferredInsertSize());
    final SamPairUtil.PairOrientation orientation = SamPairUtil.getPairOrientation(samRecord);

    return new InsertSizeCollectorArgs(insertSize, orientation);
}