use of com.github.lindenb.jvarkit.util.iterator.FilterIterator in project jvarkit by lindenb.
the class SamTranslocations method doWork.
@Override
public int doWork(final List<String> args) {
if (this.max_distance <= 0) {
LOG.error("max_distance <=0 (" + max_distance + ")");
return -1;
}
PrintWriter pw = null;
ConcatSam.ConcatSamIterator samIter = null;
ForwardPeekIterator<SAMRecord> forwardPeekIterator = null;
try {
samIter = new ConcatSam.Factory().addInterval(this.region_str).open(args);
final SAMSequenceDictionary dict = samIter.getFileHeader().getSequenceDictionary();
if (dict.size() < 2) {
LOG.error("Not enough contigs in sequence dictionary. Expected at least 2.");
return -1;
}
forwardPeekIterator = new ForwardPeekIteratorImpl<>(new FilterIterator<>(samIter, rec -> {
if (rec.getReadUnmappedFlag())
return false;
if (!rec.getReadPairedFlag())
return false;
if (rec.getMateUnmappedFlag())
return false;
final int tid1 = rec.getReferenceIndex();
final int tid2 = rec.getMateReferenceIndex();
if (tid1 == tid2)
return false;
if (this.samRecordFilter.filterOut(rec))
return false;
return true;
}));
pw = openFileOrStdoutAsPrintWriter(this.outputFile);
pw.print("#chrom1");
pw.print('\t');
pw.print("chrom1-start");
pw.print('\t');
pw.print("chrom1-end");
pw.print('\t');
pw.print("middle1\tstrand_plus_before_mid1_count\tstrand_plus_before_mid1_percent\tstrand_minus_after_mid1_count\tstrand_minus_after_mid1_percent");
pw.print('\t');
pw.print("chrom2");
pw.print('\t');
pw.print("chrom2-start");
pw.print('\t');
pw.print("chrom2-end");
pw.print('\t');
pw.print("middle2\tstrand_plus_before_mid2_count\tstrand_plus_before_mid2_percent\tstrand_minus_after_mid2_count\tstrand_minus_after_mid2_percent");
pw.print('\t');
pw.print("count-reads");
pw.print('\t');
pw.print("count-clipped");
pw.print('\t');
pw.print("partition");
pw.println();
final Map<String, PartitionState> partition2state = new HashMap<>();
final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict).logger(LOG);
while (forwardPeekIterator.hasNext()) {
final SAMRecord rec = progress.watch(forwardPeekIterator.next());
final String partition = this.samRecordPartition.getPartion(rec, "N/A");
PartitionState partitionState = partition2state.get(partition);
if (partitionState == null) {
partitionState = new PartitionState(partition);
partition2state.put(partition, partitionState);
}
if (partitionState.last_rec != null) {
if (partitionState.last_rec == rec) {
// reset last
partitionState.last_rec = null;
}
continue;
}
// searching for -->
if (rec.getReadNegativeStrandFlag())
continue;
// find forward records
final List<SAMRecord> positiveStrandList = new ArrayList<>();
positiveStrandList.add(rec);
int x = 0;
for (; ; ) {
final SAMRecord rec2 = forwardPeekIterator.peek(x);
if (rec2 == null) {
break;
}
if (!rec.getReferenceIndex().equals(rec2.getReferenceIndex())) {
break;
}
if (rec2.getStart() - rec.getEnd() > this.max_distance) {
break;
}
if (rec2.getReadNegativeStrandFlag() || Math.abs(rec2.getEnd() - rec.getEnd()) > this.fully_distance || !rec2.getMateReferenceName().equals(rec.getMateReferenceName())) {
++x;
continue;
}
positiveStrandList.add(rec2);
partitionState.last_rec = rec2;
++x;
}
if (positiveStrandList.size() < this.min_count_forward) {
partitionState.last_rec = null;
continue;
}
// find reverse records
x = 0;
final List<SAMRecord> negativeStrandList = new ArrayList<>();
for (; ; ) {
final SAMRecord rec2 = forwardPeekIterator.peek(x);
if (rec2 == null) {
break;
}
if (!rec.getReferenceIndex().equals(rec2.getReferenceIndex())) {
break;
}
if (rec2.getStart() - rec.getEnd() > this.max_distance) {
break;
}
if (!rec2.getReadNegativeStrandFlag() || (rec2.getStart() < rec.getEnd() && (rec.getEnd() - rec2.getStart()) > this.fully_distance) || !rec2.getMateReferenceName().equals(rec.getMateReferenceName())) {
++x;
continue;
}
negativeStrandList.add(rec2);
++x;
}
if (negativeStrandList.size() < this.min_count_forward) {
partitionState.last_rec = null;
continue;
}
final int start1 = (int) Percentile.median().evaluate(positiveStrandList.stream().mapToInt(R -> R.getEnd()));
final int end1 = (int) Percentile.median().evaluate(negativeStrandList.stream().mapToInt(R -> R.getStart()));
final int mid1 = (start1 + end1) / 2;
final int mid2 = (int) Percentile.median().evaluate(negativeStrandList.stream().mapToInt(R -> R.getMateAlignmentStart()));
pw.print(rec.getReferenceName());
pw.print('\t');
pw.print(positiveStrandList.stream().mapToInt(R -> R.getEnd()).max().getAsInt());
pw.print('\t');
pw.print(negativeStrandList.stream().mapToInt(R -> R.getStart()).min().getAsInt());
pw.print('\t');
pw.print(mid1);
pw.print('\t');
pw.print(positiveStrandList.stream().filter(R -> R.getEnd() <= mid1).count());
pw.print('\t');
pw.print((int) (100.0 * (positiveStrandList.stream().filter(R -> R.getEnd() <= mid1).count()) / positiveStrandList.size()));
pw.print('\t');
pw.print(negativeStrandList.stream().filter(R -> R.getStart() >= mid1).count());
pw.print('\t');
pw.print((int) (100.0 * (negativeStrandList.stream().filter(R -> R.getStart() >= mid1).count()) / positiveStrandList.size()));
pw.print('\t');
pw.print(rec.getMateReferenceName());
pw.print('\t');
pw.print(negativeStrandList.stream().mapToInt(R -> R.getMateAlignmentStart()).min().getAsInt());
pw.print('\t');
pw.print(negativeStrandList.stream().mapToInt(R -> R.getMateAlignmentStart()).max().getAsInt());
pw.print('\t');
pw.print(mid2);
pw.print('\t');
pw.print(negativeStrandList.stream().filter(R -> R.getMateAlignmentStart() <= mid2).count());
pw.print('\t');
pw.print((int) (100.0 * (negativeStrandList.stream().filter(R -> R.getMateAlignmentStart() <= mid2).count()) / negativeStrandList.size()));
pw.print('\t');
pw.print(negativeStrandList.stream().filter(R -> R.getMateAlignmentStart() >= mid2).count());
pw.print('\t');
pw.print((int) (100.0 * (negativeStrandList.stream().filter(R -> R.getMateAlignmentStart() >= mid2).count()) / negativeStrandList.size()));
pw.print('\t');
pw.print(positiveStrandList.size() + negativeStrandList.size());
pw.print('\t');
pw.print(positiveStrandList.stream().filter(R -> R.getCigar() != null && R.getCigar().isClipped()).count() + negativeStrandList.stream().filter(R -> R.getCigar() != null && R.getCigar().isClipped()).count());
pw.print('\t');
pw.print(partition);
pw.println();
}
progress.finish();
forwardPeekIterator.close();
forwardPeekIterator = null;
samIter.close();
samIter = null;
pw.flush();
pw.close();
pw = null;
return 0;
} catch (final Throwable err) {
LOG.error(err);
return -1;
} finally {
CloserUtil.close(forwardPeekIterator);
CloserUtil.close(samIter);
CloserUtil.close(pw);
}
}
use of com.github.lindenb.jvarkit.util.iterator.FilterIterator in project jvarkit by lindenb.
the class Biostar78285 method doWork.
@Override
public int doWork(final List<String> args) {
if (this.gc_percent_window < 1) {
LOG.error("Bad GC% window size:" + this.gc_percent_window);
return -1;
}
final List<File> bamFiles = IOUtil.unrollFiles(args.stream().map(F -> new File(F)).collect(Collectors.toCollection(HashSet::new)), ".bam");
SAMSequenceDictionary dict = null;
final List<SamReader> samReaders = new ArrayList<>();
final List<CloseableIterator<SAMRecord>> samIterators = new ArrayList<>();
final TreeSet<String> samples = new TreeSet<>();
final String DEFAULT_PARTITION = "UNDEFINED_PARTITION";
IndexedFastaSequenceFile indexedFastaSequenceFile = null;
VariantContextWriter out = null;
try {
final SamReaderFactory samReaderFactory = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
for (final File bamFile : bamFiles) {
LOG.info("Opening " + bamFile);
final SamReader samReader = samReaderFactory.open(bamFile);
samReaders.add(samReader);
final SAMFileHeader header = samReader.getFileHeader();
if (header == null) {
LOG.error("No header in " + bamFile);
return -1;
}
if (header.getSortOrder() != SortOrder.coordinate) {
LOG.error("Sam file " + bamFile + " is not sorted on coordinate :" + header.getSortOrder());
return -1;
}
samples.addAll(header.getReadGroups().stream().map(RG -> this.partition.apply(RG, DEFAULT_PARTITION)).collect(Collectors.toSet()));
final SAMSequenceDictionary currDict = header.getSequenceDictionary();
if (currDict == null) {
LOG.error("SamFile doesn't contain a SAMSequenceDictionary : " + bamFile);
return -1;
}
if (dict == null) {
dict = currDict;
} else if (!SequenceUtil.areSequenceDictionariesEqual(dict, currDict)) {
LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict, currDict));
return -1;
}
}
if (samReaders.isEmpty()) {
LOG.error("no bam");
return -1;
}
if (dict == null) {
LOG.error("no dictionary");
return -1;
}
final QueryInterval[] intervals;
if (this.captureBed != null) {
LOG.info("Opening " + this.captureBed);
ContigNameConverter.setDefaultAliases(dict);
final List<QueryInterval> L = new ArrayList<>();
final BedLineCodec codec = new BedLineCodec();
final LineIterator li = IOUtils.openFileForLineIterator(this.captureBed);
while (li.hasNext()) {
final BedLine bed = codec.decode(li.next());
if (bed == null)
continue;
final QueryInterval q = bed.toQueryInterval(dict);
L.add(q);
}
CloserUtil.close(li);
intervals = QueryInterval.optimizeIntervals(L.toArray(new QueryInterval[L.size()]));
} else {
intervals = null;
}
for (final SamReader samReader : samReaders) {
LOG.info("querying " + samReader.getResourceDescription());
final CloseableIterator<SAMRecord> iter;
if (intervals == null) {
iter = samReader.iterator();
} else {
iter = samReader.queryOverlapping(intervals);
}
samIterators.add(new FilterIterator<SAMRecord>(iter, R -> !R.getReadUnmappedFlag() && !filter.filterOut(R)));
}
if (this.refFile != null) {
LOG.info("opening " + refFile);
indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.refFile);
final SAMSequenceDictionary refdict = indexedFastaSequenceFile.getSequenceDictionary();
ContigNameConverter.setDefaultAliases(refdict);
if (refdict == null) {
throw new JvarkitException.FastaDictionaryMissing(this.refFile);
}
if (!SequenceUtil.areSequenceDictionariesEqual(dict, refdict)) {
LOG.error(JvarkitException.DictionariesAreNotTheSame.getMessage(dict, refdict));
return -1;
}
}
out = openVariantContextWriter(this.outputFile);
final Set<VCFHeaderLine> metaData = new HashSet<>();
VCFStandardHeaderLines.addStandardFormatLines(metaData, true, VCFConstants.DEPTH_KEY, VCFConstants.GENOTYPE_KEY);
VCFStandardHeaderLines.addStandardInfoLines(metaData, true, VCFConstants.DEPTH_KEY);
metaData.add(new VCFFormatHeaderLine("DF", 1, VCFHeaderLineType.Integer, "Number of Reads on plus strand"));
metaData.add(new VCFFormatHeaderLine("DR", 1, VCFHeaderLineType.Integer, "Number of Reads on minus strand"));
metaData.add(new VCFInfoHeaderLine("AVG_DP", 1, VCFHeaderLineType.Float, "Mean depth"));
metaData.add(new VCFInfoHeaderLine("MEDIAN_DP", 1, VCFHeaderLineType.Float, "Median depth"));
metaData.add(new VCFInfoHeaderLine("MIN_DP", 1, VCFHeaderLineType.Integer, "Min depth"));
metaData.add(new VCFInfoHeaderLine("MAX_DP", 1, VCFHeaderLineType.Integer, "Max depth"));
metaData.add(new VCFHeaderLine(Biostar78285.class.getSimpleName() + ".SamFilter", this.filter.toString()));
for (final Integer treshold : this.minDepthTresholds) {
metaData.add(new VCFFilterHeaderLine("DP_LT_" + treshold, "All genotypes have DP< " + treshold));
metaData.add(new VCFInfoHeaderLine("NUM_DP_LT_" + treshold, 1, VCFHeaderLineType.Integer, "Number of genotypes having DP< " + treshold));
metaData.add(new VCFInfoHeaderLine("FRACT_DP_LT_" + treshold, 1, VCFHeaderLineType.Float, "Fraction of genotypes having DP< " + treshold));
}
if (indexedFastaSequenceFile != null) {
metaData.add(new VCFInfoHeaderLine("GC_PERCENT", 1, VCFHeaderLineType.Integer, "GC% window_size:" + this.gc_percent_window));
}
final List<Allele> refAlleles = Collections.singletonList(Allele.create("N", true));
final List<Allele> NO_CALLS = Arrays.asList(Allele.NO_CALL, Allele.NO_CALL);
final VCFHeader vcfHeader = new VCFHeader(metaData, samples);
vcfHeader.setSequenceDictionary(dict);
out.writeHeader(vcfHeader);
final SAMRecordCoordinateComparator samRecordCoordinateComparator = new SAMRecordCoordinateComparator();
final PeekableIterator<SAMRecord> peekIter = new PeekableIterator<>(new MergingIterator<>((R1, R2) -> samRecordCoordinateComparator.fileOrderCompare(R1, R2), samIterators));
final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict);
for (final SAMSequenceRecord ssr : dict.getSequences()) {
final IntervalTree<Boolean> capturePos;
if (intervals != null) {
if (!Arrays.stream(intervals).anyMatch(I -> I.referenceIndex == ssr.getSequenceIndex())) {
continue;
}
capturePos = new IntervalTree<>();
Arrays.stream(intervals).filter(I -> I.referenceIndex == ssr.getSequenceIndex()).forEach(I -> capturePos.put(I.start, I.end, true));
;
} else {
capturePos = null;
}
final GenomicSequence genomicSequence;
if (indexedFastaSequenceFile != null && indexedFastaSequenceFile.getSequenceDictionary().getSequence(ssr.getSequenceName()) != null) {
genomicSequence = new GenomicSequence(indexedFastaSequenceFile, ssr.getSequenceName());
} else {
genomicSequence = null;
}
final List<SAMRecord> buffer = new ArrayList<>();
for (int ssr_pos = 1; ssr_pos <= ssr.getSequenceLength(); ++ssr_pos) {
if (capturePos != null && !capturePos.overlappers(ssr_pos, ssr_pos).hasNext())
continue;
progress.watch(ssr.getSequenceName(), ssr_pos);
while (peekIter.hasNext()) {
final SAMRecord rec = peekIter.peek();
if (rec.getReadUnmappedFlag()) {
// consumme
peekIter.next();
continue;
}
if (this.filter.filterOut(rec)) {
// consumme
peekIter.next();
continue;
}
if (rec.getReferenceIndex() < ssr.getSequenceIndex()) {
throw new IllegalStateException("should not happen");
}
if (rec.getReferenceIndex() > ssr.getSequenceIndex()) {
break;
}
if (rec.getAlignmentEnd() < ssr_pos) {
throw new IllegalStateException("should not happen");
}
if (rec.getAlignmentStart() > ssr_pos) {
break;
}
buffer.add(peekIter.next());
}
int x = 0;
while (x < buffer.size()) {
final SAMRecord R = buffer.get(x);
if (R.getReferenceIndex() != ssr.getSequenceIndex() || R.getAlignmentEnd() < ssr_pos) {
buffer.remove(x);
} else {
x++;
}
}
final Map<String, PosInfo> count = samples.stream().map(S -> new PosInfo(S)).collect(Collectors.toMap(P -> P.sample, Function.identity()));
for (final SAMRecord rec : buffer) {
if (rec.getReferenceIndex() != ssr.getSequenceIndex())
throw new IllegalStateException("should not happen");
if (rec.getAlignmentEnd() < ssr_pos)
continue;
if (rec.getAlignmentStart() > ssr_pos)
continue;
final Cigar cigar = rec.getCigar();
if (cigar == null)
continue;
int refpos = rec.getAlignmentStart();
final String sample = this.partition.getPartion(rec, DEFAULT_PARTITION);
for (final CigarElement ce : cigar.getCigarElements()) {
if (refpos > ssr_pos)
break;
final CigarOperator op = ce.getOperator();
if (op.consumesReferenceBases()) {
if (op.consumesReadBases()) {
if (refpos <= ssr_pos && ssr_pos <= refpos + ce.getLength()) {
final PosInfo posInfo = count.get(sample);
if (posInfo != null) {
posInfo.dp++;
if (rec.getReadNegativeStrandFlag()) {
posInfo.negative_strand++;
}
}
break;
}
}
refpos += ce.getLength();
}
}
}
final VariantContextBuilder vcb = new VariantContextBuilder();
final Set<String> filters = new HashSet<>();
vcb.chr(ssr.getSequenceName());
vcb.start(ssr_pos);
vcb.stop(ssr_pos);
if (genomicSequence == null) {
vcb.alleles(refAlleles);
} else {
vcb.alleles(Collections.singletonList(Allele.create((byte) genomicSequence.charAt(ssr_pos - 1), true)));
final GenomicSequence.GCPercent gcp = genomicSequence.getGCPercent(Math.max((ssr_pos - 1) - this.gc_percent_window, 0), Math.min(ssr_pos + this.gc_percent_window, ssr.getSequenceLength()));
if (!gcp.isEmpty()) {
vcb.attribute("GC_PERCENT", gcp.getGCPercentAsInteger());
}
}
vcb.attribute(VCFConstants.DEPTH_KEY, (int) count.values().stream().mapToInt(S -> S.dp).sum());
vcb.genotypes(count.values().stream().map(C -> new GenotypeBuilder(C.sample, NO_CALLS).DP((int) C.dp).attribute("DR", C.negative_strand).attribute("DF", C.dp - C.negative_strand).make()).collect(Collectors.toList()));
for (final Integer treshold : this.minDepthTresholds) {
final int count_lt = (int) count.values().stream().filter(S -> S.dp < treshold).count();
if (count_lt == samples.size()) {
filters.add("DP_LT_" + treshold);
}
vcb.attribute("NUM_DP_LT_" + treshold, count_lt);
if (!samples.isEmpty()) {
vcb.attribute("FRACT_DP_LT_" + treshold, count_lt / (float) samples.size());
}
}
if (!samples.isEmpty()) {
final int[] array = count.values().stream().mapToInt(S -> S.dp).toArray();
vcb.attribute("AVG_DP", Percentile.average().evaluate(array));
vcb.attribute("MEDIAN_DP", Percentile.median().evaluate(array));
vcb.attribute("MIN_DP", (int) Percentile.min().evaluate(array));
vcb.attribute("MAX_DP", (int) Percentile.max().evaluate(array));
}
if (filters.isEmpty()) {
vcb.passFilters();
} else {
vcb.filters(filters);
}
out.add(vcb.make());
}
}
progress.finish();
peekIter.close();
out.close();
out = null;
return 0;
} catch (final Exception err) {
LOG.error(err);
return -1;
} finally {
CloserUtil.close(out);
CloserUtil.close(samIterators);
CloserUtil.close(samReaders);
CloserUtil.close(indexedFastaSequenceFile);
}
}
Aggregations