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