Java SAMRecordIterator类代码示例

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
@Nullable
private Read getReadFromBamFile(ReadQuery query, Chromosome chromosome, BamFile bamFile) throws IOException {
    try (SamReader reader = bamHelper.makeSamReader(bamFile, Collections.singletonList(chromosome),
                                                    chromosome.getReferenceId())) {
        String chromosomeName = chromosome.getName();
        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }

        SAMRecordIterator iterator = reader.query(chromosomeName, query.getStartIndex(), query.getEndIndex(),
                                                  true);
        while (iterator.hasNext()) {
            final SAMRecord samRecord = iterator.next();
            if (samRecord.getReadName().equals(query.getName())) {
                BamQueryOption option = new BamQueryOption();
                option.setRefID(chromosome.getReferenceId());
                option.setChromosomeName(chromosome.getName());
                SAMRecordHandler recordHandler = new SAMRecordHandler(query.getStartIndex(), query.getEndIndex(),
                                                                      referenceManager, null, option);
                List<BasePosition> diffBase = recordHandler.computeDifferentBase(samRecord);
                return BamUtil.createExtendedRead(samRecord, diffBase);
            }
        }
    }
    return null;
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private void parseBam(final SamReader reader) {
    Assert.notNull(reader, getMessage(RESOURCE_NOT_FOUND));
    final SAMFileHeader samFileHeader = reader.getFileHeader();

    //check we can read this Bam-file
    //get list of chromosome
    final List<SAMSequenceRecord> list = samFileHeader.getSequenceDictionary().getSequences();
    Assert.notEmpty(list, getMessage(MessagesConstants.WRONG_HEADER_BAM_FILE_EMPTY_FILE));

    //get first chromosome and make a request to the file with this chromosome
    final SAMSequenceRecord samSequenceRecord = list.get(0);
    SAMRecordIterator iterator = reader.query(samSequenceRecord.getSequenceName(),
            Constants.BAM_START_INDEX_TEST, Math.min(Constants.MAX_BAM_END_INDEX_TEST,
                    samSequenceRecord.getSequenceLength()), false);
    Assert.notNull(iterator);
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private static File queryNameSortedBAM(final SamReader reader, final QueryInterval[] intervals, final String name) throws IOException {

        final SAMFileHeader header = reader.getFileHeader().clone();
        header.setSortOrder(SAMFileHeader.SortOrder.queryname);

        final File file = File.createTempFile(name, ".bam");
        final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(header, false, file);

        final SAMRecordIterator iterator = reader.queryOverlapping(intervals);
        while (iterator.hasNext()) {
            writer.addAlignment(iterator.next());
        }

        iterator.close();
        writer.close();

        return file;
    }
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
/**
 * Given a SamReader, find all reads within a region (ref:start-end) that span a splice junction and tally the number of times
 * each junction is spanned. Optionally tally the average values of "tallyTagName" tags (eg. NM).
 * 
 * @param reader
 * @param ref
 * @param start
 * @param end
 * @param orient
 * @param minOverlap
 * @param tallyTagName
 * @param separateReadCounts
 * @return
 */
public static SortedMap<GenomeSpan, MappedReadCounter> countJunctions(SamReader reader, String ref, int start, int end, Orientation orient, int minOverlap, String tallyTagName) {
    SAMRecordIterator it = reader.query(ref, start, end, true);
    SortedMap<GenomeSpan, MappedReadCounter> counters = new TreeMap<GenomeSpan, MappedReadCounter>();
    
    while (it.hasNext()) {
        SAMRecord read = it.next();
        if (read.isSecondaryOrSupplementary() || read.getDuplicateReadFlag() || read.getNotPrimaryAlignmentFlag() || read.getReadUnmappedFlag() || read.getSupplementaryAlignmentFlag()) {
            // skip all secondary / duplicate / unmapped reads
            continue;
        }
       
        if (isJunctionSpanning(read)) {
            for (GenomeSpan junction: getJunctionsForRead(read, orient, minOverlap)) {
                if (!counters.containsKey(junction)) {
                    counters.put(junction, new MappedReadCounter(tallyTagName));
                }
                
                counters.get(junction).addRead(read);
            }
        }
    }
    it.close();
    return counters;
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private void convertFileAndVerifyRecordCount(final int expectedCount,
                   final String fastqFilename1,
                   final String fastqFilename2,
                   final FastqQualityFormat version,
                   final boolean permissiveFormat,
                   final boolean useSequentialFastqs) throws IOException {
    final File samFile = convertFile(fastqFilename1, fastqFilename2, version, permissiveFormat, useSequentialFastqs);
    final SamReader samReader = SamReaderFactory.makeDefault().open(samFile);
    final SAMRecordIterator iterator = samReader.iterator();
    int actualCount = 0;
    while (iterator.hasNext()) {
        iterator.next();
        actualCount++;
    }
    samReader.close();
    Assert.assertEquals(actualCount, expectedCount);
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private static String compareReads(final File actualSam, final File expectedSam, final ValidationStringency validation, final File reference) throws IOException {
    try(final SamReader reader1 = getReader(actualSam, validation, reference);
        final SamReader reader2 = getReader(expectedSam, validation, reference)) {
        final SAMRecordIterator it1 = reader1.iterator();
        final SAMRecordIterator it2 = reader2.iterator();
        while (it1.hasNext() && it2.hasNext()) {
            final SAMRecord read1 = it1.next();
            final SAMRecord read2 = it2.next();
            final String eqMessage = readsEqualAllowAddingAttributes(read1, read2);
            if (eqMessage != null){
                return eqMessage;
            }
        }
        if (it1.hasNext() || it2.hasNext()) {
            //at least one has no more records (because the while loop is done) and at least one does have records. So we're not equal.
            return "Not the same number of reads";
        }
        return null;
    }
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private void scanPureJava() throws Exception
{
long start=System.currentTimeMillis();
   final MessageDigest md = MessageDigest.getInstance("MD5");
for(int N=0;N< N_LOOP;++N)
	{
	htsjdk.samtools.SamReader r=htsjdk.samtools.SamReaderFactory.makeDefault().
			validationStringency(ValidationStringency.DEFAULT_STRINGENCY).
			open(bamFile);
	SAMRecordIterator iter=r.iterator();
	while(iter.hasNext())
		{
		SAMRecord rec=iter.next();
		md.digest((rec.getReadName()+rec.getReadString()+rec.getBaseQualityString()+rec.getAlignmentStart()).getBytes());
		}
		
	r.close();
	}
System.out.println("PureJava\t\t\tmd5="+digest(md)+" time="+(System.currentTimeMillis()-start));
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
@Override
public SAMRecordIterator query(QueryInterval[] intervals, boolean contained) {
  final GA4GHQueryInterval[] myIntervals = new GA4GHQueryInterval[intervals.length];
  final GA4GHQueryInterval.ReadPositionConstraint constraint = contained ?
    GA4GHQueryInterval.ReadPositionConstraint.CONTAINED : 
    GA4GHQueryInterval.ReadPositionConstraint.OVERLAPPING;
  final SAMFileHeader header = getFileHeader();
  for (int i = 0; i < intervals.length; i++) {
    final QueryInterval interval = intervals[i];
    final String sequence = header.getSequence(
        interval.referenceIndex).getSequenceName();
    myIntervals[i] = new GA4GHQueryInterval(sequence, interval.start, 
        interval.end, constraint);
  }
  return query(myIntervals);
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private SortedBamIterator(File chunkDirectory, List<File> chunkFiles, final Comparator<SAMRecord> comparator) {
	this.chunkDirectory = chunkDirectory;
	this.comparator = comparator;
	currentSortedRecords = new TreeSet<SamRecordAndChunkIndex>(new Comparator<SamRecordAndChunkIndex>() {
		@Override
		public int compare(SamRecordAndChunkIndex o1, SamRecordAndChunkIndex o2) {
			return comparator.compare(o1.getRecord(), o2.getRecord());
		}
	});

	readers = new SamReader[chunkFiles.size()];
	iters = new SAMRecordIterator[chunkFiles.size()];
	for (int i = 0; i < chunkFiles.size(); i++) {
		SamReader currentReader = SamReaderFactory.makeDefault().open(chunkFiles.get(i));
		readers[i] = currentReader;
		SAMRecordIterator currentIter = currentReader.iterator();
		iters[i] = currentIter;
		if (currentIter.hasNext()) {
			currentSortedRecords.add(new SamRecordAndChunkIndex(currentIter.next(), i));
		}
	}

	nextRecord = getNextRecordToReturn();
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private SAMRecord getNextRecordToReturn() {
	SAMRecord nextRecord = null;
	if (currentSortedRecords.size() > 0) {
		SamRecordAndChunkIndex currentRecordAndReaderIndex = currentSortedRecords.pollFirst();
		nextRecord = currentRecordAndReaderIndex.getRecord();
		int chunkIndex = currentRecordAndReaderIndex.getChunkIndex();
		SAMRecordIterator currentIter = iters[chunkIndex];
		if (currentIter.hasNext()) {
			SAMRecord nextChunkRecord = currentIter.next();
			if (!currentSortedRecords.add(new SamRecordAndChunkIndex(nextChunkRecord, chunkIndex))) {
				// add this point we know there is a problem just collecting more info before reporting the problem.
				SamRecordAndChunkIndex alreadyContainedRecord = null;
				for (SamRecordAndChunkIndex samRecordAndChunk : currentSortedRecords) {
					if (comparator.compare(samRecordAndChunk.getRecord(), nextChunkRecord) == 0) {
						alreadyContainedRecord = samRecordAndChunk;
					}
				}
				throw new IllegalStateException("Provided comparator does not uniquely identify the records.  Comparator failed on sam record[" + nextChunkRecord.getReadName() + " isReadOne:"
						+ nextChunkRecord.getFirstOfPairFlag() + "] in chunk[" + chunkIndex + "] which matches [" + alreadyContainedRecord.getRecord().getReadName() + " isReadOne:"
						+ alreadyContainedRecord.getRecord().getFirstOfPairFlag() + "] in chunk[" + alreadyContainedRecord.getChunkIndex() + "].");
			}
		}
	}
	return nextRecord;
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private void scan(final SamReader r) 
{
final Row row=new Row();
SAMRecordIterator iter=null;
try{
	final SAMSequenceDictionaryProgress progress=new SAMSequenceDictionaryProgress(r.getFileHeader());
	iter=r.iterator();	
	while(iter.hasNext())
		{
		row.rec =progress.watch(iter.next());
		printAln(row);
		if(this.out.checkError()) break;
		}
	progress.finish();
	}
catch(final Exception err)
	{
	LOG.error("scan error:",err);
	throw new RuntimeException(String.valueOf(err.getMessage()),err);
	}
finally
	{
	CloserUtil.close(iter);
	}
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private void scan(final SamReader in) 
{
final SAMSequenceDictionary dict=in.getFileHeader().getSequenceDictionary();
if(dict==null) throw new RuntimeException("Sequence dictionary missing...");
final SAMRecordIterator iter=in.iterator();
final SAMSequenceDictionaryProgress progress=new SAMSequenceDictionaryProgress(dict);		
while(iter.hasNext() && !this.out.checkError())
	{
	final  SAMRecord rec=progress.watch(iter.next());
	if(rec.getReadUnmappedFlag()) continue;
	for(final PslAlign a:makePslAlign(rec,dict))
		{
		out.println(toString(a,rec));
		}
	}
progress.finish();
iter.close();
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private void scan(final SamReader in) throws IOException
 	{
 	SAMRecordIterator iter=in.iterator();
 	while(iter.hasNext())
{
 		final SAMRecord rec = iter.next();
 		String sampleName= this.partition.getPartion(rec);
 		if(sampleName==null || sampleName.isEmpty()) sampleName="__UNDEFINED_"+this.partition.name().toUpperCase()+"_";
 		
 		Counter<Integer> counter = this.lengths.get(sampleName);
 		if(counter==null) {
 			counter = new Counter<>();
 			 this.lengths.put(sampleName,counter);
 		}
 		final int len = rec.getReadLength();
 		counter.incr(len);
 		this.max_length=Math.max(max_length, len);
 		
}
 	iter.close();
 	}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private Map<Sequence, List<SAMRecord>> getReadsFromFile(final Chromosome chromosome, final Track<Sequence> track,
                                                        final BamFile bamFile, final List<Sequence> blocks)
        throws IOException {
    final Map<Sequence, List<SAMRecord>> records = new HashMap<>();
    try (SamReader reader = makeSamReader(bamFile, Collections.singletonList(chromosome),
            chromosome.getReferenceId())) {
        LOG.debug(getMessage(MessagesConstants.DEBUG_FILE_OPENING, bamFile.getPath()));
        String chromosomeName = chromosome.getName();
        final int startIndex = track.getStartIndex();
        final int endIndex = track.getEndIndex();

        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }
        final SAMRecordIterator iterator = reader.queryOverlapping(chromosomeName, startIndex, endIndex);

        LOG.debug(getMessage(MessagesConstants.DEBUG_GET_ITERATOR_QUERY, iterator.toString()));
        while (iterator.hasNext()) {
            final SAMRecord samRecord = iterator.next();
            //if read unmapped
            if (!samRecord.getSAMFlags().contains(SAMFlag.READ_UNMAPPED) && !samRecord.getCigar().isEmpty() &&
                    samRecord.getEnd() > samRecord.getStart()) {
                addRecordToBlock(blocks, records, samRecord);
            }
        }
    }
    return records;
}
 

import htsjdk.samtools.SAMRecordIterator; //导入依赖的package包/类
private static void writeToSlice(final String path, final SamReader reader, final QueryInterval[] intervals) {
    final File outputBAM = new File(path);
    final SAMFileWriter writer = new SAMFileWriterFactory().makeBAMWriter(reader.getFileHeader(), true, outputBAM);
    final SAMRecordIterator iterator = reader.queryOverlapping(intervals);
    while (iterator.hasNext()) {
        writer.addAlignment(iterator.next());
    }
    iterator.close();
    writer.close();
}