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;
 198     while (nextPos < next + length)
 199     {
 200       int insertPos = next + length;
 201       int insertLen = 0;
 202       if (it.hasNext())
 203       {
 204         // account for sequence already having insertion
 205         // if an override entry exists it gives the number of positions still to
 206         // be inserted at this point, after accounting for insertions driven by
 207         // the current read.
 208         // e.g. if the current read has 2 inserted positions here, and the
 209         // insertions list also has 2 inserted positions, the override will
 210         // cancel the insert for this read. If the insertions list had 3
 211         // insertions, the override would reduce this to 1 insertion.
 212         Entry<Integer, Integer> entry = it.next();
 213         int key = entry.getKey();
 214         insertLen = entry.getValue();
 215         if (override != null && key == override[0])
 216         {
 217           insertLen = override[1];
 218         }
 219         if (insertLen > 0)
 220         {
 221           insertPos = rec.getReadPositionAtReferencePosition(key) - 1;
 222         }
 223       }
 224       newRead.append(read.substring(nextPos, insertPos));
 225       nextPos = insertPos;
 226       addGaps(newRead, insertLen); // add insertions
 227     }
 228     return newRead;
 229   }
 230
 231   /**
 232    * Get list of positions inserted to the reference sequence. Note: assumes all
 233    * reads are on the same chromosome.
 234    *
 235    * @param it
 236    *          Iterator over all the SAMRecords
 237    * @return sorted map of insertion positions <position, insertion size>
 238    */
 239   /*
 240   1234567
 241   ABCDEFG insert 3,4
 242
 243   R1: insert 3 = gap before 3rd res
 244   AB.CDEFG
 245
 246   R2: insert 4 = gap before 4th res
 247   ABC.DEFG
 248
 249   R3: insert 3,4 = 2 gaps before 3rd res
 250   AB..CDEFG
 251
 252   => AB--C-DEFG
 253
 254   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)
 255   Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added
 256
 257   */
 258   public SortedMap<Integer, Integer> getInsertions(Iterator<SAMRecord> it)
 259   {
 260     SortedMap<Integer, Integer> inserts = new TreeMap<>();
 261     while (it.hasNext())
 262     {
 263       // check each record for insertions in the CIGAR string
 264       SAMRecord rec = it.next();
 265       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 266               .iterator();
 267       int next = 0;
 268       while (cit.hasNext())
 269       {
 270         CigarElement el = cit.next();
 271         switch (el.getOperator())
 272         {
 273         case I:
 274           // add to insertions list, and move along the read
 275           // location is 1 past the current position of next
 276           // (note if we try to retrieve +1 position from reference by calling
 277           // getReferencePositionAtReadPosition(next+1) we get 0 because it's an
 278           // insertion!)
 279           int refLocation = rec.getReferencePositionAtReadPosition(next)
 280                   + 1;
 281
 282           // if there's already an insertion at this location, keep the longest
 283           // insertion; if there's no insertion keep this one
 284           if (!inserts.containsKey(refLocation)
 285                   || (inserts.containsKey(refLocation)
 286                           && inserts.get(refLocation) < el.getLength()))
 287           {
 288             inserts.put(refLocation, el.getLength());
 289           }
 290           next += el.getLength();
 291           break;
 292         case M:
 293         case S:
 294           // match to reference, or soft clip, move along the read
 295           next += el.getLength();
 296           break;
 297         default:
 298           // deletions, introns etc don't consume any residues from the read
 299           break;
 300         }
 301       }
 302
 303     }
 304     return inserts;
 305   }
 306
 307   /**
 308    * Add sequence of gaps to a read
 309    *
 310    * @param read
 311    *          read to update
 312    * @param length
 313    *          number of gaps
 314    */
 315   private void addGaps(StringBuilder read, int length)
 316   {
 317     if (length > 0)
 318     {
 319       char[] gaps = new char[length];
 320       Arrays.fill(gaps, gapChar);
 321       read.append(gaps);
 322     }
 323   }
 324
 325   /**
 326    * Get the number of insertions before the start of an aligned read (i.e.
 327    * insertions in soft clipped region are counted too)
 328    *
 329    * @param rec
 330    *          the SAMRecord for the read
 331    * @param inserts
 332    *          the map of insertions
 333    * @return number of insertions before the read starts
 334    */
 335   private int countInsertionsBeforeRead(SAMRecord rec,
 336           SortedMap<Integer, Integer> inserts)
 337   {
 338     int gaps = 0;
 339
 340     // add in any insertion gaps before read
 341     // although we only need to get the submap from alignmentStart, there won't
 342     // be any insertions before that so we can just start at 0
 343     SortedMap<Integer, Integer> seqInserts = inserts.subMap(0,
 344             rec.getStart());
 345
 346     Iterator<Map.Entry<Integer, Integer>> it = seqInserts.entrySet()
 347             .iterator();
 348     while (it.hasNext())
 349     {
 350       Entry<Integer, Integer> entry = it.next();
 351       gaps += entry.getValue();
 352     }
 353     return gaps;
 354   }
 355 }