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