use of org.broadinstitute.hellbender.exceptions.UserException in project gatk by broadinstitute.
the class FeatureDataSource method getFeatureReader.
@SuppressWarnings("unchecked")
private static <T extends Feature> FeatureReader<T> getFeatureReader(final FeatureInput<T> featureInput, final Class<? extends Feature> targetFeatureType, final Function<SeekableByteChannel, SeekableByteChannel> cloudWrapper, final Function<SeekableByteChannel, SeekableByteChannel> cloudIndexWrapper, final Path reference) {
if (isGenomicsDBPath(featureInput.getFeaturePath())) {
try {
if (reference == null) {
throw new UserException.MissingReference("You must provide a reference if you want to load from GenomicsDB");
}
try {
final File referenceAsFile = reference.toFile();
return (FeatureReader<T>) getGenomicsDBFeatureReader(featureInput.getFeaturePath(), referenceAsFile);
} catch (final UnsupportedOperationException e) {
throw new UserException.BadInput("GenomicsDB requires that the reference be a local file.", e);
}
} catch (final ClassCastException e) {
throw new UserException("GenomicsDB inputs can only be used to provide VariantContexts.", e);
}
} else {
final Path featurePath = IOUtils.getPath(featureInput.getFeaturePath());
IOUtils.assertFileIsReadable(featurePath);
final FeatureCodec<T, ?> codec = (FeatureCodec<T, ?>) FeatureManager.getCodecForFile(featurePath, targetFeatureType);
return getTribbleFeatureReader(featureInput, codec, cloudWrapper, cloudIndexWrapper);
}
}
use of org.broadinstitute.hellbender.exceptions.UserException in project gatk by broadinstitute.
the class FeatureWalker method initializeDrivingFeatures.
@SuppressWarnings("unchecked")
private void initializeDrivingFeatures() {
final File drivingFile = getDrivingFeatureFile();
final FeatureCodec<? extends Feature, ?> codec = FeatureManager.getCodecForFile(drivingFile);
if (isAcceptableFeatureType(codec.getFeatureType())) {
drivingFeatures = new FeatureDataSource<>(new FeatureInput<>(drivingFile.getAbsolutePath()), FeatureDataSource.DEFAULT_QUERY_LOOKAHEAD_BASES, null, cloudPrefetchBuffer, cloudIndexPrefetchBuffer, referenceArguments.getReferencePath());
final FeatureInput<F> drivingFeaturesInput = new FeatureInput<>(drivingFile.getAbsolutePath(), "drivingFeatureFile");
features.addToFeatureSources(0, drivingFeaturesInput, codec.getFeatureType(), cloudPrefetchBuffer, cloudIndexPrefetchBuffer, referenceArguments.getReferencePath());
} else {
throw new UserException("File " + drivingFile + " contains features of the wrong type.");
}
}
use of org.broadinstitute.hellbender.exceptions.UserException in project gatk by broadinstitute.
the class SinglePassSamProgram method makeItSo.
public static void makeItSo(final File input, final File referenceSequence, final boolean assumeSorted, final long stopAfter, final Collection<SinglePassSamProgram> programs) {
// Setup the standard inputs
IOUtil.assertFileIsReadable(input);
final SamReader in = SamReaderFactory.makeDefault().referenceSequence(referenceSequence).open(input);
// Optionally load up the reference sequence and double check sequence dictionaries
final ReferenceSequenceFileWalker walker;
if (referenceSequence == null) {
walker = null;
} else {
IOUtil.assertFileIsReadable(referenceSequence);
walker = new ReferenceSequenceFileWalker(referenceSequence);
if (!in.getFileHeader().getSequenceDictionary().isEmpty()) {
SequenceUtil.assertSequenceDictionariesEqual(in.getFileHeader().getSequenceDictionary(), walker.getSequenceDictionary());
}
}
// Check on the sort order of the BAM file
{
final SortOrder sort = in.getFileHeader().getSortOrder();
if (sort != SortOrder.coordinate) {
if (assumeSorted) {
logger.warn("File reports sort order '" + sort + "', assuming it's coordinate sorted anyway.");
} else {
throw new UserException("File " + input.getAbsolutePath() + " should be coordinate sorted but " + "the header says the sort order is " + sort + ". If you believe the file " + "to be coordinate sorted you may pass ASSUME_SORTED=true");
}
}
}
// Call the abstract setup method!
boolean anyUseNoRefReads = false;
for (final SinglePassSamProgram program : programs) {
program.setup(in.getFileHeader(), input);
anyUseNoRefReads = anyUseNoRefReads || program.usesNoRefReads();
}
final ProgressLogger progress = new ProgressLogger(logger);
for (final SAMRecord rec : in) {
final ReferenceSequence ref;
if (walker == null || rec.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
ref = null;
} else {
ref = walker.get(rec.getReferenceIndex());
}
for (final SinglePassSamProgram program : programs) {
program.acceptRead(rec, ref);
}
progress.record(rec);
// See if we need to terminate early?
if (stopAfter > 0 && progress.getCount() >= stopAfter) {
break;
}
// And see if we're into the unmapped reads at the end
if (!anyUseNoRefReads && rec.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
break;
}
}
CloserUtil.close(in);
for (final SinglePassSamProgram program : programs) {
program.finish();
}
}
use of org.broadinstitute.hellbender.exceptions.UserException in project gatk by broadinstitute.
the class CollectSequencingArtifactMetrics method acceptRead.
@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
// see if the whole read should be skipped
if (recordFilter.filterOut(rec))
return;
// check read group + library
final String library = (rec.getReadGroup() == null) ? UNKNOWN_LIBRARY : getOrElse(rec.getReadGroup().getLibrary(), UNKNOWN_LIBRARY);
if (!libraries.contains(library)) {
// should never happen if SAM is valid
throw new UserException("Record contains library that is missing from header: " + library);
}
// iterate over aligned positions
for (final AlignmentBlock block : rec.getAlignmentBlocks()) {
for (int offset = 0; offset < block.getLength(); offset++) {
// remember, these are 1-based!
final int readPos = block.getReadStart() + offset;
final int refPos = block.getReferenceStart() + offset;
/**
* Skip regions outside of intervals.
*
* NB: IntervalListReferenceSequenceMask.get() has side-effects which assume
* that successive ReferenceSequence's passed to this method will be in-order
* (e.g. it will break if you call acceptRead() with chr1, then chr2, then chr1
* again). So this only works if the underlying iteration is coordinate-sorted.
*/
if (intervalMask != null && !intervalMask.get(ref.getContigIndex(), refPos))
continue;
// skip dbSNP sites
if (dbSnpMask != null && dbSnpMask.isDbSnpSite(ref.getName(), refPos))
continue;
// skip the ends of the reference
final int contextStartIndex = refPos - CONTEXT_SIZE - 1;
final int contextFullLength = 2 * CONTEXT_SIZE + 1;
if (contextStartIndex < 0 || contextStartIndex + contextFullLength > ref.length())
continue;
// skip contexts with N bases
final String context = StringUtil.bytesToString(ref.getBases(), contextStartIndex, contextFullLength).toUpperCase();
if (context.contains("N"))
continue;
// skip low BQ sites
if (failsBaseQualityCutoff(readPos, rec))
continue;
// skip N bases in read
final char readBase = Character.toUpperCase((char) rec.getReadBases()[readPos - 1]);
if (readBase == 'N')
continue;
// count the base!
artifactCounters.get(library).countRecord(context, readBase, rec);
}
}
}
use of org.broadinstitute.hellbender.exceptions.UserException in project gatk by broadinstitute.
the class CollectJumpingLibraryMetrics method doWork.
/**
* Calculates the detailed statistics about the jumping library and then generates the results.
*/
@Override
protected Object doWork() {
for (File f : INPUT) {
IOUtil.assertFileIsReadable(f);
}
IOUtil.assertFileIsWritable(OUTPUT);
Histogram<Integer> innieHistogram = new Histogram<>();
Histogram<Integer> outieHistogram = new Histogram<>();
int fragments = 0;
int innies = 0;
int outies = 0;
int innieDupes = 0;
int outieDupes = 0;
int crossChromPairs = 0;
int superSized = 0;
int tandemPairs = 0;
double chimeraSizeMinimum = Math.max(getOutieMode(), (double) CHIMERA_KB_MIN);
for (File f : INPUT) {
SamReader reader = SamReaderFactory.makeDefault().open(f);
if (reader.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate) {
throw new UserException("SAM file must " + f.getName() + " must be sorted in coordinate order");
}
for (SAMRecord sam : reader) {
// We're getting all our info from the first of each pair.
if (!sam.getFirstOfPairFlag()) {
continue;
}
// Ignore unmapped read pairs
if (sam.getReadUnmappedFlag()) {
if (!sam.getMateUnmappedFlag()) {
fragments++;
continue;
}
// If both ends are unmapped and we've hit unaligned reads we're done
if (sam.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
break;
}
continue;
}
if (sam.getMateUnmappedFlag()) {
fragments++;
continue;
}
// Ignore low-quality reads. If we don't have the mate mapping quality, assume it's OK
if ((sam.getAttribute(SAMTag.MQ.name()) != null && sam.getIntegerAttribute(SAMTag.MQ.name()) < MINIMUM_MAPPING_QUALITY) || sam.getMappingQuality() < MINIMUM_MAPPING_QUALITY) {
continue;
}
final int absInsertSize = Math.abs(sam.getInferredInsertSize());
if (absInsertSize > chimeraSizeMinimum) {
superSized++;
} else if (sam.getMateNegativeStrandFlag() == sam.getReadNegativeStrandFlag()) {
tandemPairs++;
} else if (!sam.getMateReferenceIndex().equals(sam.getReferenceIndex())) {
crossChromPairs++;
} else {
final PairOrientation pairOrientation = SamPairUtil.getPairOrientation(sam);
if (pairOrientation == PairOrientation.RF) {
outieHistogram.increment(absInsertSize);
outies++;
if (sam.getDuplicateReadFlag()) {
outieDupes++;
}
} else if (pairOrientation == PairOrientation.FR) {
innieHistogram.increment(absInsertSize);
innies++;
if (sam.getDuplicateReadFlag()) {
innieDupes++;
}
} else {
throw new IllegalStateException("This should never happen");
}
}
}
CloserUtil.close(reader);
}
MetricsFile<JumpingLibraryMetrics, Integer> metricsFile = getMetricsFile();
JumpingLibraryMetrics metrics = new JumpingLibraryMetrics();
metrics.JUMP_PAIRS = outies;
metrics.JUMP_DUPLICATE_PAIRS = outieDupes;
metrics.JUMP_DUPLICATE_PCT = outies != 0 ? outieDupes / (double) outies : 0;
metrics.JUMP_LIBRARY_SIZE = (outies > 0 && outieDupes > 0) ? DuplicationMetrics.estimateLibrarySize(outies, outies - outieDupes) : 0;
outieHistogram.trimByTailLimit(TAIL_LIMIT);
metrics.JUMP_MEAN_INSERT_SIZE = outieHistogram.getMean();
metrics.JUMP_STDEV_INSERT_SIZE = outieHistogram.getStandardDeviation();
metrics.NONJUMP_PAIRS = innies;
metrics.NONJUMP_DUPLICATE_PAIRS = innieDupes;
metrics.NONJUMP_DUPLICATE_PCT = innies != 0 ? innieDupes / (double) innies : 0;
metrics.NONJUMP_LIBRARY_SIZE = (innies > 0 && innieDupes > 0) ? DuplicationMetrics.estimateLibrarySize(innies, innies - innieDupes) : 0;
innieHistogram.trimByTailLimit(TAIL_LIMIT);
metrics.NONJUMP_MEAN_INSERT_SIZE = innieHistogram.getMean();
metrics.NONJUMP_STDEV_INSERT_SIZE = innieHistogram.getStandardDeviation();
metrics.CHIMERIC_PAIRS = crossChromPairs + superSized + tandemPairs;
metrics.FRAGMENTS = fragments;
double totalPairs = outies + innies + metrics.CHIMERIC_PAIRS;
metrics.PCT_JUMPS = totalPairs != 0 ? outies / totalPairs : 0;
metrics.PCT_NONJUMPS = totalPairs != 0 ? innies / totalPairs : 0;
metrics.PCT_CHIMERAS = totalPairs != 0 ? metrics.CHIMERIC_PAIRS / totalPairs : 0;
metricsFile.addMetric(metrics);
metricsFile.write(OUTPUT);
return null;
}
Aggregations