Java Code Examples for htsjdk.samtools.SAMRecord

The following examples show how to use htsjdk.samtools.SAMRecord. 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: TagValueFilteringIteratorTest.java    License: MIT License 6 votes vote down vote up
@Test
public void testCellBarcodeFiltering () {
	String cellBarcodeTag = "XC";

	SamReader inputSam = SamReaderFactory.makeDefault().open(IN_FILE);
	String [] desiredBarcodes = {"ATCAGGGACAGA", "TTGCCTTACGCG", "TACAATTAAGGC"};

	TagValueFilteringIterator<String> iter = new TagValueFilteringIterator<String>(inputSam.iterator(), cellBarcodeTag, Arrays.asList(desiredBarcodes));
	Set<String> barcodesFound = new HashSet<String>();
	while (iter.hasNext()) {
		SAMRecord r = iter.next();
		String v = r.getStringAttribute(cellBarcodeTag);
		barcodesFound.add(v);
	}

	// should be the same size.
	Assert.assertEquals(barcodesFound.size(), desiredBarcodes.length);
	// should contain the same answers.
	for (String s: desiredBarcodes)
		Assert.assertTrue(barcodesFound.contains(s));

}
 
Example 2
Source Project: abra2   Source File: SAMRecordUtils.java    License: MIT License 6 votes vote down vote up
public static String getTrailingClips(SAMRecord read) {
	List<CigarElement> elems = read.getCigar().getCigarElements();
	
	List<CigarElement> trailing = new ArrayList<CigarElement>();
	boolean isNonClippedReached = false;
	
	for (CigarElement elem : elems) {
		if (isClip(elem)) {
			if (isNonClippedReached) {
				trailing.add(elem);
			}
		} else {
			isNonClippedReached = true;
		}
	}

	String ret = "";
	if (trailing.size() > 0) {
		Cigar cigar = new Cigar(trailing);
		ret =  TextCigarCodec.encode(cigar);
	}
	
	return ret;
}
 
Example 3
Source Project: Drop-seq   Source File: DetectBeadSubstitutionErrorsTest.java    License: MIT License 6 votes vote down vote up
/**
  * Very lame test -- just confirms that program doesn't crash and return 0 exit status.
  */
 @Test
 public void testBasic() throws IOException {
     final DetectBeadSubstitutionErrors clp = new DetectBeadSubstitutionErrors();
     clp.INPUT = Arrays.asList(TEST_FILE);
     clp.OUTPUT_REPORT = File.createTempFile("DetectBeadSubstitutionErrorsTest.", ".substitution_report.txt");
     clp.OUTPUT_REPORT.deleteOnExit();
     clp.OUTPUT = File.createTempFile("DetectBeadSubstitutionErrorsTest.", ".sam");
     clp.OUTPUT.deleteOnExit();
     clp.OUTPUT_SUMMARY = File.createTempFile("DetectBeadSubstitutionErrorsTest.", "substitution_summary.txt");
     clp.OUTPUT_SUMMARY.deleteOnExit();
     Assert.assertEquals(clp.doWork(), 0);
     final TabbedTextFileWithHeaderParser parser = new TabbedTextFileWithHeaderParser(clp.OUTPUT_REPORT);
     int numRows = 0;
     for (final TabbedTextFileWithHeaderParser.Row row : parser) {
         ++numRows;
         Assert.assertEquals(row.getField("intended_barcode"), BIG_BARCODE);
         final String neighbor_barcode = row.getField("neighbor_barcode");
         Assert.assertTrue(SMALL_BARCODES.contains(neighbor_barcode));
         Assert.assertTrue(Boolean.parseBoolean(row.getField("repaired")));
     }
     final SamReader reader = SamReaderFactory.makeDefault().open(clp.OUTPUT);
     for (final SAMRecord rec : reader)
Assert.assertEquals(rec.getStringAttribute("XC"), BIG_BARCODE);
     CloserUtil.close(reader);
 }
 
Example 4
Source Project: cramtools   Source File: SamRecordComparision.java    License: Apache License 2.0 6 votes vote down vote up
/**
 * This is supposed to check if the mates have valid pairing flags.
 * 
 * @param r1
 * @param r2
 * @return
 */
private boolean checkMateFlags(SAMRecord r1, SAMRecord r2) {
	if (!r1.getReadPairedFlag() || !r2.getReadPairedFlag())
		return false;

	if (r1.getReadUnmappedFlag() != r2.getMateUnmappedFlag())
		return false;
	if (r1.getReadNegativeStrandFlag() != r2.getMateNegativeStrandFlag())
		return false;
	if (r1.getProperPairFlag() != r2.getProperPairFlag())
		return false;
	if (r1.getFirstOfPairFlag() && r2.getFirstOfPairFlag())
		return false;
	if (r1.getSecondOfPairFlag() && r2.getSecondOfPairFlag())
		return false;

	if (r2.getReadUnmappedFlag() != r1.getMateUnmappedFlag())
		return false;
	if (r2.getReadNegativeStrandFlag() != r1.getMateNegativeStrandFlag())
		return false;

	return true;
}
 
Example 5
Source Project: Hadoop-BAM   Source File: TestBAMInputFormat.java    License: MIT License 6 votes vote down vote up
@Test
public void testMultipleSplits() throws Exception {
  input = BAMTestUtil.writeBamFile(1000, SAMFileHeader.SortOrder.queryname)
      .getAbsolutePath();
  completeSetup();
  jobContext.getConfiguration().setInt(FileInputFormat.SPLIT_MAXSIZE, 40000);
  BAMInputFormat inputFormat = new BAMInputFormat();
  List<InputSplit> splits = inputFormat.getSplits(jobContext);
  assertEquals(2, splits.size());
  List<SAMRecord> split0Records = getSAMRecordsFromSplit(inputFormat, splits.get(0));
  List<SAMRecord> split1Records = getSAMRecordsFromSplit(inputFormat, splits.get(1));
  assertEquals(1577, split0Records.size());
  assertEquals(425, split1Records.size());
  SAMRecord lastRecordOfSplit0 = split0Records.get(split0Records.size() - 1);
  SAMRecord firstRecordOfSplit1 = split1Records.get(0);
  assertEquals(lastRecordOfSplit0.getReadName(), firstRecordOfSplit1.getReadName());
  assertTrue(lastRecordOfSplit0.getFirstOfPairFlag());
  assertTrue(firstRecordOfSplit1.getSecondOfPairFlag());
}
 
Example 6
Source Project: picard   Source File: MergeBamAlignmentTest.java    License: MIT License 6 votes vote down vote up
private void addAlignmentsForBestFragmentMapqStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final String sequence, final int[] mapqs) {
    boolean reverse = false;
    int alignmentStart = 1;
    for (final int mapq : mapqs) {
        final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
        alignedRecord.setReadName(unmappedRecord.getReadName());
        alignedRecord.setReadBases(unmappedRecord.getReadBases());
        alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
        alignedRecord.setReferenceName(sequence);
        alignedRecord.setAlignmentStart(alignmentStart);
        alignmentStart += 10; // Any old position will do
        alignedRecord.setReadNegativeStrandFlag(reverse);
        reverse = !reverse;
        alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
        alignedRecord.setMappingQuality(mapq);
        alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
        alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
        alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
        alignedRecord.setMateUnmappedFlag(true);
        writer.addAlignment(alignedRecord);
    }
}
 
Example 7
Source Project: abra2   Source File: AltContigGenerator.java    License: MIT License 6 votes vote down vote up
private boolean hasHighQualitySoftClipping(SAMRecord read, int start, int length) {
	
	int numHighQualBases = 0;
	int requiredHighQualBases = (int) (softClipFraction * length);
	
	for (int bq = start; bq < start+length; bq++) {
		if (read.getBaseQualities()[bq] >= minBaseQual) {
			numHighQualBases += 1;
			
			if (numHighQualBases >= requiredHighQualBases) {
				return true;
			}
		}
	}
	
	return false;
}
 
Example 8
Source Project: picard   Source File: SortSam.java    License: MIT License 6 votes vote down vote up
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);
    final SamReader reader = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);
    ;
    reader.getFileHeader().setSortOrder(SORT_ORDER.getSortOrder());
    final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(reader.getFileHeader(), false, OUTPUT);
    writer.setProgressLogger(
            new ProgressLogger(log, (int) 1e7, "Wrote", "records from a sorting collection"));

    final ProgressLogger progress = new ProgressLogger(log, (int) 1e7, "Read");
    for (final SAMRecord rec : reader) {
        writer.addAlignment(rec);
        progress.record(rec);
    }

    log.info("Finished reading inputs, merging and writing to output now.");

    CloserUtil.close(reader);
    writer.close();
    return 0;
}
 
Example 9
Source Project: picard   Source File: QuerySortedReadPairIteratorUtilTest.java    License: MIT License 6 votes vote down vote up
@Test
public void testBasicPairedRead() {
    SAMRecordSetBuilder builder = new SAMRecordSetBuilder(false, SAMFileHeader.SortOrder.queryname);
    builder.setReadLength(READ_LENGTH);
    builder.addPair("mapped_paired", 1, 1, 31);
    PeekableIterator<SAMRecord> iterator = new PeekableIterator<SAMRecord>(builder.iterator());

    QuerySortedReadPairIteratorUtil.ReadPair pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNotNull(pair);
    Assert.assertNotNull(pair.read1);
    Assert.assertNotNull(pair.read2);
    Assert.assertEquals("mapped_paired", pair.read1.getReadName());
    Assert.assertEquals("mapped_paired", pair.read2.getReadName());

    pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNull(pair);
}
 
Example 10
Source Project: picard   Source File: CollectDuplicateMetricsTester.java    License: MIT License 6 votes vote down vote up
@Override
public File getOutput() {
    final File output;
    try {
        output = File.createTempFile("CollectDuplicateMetrics", ".sam");
    } catch (IOException e) {
        throw new PicardException("problems creating output file", e);
    }
    output.deleteOnExit();

    try(SAMFileWriter writer = new SAMFileWriterFactory().makeWriter(samRecordSetBuilder.getHeader(), true, output.toPath(), this.fastaFiles.get(samRecordSetBuilder.getHeader()))) {
        samRecordSetBuilder.getRecords().forEach(a -> {
            final SAMRecord b = a.deepCopy();
            b.setDuplicateReadFlag(duplicateFlags.get(samRecordToDuplicatesFlagsKey(a)));
            writer.addAlignment(b);
        });
    }

    return output;
}
 
Example 11
Source Project: hmftools   Source File: ReadContextCounter.java    License: GNU General Public License v3.0 6 votes vote down vote up
@NotNull
private RealignedContext realignmentContext(boolean realign, int readIndex, SAMRecord record) {
    if (!realign) {
        return new RealignedContext(RealignedType.NONE, 0);
    }

    int index = readContext.readBasesPositionIndex();
    int leftIndex = readContext.readBasesLeftCentreIndex();
    int rightIndex = readContext.readBasesRightCentreIndex();

    int leftOffset = index - leftIndex;
    int rightOffset = rightIndex - index;

    int indelLength = indelLength(record);
    return Realigned.realignedAroundIndex(readContext,
            readIndex,
            record.getReadBases(),
            Math.max(indelLength + Math.max(leftOffset, rightOffset), Realigned.MAX_REPEAT_SIZE));
}
 
Example 12
Source Project: Drop-seq   Source File: SamRecordSortingIteratorFactory.java    License: MIT License 6 votes vote down vote up
/**
   * @param progressLogger pass null if not interested in progress.
   * @return An iterator with all the records from underlyingIterator, in order defined by comparator.
   */
  public static CloseableIterator<SAMRecord> create(final SAMFileHeader header,
                                         final Iterator<SAMRecord> underlyingIterator,
                                         final Comparator<SAMRecord> comparator,
                                         final ProgressLogger progressLogger) {
      final SortingIteratorFactory.ProgressCallback<SAMRecord> progressCallback;
      if (progressLogger != null)
	progressCallback = new SortingIteratorFactory.ProgressCallback<SAMRecord>() {
              @Override
              public void logProgress(final SAMRecord record) {
                  progressLogger.record(record);
              }
          };
else
	progressCallback = null;
      return SortingIteratorFactory.create(SAMRecord.class,
              underlyingIterator, comparator, new BAMRecordCodec(header),
              SAMFileWriterImpl.getDefaultMaxRecordsInRam(),
              progressCallback);
  }
 
Example 13
Source Project: Drop-seq   Source File: CellBarcodeFilteringIteratorTest.java    License: MIT License 6 votes vote down vote up
@Test
public void filterOut2() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> cellBarcodesToKeep = null;
	String cellBCTag ="XC";

	CellBarcodeFilteringIterator f = new CellBarcodeFilteringIterator(underlyingIterator, cellBCTag, cellBarcodesToKeep);
	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("1", 0, 1, false);
	SAMRecord r1 = b.getRecords().iterator().next();

	r1.setAttribute(cellBCTag, "CELL1");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(cellBCTag, "CELL2");
	Assert.assertFalse(f.filterOut(r1));

}
 
Example 14
Source Project: abra2   Source File: NativeAssembler.java    License: MIT License 6 votes vote down vote up
private int maxSoftClipLength(SAMRecord read, Feature region) {
	int len = 0;
	if (read.getCigarLength() > 1) {
		
		CigarElement elem = read.getCigar().getCigarElement(0);
		if (elem.getOperator() == CigarOperator.S && read.getAlignmentStart() >= region.getStart()-readLength) {
			len = elem.getLength();
		}
		
		elem = read.getCigar().getCigarElement(read.getCigarLength()-1);
		if (elem.getOperator() == CigarOperator.S && elem.getLength() > len && read.getAlignmentEnd() <= region.getEnd()+readLength) {
			len = elem.getLength();
		}			
	}
	
	return len;
}
 
Example 15
@Override
public SAMRecord next() {
  if (mNextRecord == null) {
    throw new NoSuchElementException();
  }
  final SAMRecord ret = mNextRecord;
  populateNext();
  return ret;
}
 
Example 16
Source Project: rtg-tools   Source File: OrientationTest.java    License: BSD 2-Clause "Simplified" License 5 votes vote down vote up
public static SAMRecord makeTestRecord(SAMFileHeader sfh, int flags, String sequence, String mateSequence, int pos, int matePos, String readGroup) {
  final SAMRecord sr = new SAMRecord(sfh);
  sr.setFlags(flags);
  sr.setReferenceName(sequence);
  sr.setMateReferenceName(mateSequence);
  sr.setAlignmentStart(pos);
  sr.setMateAlignmentStart(matePos);
  sr.setAttribute("RG", readGroup);
  if (pos < matePos) {
    sr.setInferredInsertSize(matePos - pos + 5);
  } else {
    sr.setInferredInsertSize(matePos - pos  - 5); //-1 * (pos - matePos + 5));
  }
  return sr;
}
 
Example 17
private static boolean matchesString(@NotNull final SAMRecord record, int index, @NotNull final String expected) {
    for (int i = 0; i < expected.length(); i++) {
        if ((byte) expected.charAt(i) != record.getReadBases()[index + i]) {
            return false;
        }
    }
    return true;
}
 
Example 18
private static String getReadString(final GATKRead read){
    final byte[] readBases = read.getBases();
    if (readBases.length == 0) {
        return SAMRecord.NULL_SEQUENCE_STRING;
    }
    return new String(readBases);
}
 
Example 19
Source Project: Drop-seq   Source File: ReadDuplicateWrapper.java    License: MIT License 5 votes vote down vote up
@Override
protected void processRecord(final SAMRecord rec) {
	// String recName = rec.getReadName();
	int coordinate = getCoordinate(rec);
	//Interval i = new Interval(rec.getReferenceName(), coordinate, coordinate);
	//String value = IntervalTagComparator.toString(i);
	// encode the string manually as an interval, don't use the full encoding as it's not needed!
	String value = rec.getReferenceName() + IntervalTagComparator.ENCODE_DELIMITER + coordinate;
	rec.setAttribute(this.tag, value);
	queueRecordForOutput(rec);
}
 
Example 20
Source Project: hmftools   Source File: FragmentSizeCalcs.java    License: GNU General Public License v3.0 5 votes vote down vote up
private boolean isCandidateRecord(final SAMRecord record)
{
    int fragmentLength = abs(record.getInferredInsertSize());
    if(fragmentLength > MAX_FRAGMENT_LENGTH)
        return false;

    // ignore translocations and inversions
    if(!record.getMateReferenceName().equals(record.getReferenceName()) || record.getMateNegativeStrandFlag() == record.getReadNegativeStrandFlag())
        return false;

    // ignore split and soft-clipped reads above the read length
    if(record.getCigar().containsOperator(CigarOperator.N) || !record.getCigar().containsOperator(CigarOperator.M))
        return false;

    int readLength = max(mMaxReadLength, mConfig.ReadLength);

    if(readLength > 0 && fragmentLength > readLength && record.getCigar().containsOperator(CigarOperator.S))
    {
        // permit small soft-clips up to a point and then none for longer fragments
        if(record.getCigar().getCigarElements().stream().anyMatch(x -> x.getOperator() == CigarOperator.S && x.getLength() > 2))
            return false;
    }

    // both reads must fall in the current gene
    int otherStartPos = record.getMateAlignmentStart();
    if(!positionWithin(otherStartPos, mCurrentGenesRange[SE_START], mCurrentGenesRange[SE_END]))
        return false;

    // reads cannot cover any part of an exon
    int posStart = record.getStart();
    int posEnd = record.getEnd();

    for(final TranscriptData transData : mCurrentTransDataList)
    {
        if(transData.exons().stream().anyMatch(x -> positionsOverlap(posStart, posEnd, x.ExonStart, x.ExonEnd)))
            return false;
    }

    return true;
}
 
Example 21
Source Project: cramtools   Source File: Cram2Bam.java    License: Apache License 2.0 5 votes vote down vote up
private static void setNextMate(CramCompressionRecord r, CramCompressionRecord next) {
	r.mateAlignmentStart = next.alignmentStart;
	r.setMateUnmapped(next.isSegmentUnmapped());
	r.setMateNegativeStrand(next.isNegativeStrand());
	r.mateSequenceID = next.sequenceId;
	if (r.mateSequenceID == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		r.mateAlignmentStart = SAMRecord.NO_ALIGNMENT_START;
}
 
Example 22
Source Project: picard   Source File: TargetMetricsCollector.java    License: MIT License 5 votes vote down vote up
@Override
protected PerUnitMetricCollector<METRIC_TYPE, Integer, SAMRecord> makeChildCollector(final String sample, final String library, final String readGroup) {
    final PerUnitTargetMetricCollector collector =  new PerUnitTargetMetricCollector(probeSetName, coverageByTargetForRead.keySet(),
            sample, library, readGroup, probeTerritory, targetTerritory, genomeSize,
            intervalToGc, minimumMappingQuality, minimumBaseQuality, clipOverlappingReads, includeIndels);
    if (this.probeSetName != null) {
        collector.setBaitSetName(probeSetName);
    }

    return collector;
}
 
Example 23
Source Project: abra2   Source File: SomaticLocusCaller.java    License: MIT License 5 votes vote down vote up
private boolean hasIndel(SAMRecord read, LocusInfo locus) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= locus.posStop && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case I:
				if (isWithin(refPosInRead, locus.posStart, locus.posStop)) {
					return true;
				}
				// Fall through to next case
			case S:
				readPos += elem.getLength();
				break;
			case D:
				if (isWithin(refPosInRead, locus.posStart, locus.posStop)) {
					return true;
				}
				// Fall through to next case
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				readPos += elem.getLength();
				refPosInRead += elem.getLength();
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return false;
}
 
Example 24
Source Project: Drop-seq   Source File: TagReadWithGeneFunctionTest.java    License: MIT License 5 votes vote down vote up
@Test
public void testIntronicCorrectCodingWrong () {
	// read on the negative strand
	SAMRecord r = getFakeRecord(testBAMFile, 91, 100, true);

	// gene with 2 exons, 1 coding from 1-10, one UTR from 91-100.  Positive strand gene.
	GeneFromGTF gene = new GeneFromGTF(r.getContig(), 1, 100, false, "A", "coding", "A", "coding", 1);
	final GeneFromGTF.TranscriptFromGTF tx = gene.addTranscript("trans1", 1, 100, 1, 100, 2, "trans1", "trans1", "coding");
	tx.addExon(1, 10);
	tx.addExon(91, 100);
	OverlapDetector<Gene> geneOverlapDetector = new OverlapDetector<>(0, 0);
	geneOverlapDetector.addLhs(gene, gene);

	// gene with 2 exons, 1 coding from 50-60, 1 coding from 150-160. Negative strand gene.
	GeneFromGTF gene2 = new GeneFromGTF(r.getContig(), 50, 160, true, "B", "coding", "B", "coding", 1);
	final GeneFromGTF.TranscriptFromGTF tx2 = gene2.addTranscript("trans2", 50, 160, 50, 150, 2, "trans2", "trans2", "coding");
	tx2.addExon(50, 60);
	tx2.addExon(150, 160);
	geneOverlapDetector.addLhs(gene2, gene2);

	TagReadWithGeneFunction tagger = new TagReadWithGeneFunction();
	List <Gene> genes = new ArrayList <> (geneOverlapDetector.getAll());
	Collections.sort(genes, TagReadWithGeneFunction.GENE_NAME_COMPARATOR);

	r=tagger.setAnnotations(r, geneOverlapDetector, false);
	String gn = r.getStringAttribute("gn");
	String gs = r.getStringAttribute("gs");
	String gf = r.getStringAttribute("gf");

	// names always come out alphabetically sorted.
	Assert.assertEquals(r.getStringAttribute("gn"), gene.getName()+","+gene2.getName());
	Assert.assertEquals(r.getStringAttribute("gs"), "+,-");
	Assert.assertEquals(r.getStringAttribute("gf"), LocusFunction.CODING.name() + "," + LocusFunction.INTRONIC.name());
	Assert.assertEquals(r.getStringAttribute("XF"), LocusFunction.INTRONIC.name());
}
 
Example 25
Source Project: chipster   Source File: SamBamUtils.java    License: MIT License 5 votes vote down vote up
public void normaliseBam(File bamFile, File normalisedBamFile) {

		// Read in a BAM file and its header
		SAMFileReader.setDefaultValidationStringency(ValidationStringency.SILENT);
		SAMFileReader reader = new SAMFileReader(IOUtil.openFileForReading(bamFile));
		SAMFileWriter writer = null;
		try {
			SAMFileHeader normalisedHeader = reader.getFileHeader();

			// Alter the chromosome names in header's SAMSequenceDictionary
			SAMSequenceDictionary normalisedDictionary = new SAMSequenceDictionary();
			for (SAMSequenceRecord sequenceRecord : normalisedHeader.getSequenceDictionary().getSequences()) {

				// Normalise chromosome
				String sequenceName = chromosomeNormaliser.normaliseChromosome(sequenceRecord.getSequenceName());
				normalisedDictionary.addSequence(new SAMSequenceRecord(sequenceName, sequenceRecord.getSequenceLength()));
			}
			normalisedHeader.setSequenceDictionary(normalisedDictionary);

			// Write new BAM file with normalised chromosome names
			writer = new SAMFileWriterFactory().makeBAMWriter(normalisedHeader, true, normalisedBamFile);
			for (final SAMRecord rec : reader) {
				rec.setHeader(normalisedHeader);
				writer.addAlignment(rec);
			}
			
		} finally {
			closeIfPossible(reader);
			closeIfPossible(writer);
		}
	}
 
Example 26
Source Project: Drop-seq   Source File: ParentEditDistanceMatcher.java    License: MIT License 5 votes vote down vote up
public TagValues getValues(final SAMRecord rec) {
    final TagValues ret = new TagValues(zeroEditDistanceIndices.length, nonZeroEditDistanceIndices.length);
    for (int i = 0; i < zeroEditDistanceIndices.length; ++i) {
        ret.zeroEditDistanceValues[i] = rec.getAttribute(countTag.get(zeroEditDistanceIndices[i]));
    }
    for (int i = 0; i < nonZeroEditDistanceIndices.length; ++i) {
        ret.nonZeroEditDistanceValues[i] = rec.getStringAttribute(countTag.get(nonZeroEditDistanceIndices[i]));
    }
    return ret;
}
 
Example 27
Source Project: abra2   Source File: NativeAssembler.java    License: MIT License 5 votes vote down vote up
private boolean hasLowQualityBase(SAMRecord read) {
	//TODO: Don't hardcode phred33
	for (int i=0; i<read.getBaseQualityString().length(); i++) {
		if ((read.getBaseQualityString().charAt(i) - '!') < 20) {
			return true;
		}
	}
	
	return false;
}
 
Example 28
private void checkReadNames( final File outputFile, final File reference, final List<String> expectedReadNames ) throws IOException {
    List<String> actualReadNames = new ArrayList<>();
    try ( final SamReader reader = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.SILENT).referenceSequence(reference).open(outputFile) ) {
        for ( SAMRecord read : reader ) {
            actualReadNames.add(read.getReadName());
        }
    }
    Assert.assertEquals(actualReadNames, expectedReadNames, "Read names in output do not match expected read names");
}
 
Example 29
Source Project: halvade   Source File: QualityEncoding.java    License: GNU General Public License v3.0 5 votes vote down vote up
public static SAMRecord fixMisencodedQuals(final SAMRecord read) throws QualityException {
    final byte[] quals = read.getBaseQualities();
    for ( int i = 0; i < quals.length; i++ ) {
        quals[i] -= fixQualityIlluminaToPhred;
        if ( quals[i] < 0 )
            throw new QualityException(quals[i]);
    }
    read.setBaseQualities(quals);
    return read;
}
 
Example 30
Source Project: picard   Source File: AdapterMarker.java    License: MIT License 5 votes vote down vote up
private void fixAlreadySeenReads() {
    //remove all the reads for the selected adapters
    Arrays.stream(adapters.get()).forEach(adapter -> preAdapterPrunedRecords.remove(adapter));
    //anything left is marked with the incorrect adapter and needs its XT tag removed
    preAdapterPrunedRecords.values().forEach(readList -> readList.parallelStream().forEach(read -> {
        Stream<SAMRecord.SAMTagAndValue> filterAttributes = read.getAttributes().stream().filter(tag -> !tag.tag.equals(ReservedTagConstants.XT));
        read.clearAttributes();
        filterAttributes.forEach(tag -> read.setAttribute(tag.tag, tag.value));
    }));
    preAdapterPrunedRecords.clear();
}