X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fdatamodel%2FCigarParser.java;h=a42c2b1eea4be5b467356391cf30f59b04966bd7;hb=429bbdae3d5229ec1b92ec59f4b5a25bcfcf13f8;hp=0abe9bfb2e634f1bad7360e77196c8751042919a;hpb=6f516287c0c9fa3e3f13c8e0c9b95dd99045d22a;p=jalview.git diff --git a/src/jalview/datamodel/CigarParser.java b/src/jalview/datamodel/CigarParser.java index 0abe9bf..a42c2b1 100644 --- a/src/jalview/datamodel/CigarParser.java +++ b/src/jalview/datamodel/CigarParser.java @@ -30,10 +30,13 @@ public class CigarParser * @param alignmentStart * start of alignment to be used as offset when padding reads with * gaps + * @param seq + * sequence to be annotated with inserts * @return string representing read with gaps, clipping etc applied */ public String parseCigarToSequence(SAMRecord rec, - SortedMap insertions, int alignmentStart) + SortedMap insertions, int alignmentStart, + SequenceI seq) { StringBuilder newRead = new StringBuilder(); Iterator it = rec.getCigar().getCigarElements() @@ -86,7 +89,7 @@ public class CigarParser Iterator> iit = seqInserts.entrySet() .iterator(); - String nextSegment = applyCigarOp(el, next, rec, iit, override); + String nextSegment = applyCigarOp(el, next, rec, iit, override, seq); newRead.append(nextSegment); if (el.getOperator().consumesReferenceBases()) @@ -124,11 +127,12 @@ public class CigarParser * an optional which can override a * in it to change the length. Set to null if not * used. + * @param seq * @return */ private String applyCigarOp(CigarElement el, int next, SAMRecord rec, Iterator> it, - int[] override) + int[] override, SequenceI seq) { StringBuilder newRead = new StringBuilder(); String read = rec.getReadString(); @@ -143,7 +147,7 @@ public class CigarParser case EQ: // matched or mismatched residues newRead = applyCigarConsumeReadAndRef(next, length, rec, it, - override); + override, seq); break; case N: // intron in RNA case D: // deletion @@ -161,11 +165,19 @@ public class CigarParser // the reference sequence and other reads should have been gapped for // this insertion, so just add in the residues newRead.append(read.substring(nextPos, nextPos + length)); + seq.addSequenceFeature(new SequenceFeature("INSERTION", "", + nextPos + 1, nextPos + length, "bamfile")); break; case S: // soft clipping - just skip this bit of the read + seq.addSequenceFeature(new SequenceFeature("SOFTCLIP", "", + nextPos + 1, nextPos + length, "bamfile")); + break; case H: // hard clipping - this stretch will not appear in the read + seq.addSequenceFeature(new SequenceFeature("HARDCLIP", "", + nextPos + 1, nextPos + length, "bamfile")); + break; default: break; } @@ -186,11 +198,12 @@ public class CigarParser * iterator over insertions in the CIGAR region * @param override * insertions which have already been accounted for + * @param seq * @return */ private StringBuilder applyCigarConsumeReadAndRef(int next, int length, SAMRecord rec, Iterator> it, - int[] override) + int[] override, SequenceI seq) { StringBuilder newRead = new StringBuilder(); String read = rec.getReadString(); @@ -262,52 +275,116 @@ public class CigarParser Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added */ - public SortedMap getInsertions(Iterator it) + public SortedMap[] getInsertions(Iterator it) { - SortedMap inserts = new TreeMap<>(); - while (it.hasNext()) + return getInsertions(it, null); + } + + public int firstAlColumn = -1, firstRposCol = -1; + + public SortedMap[] getInsertions(Iterator it, + Range[] extent) + { + SortedMap sinserts[] = new SortedMap[] { + new TreeMap<>(), new TreeMap<>() }; + Range xten = extent != null && extent[0] != null ? extent[0] : null; + do { - // check each record for insertions in the CIGAR string + // check extent of read SAMRecord rec = it.next(); + if (extent != null) + { + + int nstart = rec.getAlignmentStart(); + if (nstart != 0) + { + nstart = rec.getReferencePositionAtReadPosition(nstart); + } + int nend = rec.getAlignmentEnd(); + if (nend != 0) + { + nend = rec.getReferencePositionAtReadPosition(nend); + } + + xten = new Range( + nstart <= 0 ? xten.start : Math.min(nstart, xten.start), + nend == 0 ? xten.end : Math.max(nend, xten.end)); + } + + // check each record for insertions in the CIGAR string + + // pick the insert map to update + int insmap = rec.getReadNegativeStrandFlag() ? 1 : 0; + SortedMap inserts = sinserts[insmap]; + Iterator cit = rec.getCigar().getCigarElements() .iterator(); - int next = 1; + int refLocation = rec.getAlignmentStart(); + int rloc = 0; + int alcol = 0; + int alpos = 0; while (cit.hasNext()) { CigarElement el = cit.next(); switch (el.getOperator()) { - case I: - // add to insertions list, and move along the read - // location is the start of next CIGAR segment - // because getReferencePositionAtReadPosition returns 0 for read - // positions which are not in the reference (like insertions) - int refLocation = rec.getReferencePositionAtReadPosition( - next + el.getLength()); - - // if there's already an insertion at this location, keep the longest - // insertion; if there's no insertion keep this one - if (!inserts.containsKey(refLocation) - || (inserts.containsKey(refLocation) - && inserts.get(refLocation) < el.getLength())) - { - inserts.put(refLocation, el.getLength()); - } - next += el.getLength(); - break; - case M: case S: - // match to reference, or soft clip, move along the read - next += el.getLength(); + case H: + // ignore soft/hardclipped + break; default: - // deletions, introns etc don't consume any residues from the read - break; + if (el.getOperator().consumesReadBases() + && !el.getOperator().consumesReferenceBases()) + { + // add to insertions list, and move along the read + // location is the start of next CIGAR segment + // because getReferencePositionAtReadPosition returns 0 for read + // positions which are not in the reference (like insertions) + + // if there's already an insertion at this location, keep the + // longest + // insertion; if there's no insertion keep this one + if (!inserts.containsKey(refLocation) + || (inserts.containsKey(refLocation) + && inserts.get(refLocation) < el.getLength())) + { + inserts.put(refLocation, el.getLength()); + } + break; + } + if (el.getOperator().consumesReferenceBases()) + { + if (el.getOperator().consumesReadBases() && alcol == 0 + && refLocation + 1 > extent[0].start) + { + // first match position : last rloc + first match op + alcol = rloc + 1; + alpos = rec.getReferencePositionAtReadPosition(rloc + 1); + } + refLocation += el.getLength(); + } + + } + if (el.getOperator().consumesReadBases() + || el.getOperator().consumesReferenceBases()) + { + rloc += el.getLength(); } } - + // Wrong: hack that fails to relocate reference to correct place in + // sequence + if (firstAlColumn < 0) + { + firstAlColumn = alcol; + firstRposCol = alpos; + } + } while (it.hasNext()); + if (xten != null) + { + extent[0] = xten; } - return inserts; + return sinserts; } /**