src/jalview/datamodel/CigarParser.java

   1 package jalview.datamodel;
   2
   3 import java.util.Arrays;
   4 import java.util.Iterator;
   5 import java.util.Map;
   6 import java.util.Map.Entry;
   7 import java.util.SortedMap;
   8 import java.util.TreeMap;
   9
  10 import htsjdk.samtools.CigarElement;
  11 import htsjdk.samtools.CigarOperator;
  12 import htsjdk.samtools.SAMRecord;
  13
  14 public class CigarParser
  15 {
  16   private final char gapChar;
  17
  18   public CigarParser(char gap)
  19   {
  20     gapChar = gap;
  21   }
  22
  23   /**
  24    * Apply the CIGAR string to a read sequence and return the updated read.
  25    *
  26    * @param rec
  27    *          the SAM record for the read
  28    * @param insertions
  29    *          a map of inserted positions and lengths
  30    * @param alignmentStart
  31    *          start of alignment to be used as offset when padding reads with
  32    *          gaps
  33    * @param seq
  34    *          sequence to be annotated with inserts
  35    * @return string representing read with gaps, clipping etc applied
  36    */
  37   public String parseCigarToSequence(SAMRecord rec,
  38           SortedMap<Integer, Integer> insertions, int alignmentStart,
  39           SequenceI seq)
  40   {
  41     StringBuilder newRead = new StringBuilder();
  42     Iterator<CigarElement> it = rec.getCigar().getCigarElements()
  43             .iterator();
  44
  45     int next = 0;
  46     int refnext = 0;
  47
  48     // pad with spaces before read
  49     // first count gaps needed to pad to reference position
  50     // NB position is 1-based so number of gaps = pos-1
  51     int gaps = rec.getStart() - alignmentStart;
  52
  53     // now count gaps to pad for insertions in other reads
  54     int insertCount = countInsertionsBeforeRead(rec, insertions);
  55     addGaps(newRead, gaps + insertCount);
  56
  57     int lastinserts = 0;
  58     while (it.hasNext())
  59     {
  60       CigarElement el = it.next();
  61       int length = el.getLength();
  62
  63       // which insertions coincide with current CIGAR operation?
  64       SortedMap<Integer, Integer> seqInserts = insertions.subMap(
  65               rec.getStart() + refnext, rec.getStart() + refnext + length);
  66
  67       int[] override = null;
  68       if (lastinserts > 0)
  69       {
  70         if (!seqInserts.containsKey(rec.getStart() + refnext))
  71         {
  72           // ERROR
  73           int pos = rec.getStart() + refnext;
  74           System.out.println(
  75                   "Read " + rec.getReadName() + " has an insertion at "
  76                           + pos + " which is not in seqInserts");
  77         }
  78         else
  79         {
  80           // we just inserted something
  81           // remove these inserts now - should be very first entry
  82           int count = seqInserts.get(rec.getStart() + refnext);
  83           override = new int[] { rec.getStart() + refnext,
  84               count - lastinserts };
  85           lastinserts = 0;
  86         }
  87       }
  88
  89       Iterator<Map.Entry<Integer, Integer>> iit = seqInserts.entrySet()
  90               .iterator();
  91
  92       String nextSegment = applyCigarOp(el, next, rec, iit, override, seq);
  93       newRead.append(nextSegment);
  94
  95       if (el.getOperator().consumesReferenceBases())
  96       {
  97         refnext += length;
  98       }
  99       if (el.getOperator().consumesReadBases())
 100       {
 101         next += length;
 102       }
 103       if (el.getOperator().equals(CigarOperator.I))
 104       {
 105         // count how many insertions we made
 106         lastinserts = length;
 107       }
 108     }
 109     return newRead.toString();
 110   }
 111
 112   /**
 113    * Apply a cigar operation to a segment of a read, accounting for any
 114    * insertions to the reference from other reads in the region
 115    *
 116    * @param el
 117    *          the current CIGAR element
 118    * @param next
 119    *          location in read to start applying CIGAR op
 120    * @param length
 121    *          length of CIGAR op
 122    * @param rec
 123    *          SAMRecord corresponding to read
 124    * @param it
 125    *          iterator over insertions which coincide with rec
 126    * @param override
 127    *          an optional <insertion position, length> which can override a
 128    *          <position,length> in it to change the length. Set to null if not
 129    *          used.
 130    * @param seq
 131    * @return
 132    */
 133   private String applyCigarOp(CigarElement el, int next,
 134           SAMRecord rec, Iterator<Map.Entry<Integer, Integer>> it,
 135           int[] override, SequenceI seq)
 136   {
 137     StringBuilder newRead = new StringBuilder();
 138     String read = rec.getReadString();
 139
 140     int length = el.getLength();
 141     int nextPos = next;
 142
 143     switch (el.getOperator())
 144     {
 145     case M:
 146     case X:
 147     case EQ:
 148       // matched or mismatched residues
 149       newRead = applyCigarConsumeReadAndRef(next, length, rec, it,
 150               override, seq);
 151       break;
 152     case N: // intron in RNA
 153     case D: // deletion
 154       int insertCount = 0;
 155       while (it.hasNext())
 156       {
 157         // just count all the insertions we need to add in the deletion/intron
 158         // region
 159         Entry<Integer, Integer> entry = it.next();
 160         insertCount += entry.getValue();
 161       }
 162       addGaps(newRead, length + insertCount);
 163       break;
 164     case I:
 165       // the reference sequence and other reads should have been gapped for
 166       // this insertion, so just add in the residues
 167       newRead.append(read.substring(nextPos, nextPos + length));
 168       seq.addSequenceFeature(new SequenceFeature("INSERTION", "",
 169               nextPos + 1, nextPos + length, "bamfile"));
 170       break;
 171     case S:
 172       // soft clipping - just skip this bit of the read
 173     case H:
 174       // hard clipping - this stretch will not appear in the read
 175     default:
 176       break;
 177     }
 178
 179     return newRead.toString();
 180   }
 181
 182   /**
 183    * Applies a CIGAR operation which consumes both read and reference i.e. M,=,X
 184    *
 185    * @param next
 186    *          next position in read
 187    * @param length
 188    *          length of cigar operation
 189    * @param rec
 190    *          SAMRecord to apply the CIGAR op to
 191    * @param it
 192    *          iterator over insertions in the CIGAR region
 193    * @param override
 194    *          insertions which have already been accounted for
 195    * @param seq
 196    * @return
 197    */
 198   private StringBuilder applyCigarConsumeReadAndRef(int next, int length,
 199           SAMRecord rec, Iterator<Map.Entry<Integer, Integer>> it,
 200           int[] override, SequenceI seq)
 201   {
 202     StringBuilder newRead = new StringBuilder();
 203     String read = rec.getReadString();
 204
 205     int nextPos = next; // location of next position in read
 206     while (nextPos < next + length)
 207     {
 208       int insertPos = next + length; // location of next insert in read (or end
 209                                      // of read)
 210       int insertLen = 0; // no of gaps to insert after insertPos
 211       if (it.hasNext())
 212       {
 213         // account for sequence already having insertion
 214         // if an override entry exists it gives the number of positions still to
 215         // be inserted at this point, after accounting for insertions driven by
 216         // the current read.
 217         // e.g. if the current read has 2 inserted positions here, and the
 218         // insertions list also has 2 inserted positions, the override will
 219         // cancel the insert for this read. If the insertions list had 3
 220         // insertions, the override would reduce this to 1 insertion.
 221         Entry<Integer, Integer> entry = it.next();
 222         int key = entry.getKey();
 223         insertLen = entry.getValue();
 224         if (override != null && key == override[0])
 225         {
 226           insertLen = override[1];
 227         }
 228         if (insertLen > 0)
 229         {
 230           insertPos = rec.getReadPositionAtReferencePosition(key) - 1;
 231         }
 232         else
 233         {
 234           insertPos = nextPos; // bugfix - trigger by override + insertions in
 235                                // same M stretch
 236         }
 237       }
 238       newRead.append(read.substring(nextPos, insertPos));
 239       nextPos = insertPos;
 240       addGaps(newRead, insertLen); // add insertions
 241     }
 242     return newRead;
 243   }
 244
 245   /**
 246    * Get list of positions inserted to the reference sequence. Note: assumes all
 247    * reads are on the same chromosome.
 248    *
 249    * @param it
 250    *          Iterator over all the SAMRecords
 251    * @return sorted map of insertion positions <position, insertion size>
 252    */
 253   /*
 254   1234567
 255   ABCDEFG insert 3,4
 256
 257   R1: insert 3 = gap before 3rd res
 258   AB.CDEFG
 259
 260   R2: insert 4 = gap before 4th res
 261   ABC.DEFG
 262
 263   R3: insert 3,4 = 2 gaps before 3rd res
 264   AB..CDEFG
 265
 266   => AB--C-DEFG
 267
 268   So need to store insertions as (pos, length) giving here (3,1),(4,1),(3,2) - then can sub longest length at position so (3,2), (4,1)
 269   Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added
 270
 271   */
 272   public SortedMap<Integer, Integer> getInsertions(Iterator<SAMRecord> it)
 273   {
 274     SortedMap<Integer, Integer> inserts = new TreeMap<>();
 275     while (it.hasNext())
 276     {
 277       // check each record for insertions in the CIGAR string
 278       SAMRecord rec = it.next();
 279       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 280               .iterator();
 281       int next = 1;
 282       while (cit.hasNext())
 283       {
 284         CigarElement el = cit.next();
 285         switch (el.getOperator())
 286         {
 287         case I:
 288           // add to insertions list, and move along the read
 289           // location is the start of next CIGAR segment
 290           // because getReferencePositionAtReadPosition returns 0 for read
 291           // positions which are not in the reference (like insertions)
 292           int refLocation = rec.getReferencePositionAtReadPosition(
 293                   next + el.getLength());
 294
 295           // if there's already an insertion at this location, keep the longest
 296           // insertion; if there's no insertion keep this one
 297           if (!inserts.containsKey(refLocation)
 298                   || (inserts.containsKey(refLocation)
 299                           && inserts.get(refLocation) < el.getLength()))
 300           {
 301             inserts.put(refLocation, el.getLength());
 302           }
 303           next += el.getLength();
 304           break;
 305         case M:
 306         case S:
 307           // match to reference, or soft clip, move along the read
 308           next += el.getLength();
 309           break;
 310         default:
 311           // deletions, introns etc don't consume any residues from the read
 312           break;
 313         }
 314       }
 315
 316     }
 317     return inserts;
 318   }
 319
 320   /**
 321    * Add sequence of gaps to a read
 322    *
 323    * @param read
 324    *          read to update
 325    * @param length
 326    *          number of gaps
 327    */
 328   private void addGaps(StringBuilder read, int length)
 329   {
 330     if (length > 0)
 331     {
 332       char[] gaps = new char[length];
 333       Arrays.fill(gaps, gapChar);
 334       read.append(gaps);
 335     }
 336   }
 337
 338   /**
 339    * Get the number of insertions before the start of an aligned read (i.e.
 340    * insertions in soft clipped region are counted too)
 341    *
 342    * @param rec
 343    *          the SAMRecord for the read
 344    * @param inserts
 345    *          the map of insertions
 346    * @return number of insertions before the read starts
 347    */
 348   private int countInsertionsBeforeRead(SAMRecord rec,
 349           SortedMap<Integer, Integer> inserts)
 350   {
 351     int gaps = 0;
 352
 353     // add in any insertion gaps before read
 354     // although we only need to get the submap from alignmentStart, there won't
 355     // be any insertions before that so we can just start at 0
 356     SortedMap<Integer, Integer> seqInserts = inserts.subMap(0,
 357             rec.getStart());
 358
 359     Iterator<Map.Entry<Integer, Integer>> it = seqInserts.entrySet()
 360             .iterator();
 361     while (it.hasNext())
 362     {
 363       Entry<Integer, Integer> entry = it.next();
 364       gaps += entry.getValue();
 365     }
 366     return gaps;
 367   }
 368 }