use of htsjdk.samtools.SAMRecordIterator in project gridss by PapenfussLab.
the class SequentialCoverageAnnotator method createLookup.
private List<ReferenceCoverageLookup> createLookup(ProcessingContext context, List<SAMEvidenceSource> sources, int windowSize) {
List<ReferenceCoverageLookup> result = new ArrayList<>();
for (SAMEvidenceSource ses : sources) {
assert (ses.getSourceCategory() >= 0);
assert (ses.getSourceCategory() < context.getCategoryCount());
// one read-ahead thread per input file
SamReader reader = SamReaderFactory.makeDefault().open(ses.getFile());
SAMRecordIterator rawIterator = reader.iterator();
rawIterator.assertSorted(SortOrder.coordinate);
CloseableIterator<SAMRecord> sit = new AsyncBufferedIterator<SAMRecord>(rawIterator, ses.getFile().getName() + "-Coverage");
// close the async iterator first to prevent aysnc reading from a closed stream
toclose.add(sit);
toclose.add(rawIterator);
toclose.add(reader);
sit = new ProgressLoggingSAMRecordIterator(sit, new ProgressLogger(log, 10000000));
SequentialReferenceCoverageLookup sourceLookup = new SequentialReferenceCoverageLookup(sit, ses.getMetrics().getIdsvMetrics(), ses.getReadPairConcordanceCalculator(), windowSize, ses.getSourceCategory(), context.isFilterDuplicates());
context.registerBuffer(ses.getFile().getName(), sourceLookup);
result.add(sourceLookup);
}
return result;
}
use of htsjdk.samtools.SAMRecordIterator in project gridss by PapenfussLab.
the class ComputeSamTags method doWork.
@Override
protected int doWork() {
log.debug("Setting language-neutral locale");
java.util.Locale.setDefault(Locale.ROOT);
validateParameters();
SamReaderFactory readerFactory = SamReaderFactory.make();
SAMFileWriterFactory writerFactory = new SAMFileWriterFactory();
try {
try (SamReader reader = readerFactory.open(INPUT)) {
SAMFileHeader header = reader.getFileHeader();
if (!ASSUME_SORTED) {
if (header.getSortOrder() != SortOrder.queryname) {
log.error("INPUT is not sorted by queryname. " + "ComputeSamTags requires that reads with the same name be sorted together. " + "If the input file satisfies this constraint (the output from many aligners do)," + " this check can be disabled with the ASSUME_SORTED option.");
return -1;
}
}
try (SAMRecordIterator it = reader.iterator()) {
File tmpoutput = gridss.Defaults.OUTPUT_TO_TEMP_FILE ? FileSystemContext.getWorkingFileFor(OUTPUT, "gridss.tmp.ComputeSamTags.") : OUTPUT;
try (SAMFileWriter writer = writerFactory.makeSAMOrBAMWriter(header, true, tmpoutput)) {
compute(it, writer, getReference(), TAGS, SOFTEN_HARD_CLIPS, FIX_MATE_INFORMATION, RECALCULATE_SA_SUPPLEMENTARY, INPUT.getName() + "-");
}
if (tmpoutput != OUTPUT) {
FileHelper.move(tmpoutput, OUTPUT, true);
}
}
}
} catch (IOException e) {
log.error(e);
return -1;
}
return 0;
}
use of htsjdk.samtools.SAMRecordIterator in project gridss by PapenfussLab.
the class ByReadNameSinglePassSamProgram method makeItSo.
public static void makeItSo(final File input, final File referenceSequence, final boolean assumeSorted, final long stopAfter, final Collection<ByReadNameSinglePassSamProgram> programs) throws FileNotFoundException {
// Setup the standard inputs
IOUtil.assertFileIsReadable(input);
SamReader in = SamReaderFactory.makeDefault().referenceSequence(referenceSequence).open(input);
// Optionally load up the reference sequence and double check sequence dictionaries
final ReferenceLookup lookup;
if (referenceSequence == null) {
lookup = null;
} else {
IOUtil.assertFileIsReadable(referenceSequence);
lookup = new TwoBitBufferedReferenceSequenceFile(new IndexedFastaSequenceFile(referenceSequence));
if (!in.getFileHeader().getSequenceDictionary().isEmpty()) {
SequenceUtil.assertSequenceDictionariesEqual(in.getFileHeader().getSequenceDictionary(), lookup.getSequenceDictionary());
}
}
// Check on the sort order of the BAM file
final SortOrder sort = in.getFileHeader().getSortOrder();
if (sort != SortOrder.queryname) {
if (assumeSorted) {
log.warn("File reports sort order '" + sort + "', assuming it's queryname sorted anyway.");
} else {
throw new PicardException("File " + input.getAbsolutePath() + " should be queryname sorted but " + "the header says the sort order is " + sort + ". If you believe the file " + "to be queryname sorted you may pass ASSUME_SORTED=true");
}
}
for (final ByReadNameSinglePassSamProgram program : programs) {
program.setReference(lookup);
program.setup(in.getFileHeader(), input);
}
final ProgressLogger progress = new ProgressLogger(log);
final SAMRecordIterator rawit = in.iterator();
final CloseableIterator<SAMRecord> it = new AsyncBufferedIterator<SAMRecord>(rawit, "ByReadNameSinglePassSamProgram " + input.getName());
try {
List<SAMRecord> currentRecords = new ArrayList<>();
String currentReadName = null;
while (it.hasNext()) {
SAMRecord r = it.next();
String readname = r.getReadName();
// if read name we have to just treat it as a single read
if (readname == null || !readname.equals(currentReadName)) {
if (currentRecords.size() > 0) {
for (final ByReadNameSinglePassSamProgram program : programs) {
program.acceptFragment(currentRecords, lookup);
}
}
currentRecords.clear();
currentReadName = readname;
if (stopAfter > 0 && progress.getCount() >= stopAfter) {
break;
}
}
currentRecords.add(r);
progress.record(r);
}
if (currentRecords.size() > 0) {
for (final ByReadNameSinglePassSamProgram program : programs) {
program.acceptFragment(currentRecords, lookup);
}
}
} finally {
CloserUtil.close(it);
CloserUtil.close(rawit);
CloserUtil.close(in);
}
for (final ByReadNameSinglePassSamProgram program : programs) {
program.finish();
}
}
use of htsjdk.samtools.SAMRecordIterator in project gridss by PapenfussLab.
the class BamToBed method extract.
public static void extract(File input, File splitReads, File spanningReads, final double minLengthPortion, final int minMapq, final int minIndelSize, final int maxMappedBases, final int minIndelComponentSize) throws IOException {
SamReaderFactory factory = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
SamReader reader = factory.open(input);
BufferedWriter srwriter = null;
if (splitReads != null) {
srwriter = new BufferedWriter(new FileWriter(splitReads));
}
BufferedWriter spwriter = null;
if (spanningReads != null) {
spwriter = new BufferedWriter(new FileWriter(spanningReads));
}
SAMRecordIterator it = reader.iterator();
while (it.hasNext()) {
final SAMRecord r = it.next();
if (r.isSecondaryOrSupplementary())
continue;
if (splitReads != null) {
for (SimpleBEDFeature f : ChimericAlignment.getChimericAlignments(r).stream().map(ca -> getSplitReadDeletion(r, ca, minLengthPortion, minMapq, minIndelSize)).filter(f -> f != null).collect(Collectors.toList())) {
writeBed(srwriter, f);
}
}
if (spanningReads != null) {
for (SimpleBEDFeature f : getSpanningDeletion(r, minMapq, minIndelSize, maxMappedBases, minIndelComponentSize)) {
writeBed(spwriter, f);
}
}
}
CloserUtil.close(srwriter);
CloserUtil.close(spwriter);
}
use of htsjdk.samtools.SAMRecordIterator in project gridss by PapenfussLab.
the class LongReadSupportFinder method evaluateDeletion.
/**
* evaluate the long read support for the given deletion event
* @param chr
* @param start1
* @param end1
* @param start2
* @param end2
* @param length length of deletion
* @param softClipMargin error margin allowed for soft clip location
* @param spanningWindowMargin
* @param minDeletionSize
* @return
*/
public LongReadSupportLevel evaluateDeletion(String chr, int start1, int end1, int start2, int end2, int softClipMargin, int minSoftClipLength, int spanningWindowMargin, int minDeletionSize) {
chr = translateReference(chr);
if (chr == null)
return null;
SAMRecordIterator it = null;
try {
it = reader.query(chr, start1, end1, false);
LongReadSupportLevel support = new LongReadSupportLevel();
while (it.hasNext()) {
SAMRecord r = it.next();
if (SAMRecordUtil.getEndSoftClipLength(r) >= minSoftClipLength && IntervalUtil.overlapsClosed(r.getAlignmentEnd(), r.getAlignmentEnd(), start1 - softClipMargin, end1 + softClipMargin)) {
support.startClipLocations.add(r.getAlignmentEnd());
}
int deletions = countDeletions(r, minDeletionSize, start1, end2);
if (deletions > 0) {
support.spanningAlignments.add(deletions);
}
}
it.close();
it = null;
it = reader.query(chr, start2, end2, false);
while (it.hasNext()) {
SAMRecord r = it.next();
if (SAMRecordUtil.getStartSoftClipLength(r) >= minSoftClipLength && IntervalUtil.overlapsClosed(r.getAlignmentStart(), r.getAlignmentStart(), start2 - softClipMargin, end2 + softClipMargin)) {
support.endClipLocations.add(r.getAlignmentStart());
}
}
it.close();
it = null;
return support;
} catch (Exception e) {
e.printStackTrace();
return null;
} finally {
if (it != null) {
it.close();
}
}
}
Aggregations