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