Search in sources :

Example 41 with CigarOperator

use of htsjdk.samtools.CigarOperator in project jvarkit by lindenb.

the class TView method paint.

void paint(final PrintStream out) {
    final Colorizer colorizer;
    switch(this.formatOut) {
        case html:
            colorizer = new HtmlColorizer(out);
            break;
        case tty:
            colorizer = new AnsiColorizer(out);
            break;
        case plain:
            colorizer = new Colorizer(out);
            break;
        default:
            throw new IllegalStateException();
    }
    if (interval == null) {
        LOG.warn("No interval defined");
        return;
    }
    final GenomicSequence contigSequence;
    final Function<Integer, Character> refPosToBase;
    if (indexedFastaSequenceFile != null) {
        final SAMSequenceDictionary dict = SAMSequenceDictionaryExtractor.extractDictionary(referenceFile);
        if (dict.getSequence(this.interval.getContig()) == null) {
            LOG.warn("No interval with contig " + interval + " in REF");
            return;
        }
        contigSequence = new GenomicSequence(indexedFastaSequenceFile, interval.getContig());
        refPosToBase = POS -> {
            if (POS < 0 || POS >= contigSequence.length())
                return 'N';
            return contigSequence.charAt(POS);
        };
    } else {
        contigSequence = null;
        refPosToBase = POS -> 'N';
    }
    /**
     * test if genomic position is in interval
     */
    final Predicate<Integer> testInInterval = new Predicate<Integer>() {

        @Override
        public boolean test(final Integer pos) {
            return interval.getStart() <= pos && pos <= interval.getEnd();
        }
    };
    final int pixelWidth = this.interval.length();
    final Map<Integer, Integer> genomicpos2insertlen = new TreeMap<>();
    final Map<String, List<SAMRecord>> group2record = new TreeMap<>();
    for (final SamReader samReader : this.samReaders) {
        SAMRecordIterator iter = samReader.query(this.interval.getContig(), this.interval.getStart(), this.interval.getEnd(), false);
        while (iter.hasNext()) {
            final SAMRecord rec = iter.next();
            if (rec.getReadUnmappedFlag())
                continue;
            if (rec.getCigar() == null)
                continue;
            if (getRecordFilter().filterOut(rec))
                continue;
            if (!rec.getContig().equals(interval.getContig()))
                continue;
            if (right().apply(rec) < this.interval.getStart())
                continue;
            if (this.interval.getEnd() < left().apply(rec))
                continue;
            String group = this.groupBy.getPartion(rec);
            if (group == null || group.isEmpty()) {
                group = "undefined_" + this.groupBy.name();
            }
            List<SAMRecord> records = group2record.get(group);
            if (records == null) {
                records = new ArrayList<>();
                group2record.put(group, records);
            }
            records.add(rec);
            // loop over cigar, get the longest insert
            int refpos = rec.getAlignmentStart();
            for (final CigarElement ce : rec.getCigar().getCigarElements()) {
                if (!this.showInsertions)
                    break;
                final CigarOperator op = ce.getOperator();
                if (op.equals(CigarOperator.I) && testInInterval.test(refpos)) {
                    final Integer longestInsert = genomicpos2insertlen.get(refpos);
                    if (longestInsert == null || longestInsert.compareTo(ce.getLength()) < 0) {
                        genomicpos2insertlen.put(refpos, ce.getLength());
                    }
                }
                if (op.consumesReferenceBases()) {
                    refpos += ce.getLength();
                }
                if (refpos > interval.getEnd())
                    break;
            }
        }
        CloserUtil.close(iter);
        CloserUtil.close(samReader);
    }
    /**
     * compute where are the insertions
     */
    // LOG.debug(genomicpos2insertlen);
    final Predicate<Integer> insertIsPresentAtX = SCREENX -> {
        int x = 0;
        int ref = interval.getStart();
        while (x < pixelWidth) {
            if (x > SCREENX)
                return false;
            final Integer insertLen = genomicpos2insertlen.get(ref);
            if (insertLen == null) {
                ++x;
                ++ref;
            } else {
                if (x <= SCREENX && SCREENX < x + insertLen)
                    return true;
                // (+1) I DON'T UNDERSTAND WHY, BUT IT WORKS
                x += (insertLen + 1);
                ++ref;
            }
        }
        return false;
    };
    final Function<Character, AnsiColor> base2ansiColor = BASE -> {
        switch(Character.toUpperCase(BASE)) {
            case 'A':
                return AnsiColor.BLUE;
            case 'T':
                return AnsiColor.GREEN;
            case 'G':
                return AnsiColor.CYAN;
            case 'C':
                return AnsiColor.YELLOW;
            default:
                return null;
        }
    };
    /**
     * print interval title
     */
    out.println(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd());
    /**
     * paint base position
     */
    int ref = this.interval.getStart();
    int x = 0;
    out.print(margin("POS:"));
    while (x < pixelWidth) {
        if (insertIsPresentAtX.test(x)) {
            colorizer.pen(AnsiColor.RED).print("^");
            ++x;
        } else if ((ref - this.interval.getStart()) % 10 == 0) {
            final String f = String.format("%d", ref);
            for (int i = 0; i < f.length() && x < pixelWidth; ++i) {
                colorizer.pen(AnsiColor.GREEN).print(f.charAt(i));
                if (!insertIsPresentAtX.test(x))
                    ++ref;
                ++x;
            }
        } else {
            out.print(".");
            ++ref;
            ++x;
        }
    }
    out.println();
    /* paint ref base */
    out.print(margin("REF:"));
    ref = this.interval.getStart();
    x = 0;
    while (x < pixelWidth) {
        if (insertIsPresentAtX.test(x)) {
            colorizer.paper(AnsiColor.YELLOW).print("*");
            ++x;
        } else {
            char refBase = refPosToBase.apply(ref - 1);
            colorizer.pen(base2ansiColor.apply(refBase)).print(refBase);
            ++ref;
            ++x;
        }
    }
    out.println();
    /* loop over samples **/
    for (final String groupName : group2record.keySet()) {
        if (this.maxReadRowPerGroup == 0)
            continue;
        final ConsensusBuilder consensus = new ConsensusBuilder();
        int y_group = 0;
        final List<List<SAMRecord>> rows = new ArrayList<>();
        out.println(margin(""));
        switch(this.layoutReads) {
            case name:
                {
                    rows.addAll(group2record.get(groupName).stream().sorted((R1, R2) -> R1.getReadName().compareTo(R2.getReadName())).map(R -> Collections.singletonList(R)).collect(Collectors.toList()));
                    break;
                }
            default:
                {
                    /* pileup reads */
                    for (final SAMRecord rec : group2record.get(groupName)) {
                        int y = 0;
                        for (y = 0; y < rows.size(); ++y) {
                            final List<SAMRecord> row = rows.get(y);
                            final SAMRecord last = row.get(row.size() - 1);
                            if (right().apply(last) + this.distance_between_reads < left().apply(rec)) {
                                row.add(rec);
                                break;
                            }
                        }
                        if (y == rows.size()) {
                            final List<SAMRecord> row = new ArrayList<>();
                            row.add(rec);
                            rows.add(row);
                        }
                    }
                    break;
                }
        }
        // each row is only one read, so we need to print the groupName
        if (layoutReads == LayoutReads.name) {
            out.print(margin(groupName));
            out.println();
        }
        /* print each row */
        for (final List<SAMRecord> row : rows) {
            ++y_group;
            boolean print_this_line = (this.maxReadRowPerGroup < 0 || y_group <= this.maxReadRowPerGroup);
            if (print_this_line) {
                // each row is only one read, print the read name
                if (layoutReads == LayoutReads.name) {
                    String readName = row.get(0).getReadName();
                    if (row.get(0).getReadPairedFlag()) {
                        if (row.get(0).getFirstOfPairFlag()) {
                            readName += "/1";
                        }
                        if (row.get(0).getSecondOfPairFlag()) {
                            readName += "/2";
                        }
                    }
                    out.print(margin(readName));
                } else {
                    out.print(margin(y_group == 1 ? groupName : ""));
                }
            }
            ref = interval.getStart();
            x = 0;
            for (final SAMRecord rec : row) {
                int readRef = left().apply(rec);
                // pad before record
                while (x < pixelWidth && ref < readRef && testInInterval.test(ref)) {
                    if (!insertIsPresentAtX.test(x))
                        ++ref;
                    ++x;
                    if (print_this_line)
                        out.print(' ');
                    consensus.add(' ');
                }
                int readpos = 0;
                /* get read base function */
                final Function<Integer, Character> baseAt = new Function<Integer, Character>() {

                    @Override
                    public Character apply(final Integer readpos) {
                        final byte[] readBases = rec.getReadBases();
                        if (readBases == SAMRecord.NULL_SEQUENCE)
                            return 'N';
                        if (readpos < 0 || readpos >= rec.getReadLength())
                            return '?';
                        return (char) readBases[readpos];
                    }
                };
                for (final CigarElement ce : rec.getCigar()) {
                    final CigarOperator op = ce.getOperator();
                    if (op.equals(CigarOperator.PADDING))
                        continue;
                    /* IN INSERTION, only print if showInsertions is true */
                    if (this.showInsertions && op.equals(CigarOperator.I)) {
                        int cigarIdx = 0;
                        while (x < pixelWidth && cigarIdx < ce.getLength()) {
                            if (testInInterval.test(readRef)) {
                                final char readbase = baseAt.apply(readpos);
                                if (print_this_line)
                                    colorizer.paper(AnsiColor.RED).print(readbase);
                                consensus.add(readbase);
                                ++x;
                            }
                            ++cigarIdx;
                            ++readpos;
                        }
                        continue;
                    }
                    int cigarIdx = 0;
                    while (x < pixelWidth && cigarIdx < ce.getLength()) {
                        colorizer.clear();
                        // pad before base
                        while (x < pixelWidth && testInInterval.test(readRef) && (insertIsPresentAtX.test(x))) {
                            ++x;
                            if (print_this_line)
                                colorizer.paper(AnsiColor.YELLOW).print("*");
                            consensus.add(' ');
                            continue;
                        }
                        switch(op) {
                            case I:
                                {
                                    // if visible, processed above
                                    if (showInsertions)
                                        throw new IllegalStateException();
                                    readpos++;
                                    break;
                                }
                            case P:
                                break;
                            case H:
                                {
                                    if (showClip) {
                                        if (testInInterval.test(readRef)) {
                                            if (print_this_line)
                                                colorizer.paper(AnsiColor.YELLOW).print('N');
                                            // CLIPPED base not part of consensus
                                            consensus.add(' ');
                                            ++x;
                                        }
                                        ++readRef;
                                    }
                                    break;
                                }
                            case S:
                                {
                                    if (showClip) {
                                        if (testInInterval.test(readRef)) {
                                            final char readBase = baseAt.apply(readpos);
                                            if (print_this_line)
                                                colorizer.paper(AnsiColor.YELLOW).print(readBase);
                                            // CLIPPED base not part of consensus
                                            consensus.add(' ');
                                            ++x;
                                        }
                                        ++readpos;
                                        ++readRef;
                                    } else {
                                        readpos++;
                                    }
                                    break;
                                }
                            case D:
                            case N:
                                {
                                    if (testInInterval.test(readRef)) {
                                        if (print_this_line)
                                            colorizer.paper(AnsiColor.RED).print('-');
                                        // deletion not not part of consensus
                                        consensus.add(' ');
                                        ++x;
                                    }
                                    ++readRef;
                                    break;
                                }
                            case EQ:
                            case M:
                            case X:
                                {
                                    if (testInInterval.test(readRef)) {
                                        final char refBase = Character.toUpperCase(refPosToBase.apply(readRef - 1));
                                        char readBase = Character.toUpperCase(baseAt.apply(readpos));
                                        consensus.add(readBase);
                                        colorizer.pen(base2ansiColor.apply(readBase));
                                        if (op.equals(CigarOperator.X) || (refBase != 'N' && readBase != 'N' && readBase != refBase)) {
                                            colorizer.pen(AnsiColor.RED);
                                        } else if (hideBases) {
                                            if (rec.getReadNegativeStrandFlag()) {
                                                readBase = ',';
                                            } else {
                                                readBase = '.';
                                            }
                                        }
                                        if (showReadName) {
                                            final String readName = rec.getReadName();
                                            if (readpos < 0 || readpos >= readName.length()) {
                                                readBase = '_';
                                            } else {
                                                readBase = readName.charAt(readpos);
                                            }
                                        }
                                        if (rec.getReadNegativeStrandFlag()) {
                                            readBase = Character.toLowerCase(readBase);
                                        } else {
                                            readBase = Character.toUpperCase(readBase);
                                        }
                                        if (print_this_line)
                                            colorizer.print(readBase);
                                        ++x;
                                    }
                                    ++readpos;
                                    ++readRef;
                                    break;
                                }
                        }
                        ++cigarIdx;
                    }
                }
                // end of loop cigar
                ref = readRef;
            }
            // out.println( " "+ref+" "+row.get(0).getAlignmentStart()+" "+row.get(0).getCigarString()+" "+row.get(0).getReadString());
            while (x < pixelWidth) {
                if (print_this_line)
                    out.print(" ");
                ++x;
            }
            if (print_this_line)
                out.println();
            consensus.eol();
            if (out.checkError())
                break;
        }
        if (out.checkError())
            break;
        if (!this.hideConsensus && consensus.bases.stream().anyMatch(C -> C.getCoverage() > 0)) {
            out.print(margin(groupName + " CONSENSUS"));
            x = 0;
            ref = interval.getStart();
            while (x < consensus.bases.size()) {
                final char refBase = Character.toUpperCase(refPosToBase.apply(ref - 1));
                final char consensusBase = consensus.bases.get(x).getConsensus();
                if (Character.isWhitespace(consensusBase)) {
                // nothing
                } else if (refBase != 'N' && consensusBase != refBase) {
                    colorizer.pen(AnsiColor.RED);
                } else {
                    colorizer.pen(base2ansiColor.apply(consensusBase));
                }
                if (!insertIsPresentAtX.test(x))
                    ++ref;
                colorizer.print(consensusBase);
                ++x;
            }
            out.println();
        }
        if (this.numCoverageRows > 0) {
            int minCov = consensus.bases.stream().mapToInt(C -> C.getCoverage()).min().orElse(0);
            final int maxCov = consensus.bases.stream().mapToInt(C -> C.getCoverage()).max().orElse(0);
            for (int y = 0; maxCov > 0 && y < this.numCoverageRows; ++y) {
                if (minCov == maxCov)
                    minCov--;
                double fract = (maxCov - minCov) / ((double) this.numCoverageRows);
                int inverse_y = (this.numCoverageRows - 1) - y;
                int d0 = (int) ((fract) * inverse_y);
                // int d1 = (int)((fract) * (inverse_y+1));
                out.print(margin(y == 0 ? groupName + " " + maxCov : (y + 1 == this.numCoverageRows ? String.valueOf(minCov) : "")));
                for (x = 0; x < consensus.bases.size(); ++x) {
                    int depth = consensus.bases.get(x).getCoverage() - minCov;
                    colorizer.print(depth >= d0 ? BLACK_SQUARE : ' ');
                }
                out.println();
            }
        }
    }
    if (this.tabixKnownGene != null && this.indexedFastaSequenceFile != null) {
        final List<KnownGene> genes = this.tabixKnownGene.getItemsInInterval(this.interval);
        if (!genes.isEmpty()) {
            out.println(this.knownGeneUri);
            for (final KnownGene gene : genes) {
                final KnownGene.CodingRNA codingRna = gene.getCodingRNA(contigSequence);
                final KnownGene.Peptide peptide = codingRna.getPeptide();
                out.print(margin(gene.getName()));
                x = 0;
                int ref0 = this.interval.getStart() - 1;
                while (x < pixelWidth) {
                    if (insertIsPresentAtX.test(x)) {
                        out.print("*");
                        ++x;
                    } else {
                        char pepChar = ' ';
                        if (ref0 >= gene.getTxStart() && ref0 < gene.getTxEnd()) {
                            pepChar = (gene.isPositiveStrand() ? '>' : '<');
                            int pepIdx = peptide.convertGenomicToPeptideCoordinate(ref0);
                            if (pepIdx != -1) {
                                final String aa3 = GeneticCode.aminoAcidTo3Letters(peptide.charAt(pepIdx));
                                final int[] offset = peptide.convertToGenomicCoordinates(pepIdx);
                                if (offset != null && offset.length == 3 && aa3 != null && aa3.length() == 3) {
                                    if (offset[0] == ref0)
                                        pepChar = aa3.charAt(0);
                                    else if (offset[1] == ref0)
                                        pepChar = aa3.charAt(1);
                                    else if (offset[2] == ref0)
                                        pepChar = aa3.charAt(2);
                                    else
                                        pepChar = '?';
                                } else {
                                    pepChar = '?';
                                }
                            }
                        }
                        out.print(pepChar);
                        ++ref0;
                        ++x;
                    }
                }
                while (x < pixelWidth) {
                    out.print(" ");
                    ++x;
                }
                out.println();
            }
        }
        out.println();
    }
    /**
     * variant section
     */
    if (!this.vcfReaders.isEmpty() && !out.checkError()) {
        final Function<GenotypeType, Character> gTypeToSymbol = new Function<GenotypeType, Character>() {

            @Override
            public Character apply(final GenotypeType gt) {
                switch(gt) {
                    case NO_CALL:
                        return '?';
                    case HOM_REF:
                        return '0';
                    case HET:
                        return '1';
                    case HOM_VAR:
                        return '2';
                    case MIXED:
                        return 'm';
                    case UNAVAILABLE:
                        return 'u';
                    default:
                        return '.';
                }
            }
        };
        out.println();
        for (final VcfSource r : this.vcfReaders) {
            if (out.checkError())
                break;
            final VCFHeader header = r.vcfFileReader.getFileHeader();
            final CloseableIterator<VariantContext> iter = r.vcfFileReader.query(this.interval.getContig(), interval.getStart(), interval.getEnd());
            final List<VariantContext> variants = new ArrayList<>();
            while (iter.hasNext()) {
                variants.add(iter.next());
            }
            iter.close();
            if (variants.isEmpty())
                continue;
            out.println(r.vcfFile.getPath());
            if (header.hasGenotypingData()) {
                for (final String sample : header.getSampleNamesInOrder()) {
                    if (!variants.stream().map(V -> V.getGenotype(sample)).filter(G -> !hideNoCall || (hideNoCall && !G.isNoCall())).filter(G -> !hideHomRef || (hideHomRef && !G.isHomRef())).findAny().isPresent()) {
                        continue;
                    }
                    out.print(margin(sample));
                    ref = this.interval.getStart();
                    x = 0;
                    while (x < pixelWidth) {
                        if (insertIsPresentAtX.test(x)) {
                            out.print("*");
                            ++x;
                        } else {
                            char refBase = ' ';
                            for (final VariantContext ctx : variants) {
                                if (ctx.getStart() == ref) {
                                    final Genotype g = ctx.getGenotype(sample);
                                    if (g.isNoCall() && this.hideNoCall)
                                        continue;
                                    if (g.isHomRef() && this.hideHomRef)
                                        continue;
                                    refBase = gTypeToSymbol.apply(g.getType());
                                    break;
                                }
                            }
                            out.print(refBase);
                            ++ref;
                            ++x;
                        }
                    }
                    out.println();
                }
            } else // no genotype
            {
                for (final VariantContext ctx : variants) {
                    out.print(margin(String.valueOf(ctx.getStart()) + ":" + ctx.getReference().getDisplayString() + "/" + ctx.getAlternateAlleles().stream().map(A -> A.getDisplayString()).collect(Collectors.joining(","))));
                    ref = this.interval.getStart();
                    x = 0;
                    while (x < pixelWidth) {
                        if (insertIsPresentAtX.test(x)) {
                            out.print("*");
                            ++x;
                        } else {
                            out.print(ctx.getStart() == ref ? '+' : ' ');
                            ++ref;
                            ++x;
                        }
                    }
                    out.println();
                }
            }
        }
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) CloseableIterator(htsjdk.samtools.util.CloseableIterator) Parameter(com.beust.jcommander.Parameter) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VCFHeader(htsjdk.variant.vcf.VCFHeader) CigarElement(htsjdk.samtools.CigarElement) CigarOperator(htsjdk.samtools.CigarOperator) SAMRecordPartition(com.github.lindenb.jvarkit.util.samtools.SAMRecordPartition) GenomicSequence(com.github.lindenb.jvarkit.util.picard.GenomicSequence) Function(java.util.function.Function) ValidationStringency(htsjdk.samtools.ValidationStringency) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) Interval(htsjdk.samtools.util.Interval) Map(java.util.Map) XMLStreamException(javax.xml.stream.XMLStreamException) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) XMLStreamWriter(javax.xml.stream.XMLStreamWriter) SAMSequenceDictionaryExtractor(htsjdk.variant.utils.SAMSequenceDictionaryExtractor) GeneticCode(com.github.lindenb.jvarkit.util.bio.GeneticCode) CloserUtil(htsjdk.samtools.util.CloserUtil) PrintStream(java.io.PrintStream) Counter(com.github.lindenb.jvarkit.util.Counter) Predicate(java.util.function.Predicate) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) Logger(com.github.lindenb.jvarkit.util.log.Logger) GenotypeType(htsjdk.variant.variantcontext.GenotypeType) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Set(java.util.Set) SamFilterParser(com.github.lindenb.jvarkit.util.bio.samfilter.SamFilterParser) IOException(java.io.IOException) SamReader(htsjdk.samtools.SamReader) Collectors(java.util.stream.Collectors) KnownGene(com.github.lindenb.jvarkit.util.ucsc.KnownGene) TabixKnownGeneFileReader(com.github.lindenb.jvarkit.util.ucsc.TabixKnownGeneFileReader) File(java.io.File) SAMRecord(htsjdk.samtools.SAMRecord) XMLOutputFactory(javax.xml.stream.XMLOutputFactory) SamRecordFilter(htsjdk.samtools.filter.SamRecordFilter) List(java.util.List) SamInputResource(htsjdk.samtools.SamInputResource) TreeMap(java.util.TreeMap) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) Closeable(java.io.Closeable) VariantContext(htsjdk.variant.variantcontext.VariantContext) Collections(java.util.Collections) SamReaderFactory(htsjdk.samtools.SamReaderFactory) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) ArrayList(java.util.ArrayList) VariantContext(htsjdk.variant.variantcontext.VariantContext) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Predicate(java.util.function.Predicate) SamReader(htsjdk.samtools.SamReader) Function(java.util.function.Function) ArrayList(java.util.ArrayList) List(java.util.List) VCFHeader(htsjdk.variant.vcf.VCFHeader) GenomicSequence(com.github.lindenb.jvarkit.util.picard.GenomicSequence) Genotype(htsjdk.variant.variantcontext.Genotype) CigarOperator(htsjdk.samtools.CigarOperator) TreeMap(java.util.TreeMap) CigarElement(htsjdk.samtools.CigarElement) SAMRecord(htsjdk.samtools.SAMRecord) GenotypeType(htsjdk.variant.variantcontext.GenotypeType) KnownGene(com.github.lindenb.jvarkit.util.ucsc.KnownGene)

Example 42 with CigarOperator

use of htsjdk.samtools.CigarOperator in project jvarkit by lindenb.

the class VCFCombineTwoSnvs method doVcfToVcf.

@Override
protected int doVcfToVcf(final String inputName, File saveAs) {
    BufferedReader bufferedReader = null;
    htsjdk.variant.variantcontext.writer.VariantContextWriter w = null;
    SortingCollection<CombinedMutation> mutations = null;
    CloseableIterator<Variant> varIter = null;
    CloseableIterator<CombinedMutation> mutIter = null;
    Map<String, SamReader> sample2samReader = new HashMap<>();
    try {
        bufferedReader = inputName == null ? IOUtils.openStreamForBufferedReader(stdin()) : IOUtils.openURIForBufferedReading(inputName);
        final VCFUtils.CodecAndHeader cah = VCFUtils.parseHeader(bufferedReader);
        /* get VCF header */
        final VCFHeader header = cah.header;
        final Set<String> sampleNamesInOrder = new HashSet<>(header.getSampleNamesInOrder());
        LOG.info("opening REF:" + referenceFile);
        this.indexedFastaSequenceFile = new IndexedFastaSequenceFile(this.referenceFile);
        final SAMSequenceDictionary dict = this.indexedFastaSequenceFile.getSequenceDictionary();
        if (dict == null)
            throw new IOException("dictionary missing");
        if (this.bamIn != null) {
            /**
             * unroll and open bam file
             */
            for (final File bamFile : IOUtils.unrollFileCollection(Collections.singletonList(this.bamIn))) {
                LOG.info("opening BAM :" + this.bamIn);
                final SamReader samReader = SamReaderFactory.makeDefault().referenceSequence(this.referenceFile).validationStringency(ValidationStringency.LENIENT).open(this.bamIn);
                if (!samReader.hasIndex()) {
                    throw new IOException("Sam file is NOT indexed: " + bamFile);
                }
                final SAMFileHeader samHeader = samReader.getFileHeader();
                if (samHeader.getSequenceDictionary() == null || !SequenceUtil.areSequenceDictionariesEqual(dict, samReader.getFileHeader().getSequenceDictionary())) {
                    throw new IOException(bamFile + " and REF don't have the same Sequence Dictionary.");
                }
                /* get sample name */
                String sampleName = null;
                for (final SAMReadGroupRecord rg : samHeader.getReadGroups()) {
                    if (rg.getSample() == null)
                        continue;
                    if (sampleName != null && !sampleName.equals(rg.getSample())) {
                        samReader.close();
                        throw new IOException(bamFile + " Contains two samples " + sampleName + " " + rg.getSample());
                    }
                    sampleName = rg.getSample();
                }
                if (sampleName == null) {
                    samReader.close();
                    LOG.warn("no sample in " + bamFile);
                    continue;
                }
                if (!sampleNamesInOrder.contains(sampleName)) {
                    samReader.close();
                    LOG.warn("no sample " + sampleName + " in vcf");
                    continue;
                }
                sample2samReader.put(sampleName, samReader);
            }
        }
        loadKnownGenesFromUri();
        this.variants = SortingCollection.newInstance(Variant.class, new VariantCodec(), new VariantComparator(dict), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
        this.variants.setDestructiveIteration(true);
        SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header);
        String vcfLine = null;
        while ((vcfLine = bufferedReader.readLine()) != null) {
            final VariantContext ctx = progress.watch(cah.codec.decode(vcfLine));
            /* discard non SNV variant */
            if (!ctx.isVariant() || ctx.isIndel()) {
                continue;
            }
            /* find the overlapping genes : extend the interval of the variant to include the stop codon */
            final Collection<KnownGene> genes = new ArrayList<>();
            for (List<KnownGene> lkg : this.knownGenes.getOverlapping(new Interval(ctx.getContig(), Math.max(1, ctx.getStart() - 3), ctx.getEnd() + 3))) {
                genes.addAll(lkg);
            }
            final List<Allele> alternateAlleles = ctx.getAlternateAlleles();
            /* loop over overlapping genes */
            for (final KnownGene kg : genes) {
                /* loop over available alleles */
                for (int allele_idx = 0; allele_idx < alternateAlleles.size(); ++allele_idx) {
                    final Allele alt = alternateAlleles.get(allele_idx);
                    challenge(ctx, alt, kg, vcfLine);
                }
            }
        }
        progress.finish();
        this.variants.doneAdding();
        mutations = SortingCollection.newInstance(CombinedMutation.class, new MutationCodec(), new MutationComparator(dict), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
        mutations.setDestructiveIteration(true);
        final VCFFilterHeaderLine vcfFilterHeaderLine = new VCFFilterHeaderLine("TwoHaplotypes", "(number of reads carrying both mutation) < (reads carrying variant 1 + reads carrying variant 2) ");
        varIter = this.variants.iterator();
        progress = new SAMSequenceDictionaryProgress(header);
        final ArrayList<Variant> buffer = new ArrayList<>();
        for (; ; ) {
            Variant variant = null;
            if (varIter.hasNext()) {
                variant = varIter.next();
                progress.watch(variant.contig, variant.genomicPosition1);
            }
            if (variant == null || !(!buffer.isEmpty() && buffer.get(0).contig.equals(variant.contig) && buffer.get(0).transcriptName.equals(variant.transcriptName))) {
                if (!buffer.isEmpty()) {
                    for (int i = 0; i < buffer.size(); ++i) {
                        final Variant v1 = buffer.get(i);
                        for (int j = i + 1; j < buffer.size(); ++j) {
                            final Variant v2 = buffer.get(j);
                            if (v1.codonStart() != v2.codonStart())
                                continue;
                            if (v1.positionInCodon() == v2.positionInCodon())
                                continue;
                            if (!v1.wildCodon.equals(v2.wildCodon)) {
                                throw new IllegalStateException();
                            }
                            final StringBuilder combinedCodon = new StringBuilder(v1.wildCodon);
                            combinedCodon.setCharAt(v1.positionInCodon(), v1.mutCodon.charAt(v1.positionInCodon()));
                            combinedCodon.setCharAt(v2.positionInCodon(), v2.mutCodon.charAt(v2.positionInCodon()));
                            final String pwild = new ProteinCharSequence(v1.wildCodon).getString();
                            final String p1 = new ProteinCharSequence(v1.mutCodon).getString();
                            final String p2 = new ProteinCharSequence(v2.mutCodon).getString();
                            final String pCombined = new ProteinCharSequence(combinedCodon).getString();
                            final String combinedSO;
                            final String combinedType;
                            /* both AA are synonymous, while combined is not */
                            if (!pCombined.equals(pwild) && p1.equals(pwild) && p2.equals(pwild)) {
                                combinedType = "combined_is_nonsynonymous";
                                if (pCombined.equals("*")) {
                                    /* http://www.sequenceontology.org/browser/current_svn/term/SO:0001587 */
                                    combinedSO = "stop_gained";
                                } else if (pwild.equals("*")) {
                                    /* http://www.sequenceontology.org/browser/current_svn/term/SO:0002012 */
                                    combinedSO = "stop_lost";
                                } else {
                                    /* http://www.sequenceontology.org/miso/current_svn/term/SO:0001992 */
                                    combinedSO = "nonsynonymous_variant";
                                }
                            } else if (!pCombined.equals(p1) && !pCombined.equals(p2) && !pCombined.equals(pwild)) {
                                combinedType = "combined_is_new";
                                if (pCombined.equals("*")) {
                                    /* http://www.sequenceontology.org/browser/current_svn/term/SO:0001587 */
                                    combinedSO = "stop_gained";
                                } else {
                                    /* http://www.sequenceontology.org/miso/current_svn/term/SO:0001992 */
                                    combinedSO = "nonsynonymous_variant";
                                }
                            } else {
                                combinedType = null;
                                combinedSO = null;
                            }
                            /**
                             * ok, there is something interesting here ,
                             * create two new Mutations carrying the
                             * two variants
                             */
                            if (combinedSO != null) {
                                /**
                                 * grantham score is max found combined vs (p1/p2/wild)
                                 */
                                int grantham_score = GranthamScore.score(pCombined.charAt(0), pwild.charAt(0));
                                grantham_score = Math.max(grantham_score, GranthamScore.score(pCombined.charAt(0), p1.charAt(0)));
                                grantham_score = Math.max(grantham_score, GranthamScore.score(pCombined.charAt(0), p2.charAt(0)));
                                /**
                                 * info that will be displayed in the vcf
                                 */
                                final Map<String, Object> info1 = v1.getInfo(v2);
                                final Map<String, Object> info2 = v2.getInfo(v1);
                                // filter for this combined: default it fails the filter
                                String filter = vcfFilterHeaderLine.getID();
                                final Map<String, Object> combinedMap = new LinkedHashMap<>();
                                combinedMap.put("CombinedCodon", combinedCodon);
                                combinedMap.put("CombinedAA", pCombined);
                                combinedMap.put("CombinedSO", combinedSO);
                                combinedMap.put("CombinedType", combinedType);
                                combinedMap.put("GranthamScore", grantham_score);
                                info1.putAll(combinedMap);
                                info2.putAll(combinedMap);
                                final Map<String, CoverageInfo> sample2coverageInfo = new HashMap<>(sample2samReader.size());
                                final int chromStart = Math.min(v1.genomicPosition1, v2.genomicPosition1);
                                final int chromEnd = Math.max(v1.genomicPosition1, v2.genomicPosition1);
                                /* get phasing info for each sample*/
                                for (final String sampleName : sample2samReader.keySet()) {
                                    final SamReader samReader = sample2samReader.get(sampleName);
                                    final CoverageInfo covInfo = new CoverageInfo();
                                    sample2coverageInfo.put(sampleName, covInfo);
                                    SAMRecordIterator iter = null;
                                    try {
                                        iter = samReader.query(v1.contig, chromStart, chromEnd, false);
                                        while (iter.hasNext()) {
                                            final SAMRecord rec = iter.next();
                                            if (rec.getReadUnmappedFlag())
                                                continue;
                                            if (rec.isSecondaryOrSupplementary())
                                                continue;
                                            if (rec.getDuplicateReadFlag())
                                                continue;
                                            if (rec.getReadFailsVendorQualityCheckFlag())
                                                continue;
                                            // get DEPTh for variant 1
                                            if (rec.getAlignmentStart() <= v1.genomicPosition1 && v1.genomicPosition1 <= rec.getAlignmentEnd()) {
                                                covInfo.depth1++;
                                            }
                                            // get DEPTh for variant 2
                                            if (rec.getAlignmentStart() <= v2.genomicPosition1 && v2.genomicPosition1 <= rec.getAlignmentEnd()) {
                                                covInfo.depth2++;
                                            }
                                            if (rec.getAlignmentEnd() < chromEnd)
                                                continue;
                                            if (rec.getAlignmentStart() > chromStart)
                                                continue;
                                            final Cigar cigar = rec.getCigar();
                                            if (cigar == null)
                                                continue;
                                            final byte[] bases = rec.getReadBases();
                                            if (bases == null)
                                                continue;
                                            int refpos1 = rec.getAlignmentStart();
                                            int readpos = 0;
                                            boolean found_variant1_on_this_read = false;
                                            boolean found_variant2_on_this_read = false;
                                            /**
                                             * loop over cigar
                                             */
                                            for (final CigarElement ce : cigar.getCigarElements()) {
                                                final CigarOperator op = ce.getOperator();
                                                switch(op) {
                                                    case P:
                                                        continue;
                                                    case S:
                                                    case I:
                                                        readpos += ce.getLength();
                                                        break;
                                                    case D:
                                                    case N:
                                                        refpos1 += ce.getLength();
                                                        break;
                                                    case H:
                                                        continue;
                                                    case EQ:
                                                    case M:
                                                    case X:
                                                        for (int x = 0; x < ce.getLength(); ++x) {
                                                            if (refpos1 == v1.genomicPosition1 && same(bases[readpos], v1.altAllele)) {
                                                                found_variant1_on_this_read = true;
                                                            } else if (refpos1 == v2.genomicPosition1 && same(bases[readpos], v2.altAllele)) {
                                                                found_variant2_on_this_read = true;
                                                            }
                                                            refpos1++;
                                                            readpos++;
                                                        }
                                                        break;
                                                    default:
                                                        throw new IllegalStateException(op.name());
                                                }
                                                /* skip remaining bases after last variant */
                                                if (refpos1 > chromEnd)
                                                    break;
                                            }
                                            /* sum-up what we found */
                                            if (found_variant1_on_this_read && found_variant2_on_this_read) {
                                                covInfo.count_reads_having_both_variants++;
                                            } else if (!found_variant1_on_this_read && !found_variant2_on_this_read) {
                                                covInfo.count_reads_having_no_variants++;
                                            } else if (found_variant1_on_this_read) {
                                                covInfo.count_reads_having_variant1++;
                                            } else if (found_variant2_on_this_read) {
                                                covInfo.count_reads_having_variant2++;
                                            }
                                        }
                                    /* end of loop over reads */
                                    } finally {
                                        iter.close();
                                        iter = null;
                                    }
                                    info1.put("N_READS_BOTH_VARIANTS_" + sampleName, covInfo.count_reads_having_both_variants);
                                    info2.put("N_READS_BOTH_VARIANTS_" + sampleName, covInfo.count_reads_having_both_variants);
                                    info1.put("N_READS_NO_VARIANTS_" + sampleName, covInfo.count_reads_having_no_variants);
                                    info2.put("N_READS_NO_VARIANTS_" + sampleName, covInfo.count_reads_having_no_variants);
                                    info1.put("N_READS_TOTAL_" + sampleName, covInfo.count_reads_having_both_variants + covInfo.count_reads_having_no_variants + covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2);
                                    info2.put("N_READS_TOTAL_" + sampleName, covInfo.count_reads_having_both_variants + covInfo.count_reads_having_no_variants + covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2);
                                    // count for variant 1
                                    info1.put("N_READS_ONLY_1_" + sampleName, covInfo.count_reads_having_variant1);
                                    info1.put("N_READS_ONLY_2_" + sampleName, covInfo.count_reads_having_variant2);
                                    info1.put("DEPTH_1_" + sampleName, covInfo.depth1);
                                    // inverse previous count
                                    info2.put("N_READS_ONLY_1_" + sampleName, covInfo.count_reads_having_variant2);
                                    info2.put("N_READS_ONLY_2_" + sampleName, covInfo.count_reads_having_variant1);
                                    info2.put("DEPTH_2_" + sampleName, covInfo.depth2);
                                    /* number of reads with both variant is greater than
									 * reads carrying only one variant: reset the filter 
									 */
                                    if (2 * covInfo.count_reads_having_both_variants > (covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2)) {
                                        /* reset filter */
                                        filter = VCFConstants.UNFILTERED;
                                        info1.put("FILTER_1_" + sampleName, ".");
                                        info2.put("FILTER_2_" + sampleName, ".");
                                    } else {
                                        info1.put("FILTER_1_" + sampleName, vcfFilterHeaderLine.getID());
                                        info2.put("FILTER_2_" + sampleName, vcfFilterHeaderLine.getID());
                                    }
                                }
                                /* end of loop over bams */
                                final CombinedMutation m1 = new CombinedMutation();
                                m1.contig = v1.contig;
                                m1.genomicPosition1 = v1.genomicPosition1;
                                m1.id = v1.id;
                                m1.refAllele = v1.refAllele;
                                m1.altAllele = v1.altAllele;
                                m1.vcfLine = v1.vcfLine;
                                m1.info = mapToString(info1);
                                m1.filter = filter;
                                m1.grantham_score = grantham_score;
                                m1.sorting_id = ID_GENERATOR++;
                                mutations.add(m1);
                                final CombinedMutation m2 = new CombinedMutation();
                                m2.contig = v2.contig;
                                m2.genomicPosition1 = v2.genomicPosition1;
                                m2.id = v2.id;
                                m2.refAllele = v2.refAllele;
                                m2.altAllele = v2.altAllele;
                                m2.vcfLine = v2.vcfLine;
                                m2.info = mapToString(info2);
                                m2.filter = filter;
                                m2.grantham_score = grantham_score;
                                m2.sorting_id = ID_GENERATOR++;
                                mutations.add(m2);
                            }
                        }
                    }
                }
                buffer.clear();
                if (variant == null)
                    break;
            }
            buffer.add(variant);
        }
        progress.finish();
        mutations.doneAdding();
        varIter.close();
        varIter = null;
        variants.cleanup();
        variants = null;
        final ArrayList<CombinedMutation> mBuffer = new ArrayList<>();
        final VCFHeader header2 = new VCFHeader(header);
        header2.addMetaDataLine(new VCFHeaderLine(getProgramName() + "AboutQUAL", "QUAL is filled with Grantham Score  http://www.ncbi.nlm.nih.gov/pubmed/4843792"));
        final StringBuilder infoDesc = new StringBuilder("Variant affected by two distinct mutation. Format is defined in the INFO column. ");
        final VCFInfoHeaderLine infoHeaderLine = new VCFInfoHeaderLine("CodonVariant", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, infoDesc.toString());
        super.addMetaData(header2);
        header2.addMetaDataLine(infoHeaderLine);
        if (!sample2samReader.isEmpty()) {
            header2.addMetaDataLine(vcfFilterHeaderLine);
        }
        w = super.openVariantContextWriter(saveAs);
        w.writeHeader(header2);
        progress = new SAMSequenceDictionaryProgress(header);
        mutIter = mutations.iterator();
        for (; ; ) {
            CombinedMutation mutation = null;
            if (mutIter.hasNext()) {
                mutation = mutIter.next();
                progress.watch(mutation.contig, mutation.genomicPosition1);
            }
            if (mutation == null || !(!mBuffer.isEmpty() && mBuffer.get(0).contig.equals(mutation.contig) && mBuffer.get(0).genomicPosition1 == mutation.genomicPosition1 && mBuffer.get(0).refAllele.equals(mutation.refAllele))) {
                if (!mBuffer.isEmpty()) {
                    // default grantham score used in QUAL
                    int grantham_score = -1;
                    // default filter fails
                    String filter = vcfFilterHeaderLine.getID();
                    final CombinedMutation first = mBuffer.get(0);
                    final Set<String> info = new HashSet<>();
                    final VariantContext ctx = cah.codec.decode(first.vcfLine);
                    final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
                    vcb.chr(first.contig);
                    vcb.start(first.genomicPosition1);
                    vcb.stop(first.genomicPosition1 + first.refAllele.length() - 1);
                    if (!first.id.equals(VCFConstants.EMPTY_ID_FIELD))
                        vcb.id(first.id);
                    for (final CombinedMutation m : mBuffer) {
                        info.add(m.info);
                        grantham_score = Math.max(grantham_score, m.grantham_score);
                        if (VCFConstants.UNFILTERED.equals(m.filter)) {
                            // at least one SNP is ok one this line
                            filter = null;
                        }
                    }
                    vcb.unfiltered();
                    if (filter != null && !sample2samReader.isEmpty()) {
                        vcb.filter(filter);
                    } else {
                        vcb.passFilters();
                    }
                    vcb.attribute(infoHeaderLine.getID(), new ArrayList<String>(info));
                    if (grantham_score > 0) {
                        vcb.log10PError(grantham_score / -10.0);
                    } else {
                        vcb.log10PError(VariantContext.NO_LOG10_PERROR);
                    }
                    w.add(vcb.make());
                }
                mBuffer.clear();
                if (mutation == null)
                    break;
            }
            mBuffer.add(mutation);
        }
        progress.finish();
        mutIter.close();
        mutations.cleanup();
        mutations = null;
        return RETURN_OK;
    } catch (Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(this.indexedFastaSequenceFile);
        CloserUtil.close(mutIter);
        CloserUtil.close(varIter);
        if (this.variants != null)
            this.variants.cleanup();
        if (mutations != null)
            mutations.cleanup();
        this.variants = null;
        for (SamReader r : sample2samReader.values()) CloserUtil.close(r);
        CloserUtil.close(w);
        CloserUtil.close(bufferedReader);
    }
}
Also used : VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) HashMap(java.util.HashMap) LinkedHashMap(java.util.LinkedHashMap) ArrayList(java.util.ArrayList) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) LinkedHashMap(java.util.LinkedHashMap) HashSet(java.util.HashSet) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) SAMRecord(htsjdk.samtools.SAMRecord) KnownGene(com.github.lindenb.jvarkit.util.ucsc.KnownGene) SAMFileHeader(htsjdk.samtools.SAMFileHeader) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) File(java.io.File) Interval(htsjdk.samtools.util.Interval) VCFUtils(com.github.lindenb.jvarkit.util.vcf.VCFUtils) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) VariantContext(htsjdk.variant.variantcontext.VariantContext) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) SamReader(htsjdk.samtools.SamReader) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) VCFHeader(htsjdk.variant.vcf.VCFHeader) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) IOException(java.io.IOException) VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine) IOException(java.io.IOException) Allele(htsjdk.variant.variantcontext.Allele) Cigar(htsjdk.samtools.Cigar) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BufferedReader(java.io.BufferedReader)

Example 43 with CigarOperator

use of htsjdk.samtools.CigarOperator in project jvarkit by lindenb.

the class CigarOpPerPositionChartFactory method build.

@Override
public StackedBarChart<String, Number> build() {
    final CategoryAxis xAxis = new CategoryAxis();
    xAxis.setLabel("Position in Read");
    final NumberAxis yAxis = new NumberAxis();
    yAxis.setLabel("Count");
    final List<XYChart.Series<String, Number>> base2count = new ArrayList<>();
    for (final CigarOperator cigarop : CigarOperator.values()) {
        if (cigarop == CigarOperator.P)
            continue;
        final XYChart.Series<String, Number> serie = new XYChart.Series<String, Number>();
        serie.setName(cigarop.name());
        base2count.add(serie);
        for (int i = 0; i < this.cigar2pos2count.size(); ++i) {
            serie.getData().add(new XYChart.Data<String, Number>(String.valueOf(i + 1), this.cigar2pos2count.get(i).count(cigarop)));
        }
    }
    final StackedBarChart<String, Number> sbc = new StackedBarChart<String, Number>(xAxis, yAxis);
    sbc.setTitle(getName());
    sbc.getData().addAll(base2count);
    sbc.setCategoryGap(0.2);
    return sbc;
}
Also used : NumberAxis(javafx.scene.chart.NumberAxis) ArrayList(java.util.ArrayList) StackedBarChart(javafx.scene.chart.StackedBarChart) CigarOperator(htsjdk.samtools.CigarOperator) CategoryAxis(javafx.scene.chart.CategoryAxis) XYChart(javafx.scene.chart.XYChart)

Example 44 with CigarOperator

use of htsjdk.samtools.CigarOperator 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);
    }
}
Also used : PrintStream(java.io.PrintStream) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) GenomicSequence(com.github.lindenb.jvarkit.util.picard.GenomicSequence) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) SamReader(htsjdk.samtools.SamReader) Cigar(htsjdk.samtools.Cigar) SAMRecord(htsjdk.samtools.SAMRecord)

Example 45 with CigarOperator

use of htsjdk.samtools.CigarOperator in project jvarkit by lindenb.

the class Biostar214299 method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.positionFile == null) {
        LOG.error("position File is not defined.");
        return -1;
    }
    final String UNAFFECTED_SAMPLE = "UNAFFECTED";
    final String AMBIGOUS_SAMPLE = "AMBIGOUS";
    final String UNMAPPED = "UNMAPPED";
    SamReader sfr = null;
    SAMFileWriter sfw = null;
    final IntervalTreeMap<Position> positionsTreeMap = new IntervalTreeMap<>();
    final Set<String> samples = new HashSet<>();
    try {
        sfr = openSamReader(oneFileOrNull(args));
        final SAMFileHeader header = sfr.getFileHeader();
        final SAMSequenceDictionary dict = header.getSequenceDictionary();
        if (dict == null) {
            LOG.error("Dictionary missing in input sam");
            return -1;
        }
        try (BufferedReader br = IOUtils.openFileForBufferedReading(this.positionFile)) {
            String line;
            while ((line = br.readLine()) != null) {
                if (line.trim().isEmpty() || line.startsWith("#"))
                    continue;
                final String[] tokens = line.split("[\t]");
                if (tokens.length < 4) {
                    LOG.error("Not enough columns in " + line);
                    return -1;
                }
                final String contig = tokens[0];
                if (dict.getSequence(contig) == null) {
                    LOG.error("No such contig in input's sam dictionary: " + contig);
                    return -1;
                }
                final int refpos = Integer.parseInt(tokens[1]);
                final Interval interval = new Interval(contig, refpos, refpos);
                Position position = positionsTreeMap.get(interval);
                if (position == null) {
                    position = new Position();
                    // position.contig = contig;
                    position.refpos = refpos;
                    positionsTreeMap.put(interval, position);
                }
                final String bases = tokens[2].toUpperCase();
                if (bases.length() != 1 || !bases.matches("[ATGC]")) {
                    LOG.error("in " + line + " bases should be one letter an ATGC");
                    return -1;
                }
                if (position.base2sample.containsKey(bases.charAt(0))) {
                    LOG.error("in " + line + " bases already defined for this position");
                    return -1;
                }
                final String sampleName = tokens[3].trim();
                if (sampleName.isEmpty()) {
                    LOG.error("sample name cannot be empty");
                    return -1;
                }
                samples.add(sampleName);
                position.base2sample.put(bases.charAt(0), sampleName);
            }
        } catch (final IOException err) {
            LOG.error(err);
            return -1;
        }
        if (samples.contains(UNAFFECTED_SAMPLE)) {
            LOG.error("Sample cannot be named " + UNAFFECTED_SAMPLE);
            return -1;
        }
        if (samples.contains(AMBIGOUS_SAMPLE)) {
            LOG.error("Sample cannot be named " + AMBIGOUS_SAMPLE);
            return -1;
        }
        if (samples.contains(UNMAPPED)) {
            LOG.error("Sample cannot be named " + UNMAPPED);
            return -1;
        }
        samples.add(UNAFFECTED_SAMPLE);
        samples.add(AMBIGOUS_SAMPLE);
        samples.add(UNMAPPED);
        final SAMFileHeader newHeader = new SAMFileHeader();
        newHeader.setSortOrder(header.getSortOrder());
        newHeader.setSequenceDictionary(dict);
        newHeader.addComment("generated with " + getProgramName() + " " + getVersion() + " Pierre Lindenbaum : " + getProgramCommandLine());
        /* create groups */
        for (final String sample : samples) {
            final SAMReadGroupRecord rg = new SAMReadGroupRecord(sample);
            rg.setSample(sample);
            rg.setLibrary(sample);
            newHeader.addReadGroup(rg);
        }
        sfw = this.writingBamArgs.openSAMFileWriter(this.outputFile, newHeader, true);
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header).logger(LOG);
        final SAMRecordIterator iter = sfr.iterator();
        while (iter.hasNext()) {
            final SAMRecord rec = progress.watch(iter.next());
            rec.setAttribute("RG", null);
            if (rec.getReadUnmappedFlag()) {
                rec.setAttribute("RG", UNMAPPED);
                sfw.addAlignment(rec);
                continue;
            }
            final Cigar cigar = rec.getCigar();
            final Collection<Position> snps = positionsTreeMap.getContained(new Interval(rec.getContig(), rec.getUnclippedStart(), rec.getUnclippedEnd()));
            if (snps == null || snps.isEmpty()) {
                rec.setAttribute("RG", UNAFFECTED_SAMPLE);
                sfw.addAlignment(rec);
                continue;
            }
            final Map<Integer, Position> index2pos = snps.stream().collect(Collectors.toMap(P -> P.refpos, P -> P));
            final Set<String> selectedSamples = new HashSet<>();
            final byte[] bases = rec.getReadBases();
            if (bases == null || bases.equals(SAMRecord.NULL_SEQUENCE)) {
                LOG.error("Bases missing in read " + rec);
                return -1;
            }
            int refPos1 = rec.getUnclippedStart();
            int readPos0 = 0;
            for (final CigarElement ce : cigar.getCigarElements()) {
                final CigarOperator op = ce.getOperator();
                final boolean consummeReadBaseOrSoftClip = op.consumesReadBases() || op.equals(CigarOperator.S);
                if (op.consumesReferenceBases() && consummeReadBaseOrSoftClip) {
                    for (int i = 0; i < ce.getLength(); ++i) {
                        final int nowRefPos1 = (refPos1 + i);
                        final int nowReadPos0 = (readPos0 + i);
                        final Position position = index2pos.get(nowRefPos1);
                        if (position == null)
                            continue;
                        if (nowReadPos0 >= bases.length)
                            continue;
                        final char base = (char) Character.toUpperCase(bases[nowReadPos0]);
                        final String sample = position.base2sample.get(base);
                        if (sample == null)
                            continue;
                        selectedSamples.add(sample);
                        index2pos.remove(nowRefPos1);
                        if (index2pos.isEmpty())
                            break;
                    }
                }
                if (op.consumesReferenceBases())
                    refPos1 += ce.getLength();
                if (consummeReadBaseOrSoftClip || op.equals(CigarOperator.H)) {
                    readPos0 += ce.getLength();
                }
            }
            if (selectedSamples.isEmpty()) {
                rec.setAttribute("RG", UNAFFECTED_SAMPLE);
            } else if (selectedSamples.size() == 1) {
                rec.setAttribute("RG", selectedSamples.iterator().next());
            } else {
                rec.setAttribute("RG", AMBIGOUS_SAMPLE);
            }
            sfw.addAlignment(rec);
        }
        progress.finish();
        LOG.info("done");
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(sfr);
        CloserUtil.close(sfw);
    }
}
Also used : Cigar(htsjdk.samtools.Cigar) Program(com.github.lindenb.jvarkit.util.jcommander.Program) Parameter(com.beust.jcommander.Parameter) CigarElement(htsjdk.samtools.CigarElement) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) CigarOperator(htsjdk.samtools.CigarOperator) HashMap(java.util.HashMap) Term(com.github.lindenb.semontology.Term) SAMFileHeader(htsjdk.samtools.SAMFileHeader) ParametersDelegate(com.beust.jcommander.ParametersDelegate) HashSet(java.util.HashSet) Interval(htsjdk.samtools.util.Interval) Map(java.util.Map) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) CloserUtil(htsjdk.samtools.util.CloserUtil) IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) Collection(java.util.Collection) Logger(com.github.lindenb.jvarkit.util.log.Logger) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Set(java.util.Set) IOException(java.io.IOException) SAMFileWriter(htsjdk.samtools.SAMFileWriter) SamReader(htsjdk.samtools.SamReader) Collectors(java.util.stream.Collectors) File(java.io.File) SAMRecord(htsjdk.samtools.SAMRecord) List(java.util.List) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) BufferedReader(java.io.BufferedReader) SAMRecordIterator(htsjdk.samtools.SAMRecordIterator) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) SamReader(htsjdk.samtools.SamReader) HashSet(java.util.HashSet) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) SAMFileWriter(htsjdk.samtools.SAMFileWriter) IOException(java.io.IOException) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) IOException(java.io.IOException) Cigar(htsjdk.samtools.Cigar) SAMRecord(htsjdk.samtools.SAMRecord) BufferedReader(java.io.BufferedReader) SAMFileHeader(htsjdk.samtools.SAMFileHeader) IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap) Interval(htsjdk.samtools.util.Interval)

Aggregations

CigarOperator (htsjdk.samtools.CigarOperator)48 CigarElement (htsjdk.samtools.CigarElement)31 Cigar (htsjdk.samtools.Cigar)24 SAMRecord (htsjdk.samtools.SAMRecord)22 ArrayList (java.util.ArrayList)17 SAMRecordIterator (htsjdk.samtools.SAMRecordIterator)15 SamReader (htsjdk.samtools.SamReader)15 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)13 SAMFileHeader (htsjdk.samtools.SAMFileHeader)12 File (java.io.File)12 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)11 List (java.util.List)11 Parameter (com.beust.jcommander.Parameter)9 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)9 Logger (com.github.lindenb.jvarkit.util.log.Logger)9 SAMReadGroupRecord (htsjdk.samtools.SAMReadGroupRecord)9 CloserUtil (htsjdk.samtools.util.CloserUtil)9 Program (com.github.lindenb.jvarkit.util.jcommander.Program)8 IndexedFastaSequenceFile (htsjdk.samtools.reference.IndexedFastaSequenceFile)8 HashMap (java.util.HashMap)8