Examples with GenomicSequence com.github.lindenb.jvarkit.util.picard.GenomicSequence used on opensource projects

Example 21 with GenomicSequence

use of com.github.lindenb.jvarkit.util.picard.GenomicSequence in project jvarkit by lindenb.

the class IndexedGenome method _get.

private GenomicSequence _get(final String contigName) {
    final SAMSequenceRecord rec = this.dict.getSequence(contigName);
    if (rec == null)
        throw new NoSuchElementException(contigName);
    if (this.cached != null && this.cached.getChrom().equals(contigName)) {
        return this.cached;
    }
    this.cached = new GenomicSequence(this.indexedFastaSequenceFile, rec.getSequenceName());
    return this.cached;
}

Example 22 with GenomicSequence

use of com.github.lindenb.jvarkit.util.picard.GenomicSequence in project jvarkit by lindenb.

the class Biostar170742 method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.faidx == null) {
        LOG.error("Reference sequence was not defined");
        return -1;
    }
    PrintStream out = null;
    SamReader sfr = null;
    SAMRecordIterator iter = null;
    GenomicSequence genomicSequence = null;
    IndexedFastaSequenceFile indexedFastaSequenceFile = null;
    try {
        indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.faidx);
        long align_id = 0;
        sfr = openSamReader(oneFileOrNull(args));
        out = super.openFileOrStdoutAsPrintStream(this.outputFile);
        final StringBuilder refseq = new StringBuilder();
        final StringBuilder readseq = new StringBuilder();
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(sfr.getFileHeader());
        iter = sfr.iterator();
        while (iter.hasNext()) {
            final SAMRecord rec = progress.watch(iter.next());
            if (rec.getReadUnmappedFlag())
                continue;
            final Cigar cigar = rec.getCigar();
            if (cigar == null)
                continue;
            final byte[] readbases = rec.getReadBases();
            if (readbases == null)
                continue;
            if (genomicSequence == null || !rec.getReferenceName().equals(genomicSequence.getChrom())) {
                genomicSequence = new GenomicSequence(indexedFastaSequenceFile, rec.getReferenceName());
            }
            int refpos1 = rec.getAlignmentStart();
            int readpos = 0;
            refseq.setLength(0);
            readseq.setLength(0);
            for (final CigarElement ce : cigar.getCigarElements()) {
                final CigarOperator op = ce.getOperator();
                if (op.equals(CigarOperator.S)) {
                    readpos += ce.getLength();
                    continue;
                }
                if (op.equals(CigarOperator.H)) {
                    continue;
                }
                for (int i = 0; i < ce.getLength(); ++i) {
                    if (op.consumesReferenceBases() && op.consumesReadBases()) {
                        refseq.append(genomicSequence.charAt(refpos1 - 1));
                        readseq.append((char) readbases[readpos]);
                        readpos++;
                        refpos1++;
                    } else if (op.consumesReferenceBases()) {
                        refseq.append(genomicSequence.charAt(refpos1 - 1));
                        readseq.append('-');
                        refpos1++;
                    } else if (op.consumesReadBases()) {
                        refseq.append('-');
                        readseq.append((char) readbases[readpos]);
                        readpos++;
                    }
                }
            }
            out.print(align_id);
            out.print(' ');
            out.print(rec.getReferenceName());
            out.print(' ');
            out.print(rec.getAlignmentStart());
            out.print(' ');
            out.print(rec.getAlignmentEnd());
            out.print(' ');
            out.print(rec.getReadName());
            if (rec.getReadPairedFlag()) {
                if (rec.getFirstOfPairFlag()) {
                    out.print("/1");
                } else if (rec.getSecondOfPairFlag()) {
                    out.print("/2");
                }
            }
            out.print(' ');
            out.print(1 + rec.getAlignmentStart() - rec.getUnclippedStart());
            out.print(' ');
            out.print(rec.getReadLength() - (rec.getUnclippedEnd() - rec.getAlignmentEnd()));
            out.print(' ');
            out.print(rec.getReadNegativeStrandFlag() ? "-" : "+");
            out.print(' ');
            out.print(rec.getMappingQuality());
            out.println();
            out.println(refseq);
            out.println(readseq);
            out.println();
            ++align_id;
        }
        progress.finish();
        iter.close();
        out.flush();
        LOG.info("done");
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(out);
        CloserUtil.close(iter);
        CloserUtil.close(sfr);
        CloserUtil.close(indexedFastaSequenceFile);
    }
}

Example 23 with GenomicSequence

use of com.github.lindenb.jvarkit.util.picard.GenomicSequence 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);
    }
}

Example 24 with GenomicSequence

use of com.github.lindenb.jvarkit.util.picard.GenomicSequence in project jvarkit by lindenb.

the class BimToVcf method doWork.

@Override
public int doWork(List<String> args) {
    VariantContextWriter w = null;
    BufferedReader r = null;
    IndexedFastaSequenceFile faidx = null;
    GenomicSequence genomic = null;
    try {
        if (this.REF == null) {
            LOG.error("Reference -R missing.");
            return -1;
        }
        faidx = new IndexedFastaSequenceFile(this.REF);
        final SAMSequenceDictionary dict = faidx.getSequenceDictionary();
        if (dict == null) {
            LOG.error("No dictionary in " + this.REF);
            return -1;
        }
        r = super.openBufferedReader(oneFileOrNull(args));
        final Set<VCFHeaderLine> headerLines = new HashSet<>();
        final VCFInfoHeaderLine morgan = new VCFInfoHeaderLine("MORGAN", 1, VCFHeaderLineType.Float, "Centimorgan");
        final VCFInfoHeaderLine svtype = new VCFInfoHeaderLine("SVTYPE", 1, VCFHeaderLineType.String, "Variation type");
        VCFStandardHeaderLines.addStandardInfoLines(headerLines, false, "");
        // super.addMetaData(headerLines);
        headerLines.add(morgan);
        headerLines.add(svtype);
        final List<String> genotypeSampleNames = Collections.emptyList();
        final VCFHeader header = new VCFHeader(headerLines, genotypeSampleNames);
        header.setSequenceDictionary(dict);
        w = super.openVariantContextWriter(this.outputFile);
        w.writeHeader(header);
        final Pattern tab = Pattern.compile("[\t]");
        String line;
        final Pattern iupacATGC = Pattern.compile("[atgcATGC]");
        while ((line = r.readLine()) != null) {
            String[] tokens = tab.split(line);
            if (tokens.length != 6) {
                LOG.error("expected 6 column in " + line);
                return -1;
            }
            Allele a1 = null;
            Allele a2 = null;
            Allele ref = null;
            String contig = tokens[0];
            SAMSequenceRecord ssr = null;
            ssr = dict.getSequence(contig);
            // ugly below !!
            if (ssr == null && contig.equals("23")) {
                ssr = dict.getSequence("X");
            }
            if (ssr == null && contig.equals("23")) {
                ssr = dict.getSequence("chrX");
            }
            if (ssr == null && contig.equals("24")) {
                ssr = dict.getSequence("Y");
            }
            if (ssr == null && contig.equals("24")) {
                ssr = dict.getSequence("chrY");
            }
            if (ssr == null && contig.equals("26")) {
                ssr = dict.getSequence("chrM");
            }
            if (ssr == null && contig.equals("26")) {
                ssr = dict.getSequence("MT");
            }
            if (ssr == null && contig.equals("25")) {
                LOG.warn("ignoring " + line);
                continue;
            }
            if (ssr == null) {
                LOG.error("unknown chrom in " + line);
                return -1;
            }
            if (genomic == null || !ssr.getSequenceName().equals(genomic.getChrom())) {
                genomic = new GenomicSequence(faidx, ssr.getSequenceName());
            }
            int pos1 = Integer.parseInt(tokens[3]);
            if (tokens[4].equals("0"))
                tokens[4] = tokens[5];
            if (tokens[5].equals("0"))
                tokens[5] = tokens[4];
            final VariantContextBuilder vcb = new VariantContextBuilder();
            vcb.chr(ssr.getSequenceName());
            vcb.attribute(morgan.getID(), Float.parseFloat(tokens[2]));
            if (iupacATGC.matcher(tokens[4]).matches() && iupacATGC.matcher(tokens[5]).matches()) {
                String refBase = String.valueOf(genomic.charAt(pos1 - 1));
                ref = Allele.create(refBase, true);
                a1 = refBase.equalsIgnoreCase(tokens[4]) ? ref : Allele.create(tokens[4], false);
                a2 = refBase.equalsIgnoreCase(tokens[5]) ? ref : Allele.create(tokens[5], false);
                vcb.attribute(svtype.getID(), a1.isReference() && a2.isReference() ? "NOVARIATION" : "SNV");
            } else if ((tokens[4].equals("-") && iupacATGC.matcher(tokens[5]).matches()) || (tokens[5].equals("-") && iupacATGC.matcher(tokens[4]).matches())) {
                // shift left
                pos1--;
                String refBase = String.valueOf(genomic.charAt(pos1 - 1));
                a1 = Allele.create(refBase, false);
                ref = Allele.create(refBase + tokens[tokens[4].equals("-") ? 5 : 4], true);
                a2 = a1;
                vcb.attribute(svtype.getID(), "DEL");
            } else if (tokens[4].equals("-") && tokens[5].equals("-")) {
                // shift left
                pos1--;
                String refBase = String.valueOf(genomic.charAt(pos1 - 1));
                a1 = Allele.create(refBase, false);
                ref = Allele.create(refBase + genomic.charAt(pos1), true);
                a2 = a1;
                vcb.attribute(svtype.getID(), "DEL");
            } else {
                LOG.error("not handled: " + line);
                return -1;
            }
            final Set<Allele> alleles = new HashSet<>();
            alleles.add(ref);
            alleles.add(a1);
            alleles.add(a2);
            vcb.start(pos1);
            vcb.stop(pos1 + ref.length() - 1);
            if (!tokens[1].isEmpty())
                vcb.id(tokens[1]);
            vcb.alleles(alleles);
            w.add(vcb.make());
        }
        r.close();
        r = null;
        w.close();
        w = null;
        return RETURN_OK;
    } catch (final Exception e) {
        LOG.error(e);
        return -1;
    } finally {
        CloserUtil.close(faidx);
        CloserUtil.close(w);
        CloserUtil.close(r);
    }
}