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       seq.addSequenceFeature(new SequenceFeature("SOFTCLIP", "",
 174               nextPos + 1, nextPos + length, "bamfile"));
 175       break;
 176     case H:
 177       // hard clipping - this stretch will not appear in the read
 178       seq.addSequenceFeature(new SequenceFeature("HARDCLIP", "",
 179               nextPos + 1, nextPos + length, "bamfile"));
 180       break;
 181     default:
 182       break;
 183     }
 184
 185     return newRead.toString();
 186   }
 187
 188   /**
 189    * Applies a CIGAR operation which consumes both read and reference i.e. M,=,X
 190    *
 191    * @param next
 192    *          next position in read
 193    * @param length
 194    *          length of cigar operation
 195    * @param rec
 196    *          SAMRecord to apply the CIGAR op to
 197    * @param it
 198    *          iterator over insertions in the CIGAR region
 199    * @param override
 200    *          insertions which have already been accounted for
 201    * @param seq
 202    * @return
 203    */
 204   private StringBuilder applyCigarConsumeReadAndRef(int next, int length,
 205           SAMRecord rec, Iterator<Map.Entry<Integer, Integer>> it,
 206           int[] override, SequenceI seq)
 207   {
 208     StringBuilder newRead = new StringBuilder();
 209     String read = rec.getReadString();
 210
 211     int nextPos = next; // location of next position in read
 212     while (nextPos < next + length)
 213     {
 214       int insertPos = next + length; // location of next insert in read (or end
 215                                      // of read)
 216       int insertLen = 0; // no of gaps to insert after insertPos
 217       if (it.hasNext())
 218       {
 219         // account for sequence already having insertion
 220         // if an override entry exists it gives the number of positions still to
 221         // be inserted at this point, after accounting for insertions driven by
 222         // the current read.
 223         // e.g. if the current read has 2 inserted positions here, and the
 224         // insertions list also has 2 inserted positions, the override will
 225         // cancel the insert for this read. If the insertions list had 3
 226         // insertions, the override would reduce this to 1 insertion.
 227         Entry<Integer, Integer> entry = it.next();
 228         int key = entry.getKey();
 229         insertLen = entry.getValue();
 230         if (override != null && key == override[0])
 231         {
 232           insertLen = override[1];
 233         }
 234         if (insertLen > 0)
 235         {
 236           insertPos = rec.getReadPositionAtReferencePosition(key) - 1;
 237         }
 238         else
 239         {
 240           insertPos = nextPos; // bugfix - trigger by override + insertions in
 241                                // same M stretch
 242         }
 243       }
 244       newRead.append(read.substring(nextPos, insertPos));
 245       nextPos = insertPos;
 246       addGaps(newRead, insertLen); // add insertions
 247     }
 248     return newRead;
 249   }
 250
 251   /**
 252    * Get list of positions inserted to the reference sequence. Note: assumes all
 253    * reads are on the same chromosome.
 254    *
 255    * @param it
 256    *          Iterator over all the SAMRecords
 257    * @return sorted map of insertion positions <position, insertion size>
 258    */
 259   /*
 260   1234567
 261   ABCDEFG insert 3,4
 262
 263   R1: insert 3 = gap before 3rd res
 264   AB.CDEFG
 265
 266   R2: insert 4 = gap before 4th res
 267   ABC.DEFG
 268
 269   R3: insert 3,4 = 2 gaps before 3rd res
 270   AB..CDEFG
 271
 272   => AB--C-DEFG
 273
 274   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)
 275   Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added
 276
 277   */
 278   public SortedMap<Integer, Integer>[] getInsertions(Iterator<SAMRecord> it)
 279   {
 280     return getInsertions(it, null);
 281   }
 282
 283   public int firstAlColumn = -1, firstRposCol = -1;
 284
 285   public SortedMap<Integer, Integer>[] getInsertions(Iterator<SAMRecord> it,
 286           Range[] extent)
 287   {
 288     SortedMap<Integer, Integer> sinserts[] = new SortedMap[] {
 289         new TreeMap<>(), new TreeMap<>() };
 290     Range xten = extent != null && extent[0] != null ? extent[0] : null;
 291     do
 292     {
 293       // check extent of read
 294       SAMRecord rec = it.next();
 295       if (extent != null)
 296       {
 297
 298         int nstart = rec.getAlignmentStart();
 299         if (nstart != 0)
 300         {
 301           nstart = rec.getReferencePositionAtReadPosition(nstart);
 302         }
 303         int nend = rec.getAlignmentEnd();
 304         if (nend != 0)
 305         {
 306           nend = rec.getReferencePositionAtReadPosition(nend);
 307         }
 308
 309         xten = new Range(
 310                 nstart <= 0 ? xten.start : Math.min(nstart, xten.start),
 311                 nend == 0 ? xten.end : Math.max(nend, xten.end));
 312       }
 313
 314       // check each record for insertions in the CIGAR string
 315
 316       // pick the insert map to update
 317       int insmap = rec.getReadNegativeStrandFlag() ? 1 : 0;
 318       SortedMap<Integer, Integer> inserts = sinserts[insmap];
 319
 320       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 321               .iterator();
 322       int refLocation = rec.getAlignmentStart();
 323       int rloc = 0;
 324       int alcol = 0;
 325       int alpos = 0;
 326       while (cit.hasNext())
 327       {
 328         CigarElement el = cit.next();
 329         switch (el.getOperator())
 330         {
 331         case S:
 332         case H:
 333           // ignore soft/hardclipped
 334
 335           break;
 336         default:
 337           if (el.getOperator().consumesReadBases()
 338                   && !el.getOperator().consumesReferenceBases())
 339           {
 340             // add to insertions list, and move along the read
 341             // location is the start of next CIGAR segment
 342             // because getReferencePositionAtReadPosition returns 0 for read
 343             // positions which are not in the reference (like insertions)
 344
 345             // if there's already an insertion at this location, keep the
 346             // longest
 347             // insertion; if there's no insertion keep this one
 348             if (!inserts.containsKey(refLocation)
 349                     || (inserts.containsKey(refLocation)
 350                             && inserts.get(refLocation) < el.getLength()))
 351             {
 352               inserts.put(refLocation, el.getLength());
 353             }
 354             break;
 355           }
 356           if (el.getOperator().consumesReferenceBases())
 357           {
 358             if (el.getOperator().consumesReadBases() && alcol == 0
 359                     && refLocation + 1 > extent[0].start)
 360             {
 361               // first match position : last rloc + first match op
 362               alcol = rloc + 1;
 363               alpos = rec.getReferencePositionAtReadPosition(rloc + 1);
 364             }
 365             refLocation += el.getLength();
 366           }
 367
 368         }
 369         if (el.getOperator().consumesReadBases()
 370                 || el.getOperator().consumesReferenceBases())
 371         {
 372           rloc += el.getLength();
 373         }
 374       }
 375       // Wrong: hack that fails to relocate reference to correct place in
 376       // sequence
 377       if (firstAlColumn < 0)
 378       {
 379         firstAlColumn = alcol;
 380         firstRposCol = alpos;
 381       }
 382     } while (it.hasNext());
 383     if (xten != null)
 384     {
 385       extent[0] = xten;
 386     }
 387     return sinserts;
 388   }
 389
 390   /**
 391    * Add sequence of gaps to a read
 392    *
 393    * @param read
 394    *          read to update
 395    * @param length
 396    *          number of gaps
 397    */
 398   private void addGaps(StringBuilder read, int length)
 399   {
 400     if (length > 0)
 401     {
 402       char[] gaps = new char[length];
 403       Arrays.fill(gaps, gapChar);
 404       read.append(gaps);
 405     }
 406   }
 407
 408   /**
 409    * Get the number of insertions before the start of an aligned read (i.e.
 410    * insertions in soft clipped region are counted too)
 411    *
 412    * @param rec
 413    *          the SAMRecord for the read
 414    * @param inserts
 415    *          the map of insertions
 416    * @return number of insertions before the read starts
 417    */
 418   private int countInsertionsBeforeRead(SAMRecord rec,
 419           SortedMap<Integer, Integer> inserts)
 420   {
 421     int gaps = 0;
 422
 423     // add in any insertion gaps before read
 424     // although we only need to get the submap from alignmentStart, there won't
 425     // be any insertions before that so we can just start at 0
 426     SortedMap<Integer, Integer> seqInserts = inserts.subMap(0,
 427             rec.getStart());
 428
 429     Iterator<Map.Entry<Integer, Integer>> it = seqInserts.entrySet()
 430             .iterator();
 431     while (it.hasNext())
 432     {
 433       Entry<Integer, Integer> entry = it.next();
 434       gaps += entry.getValue();
 435     }
 436     return gaps;
 437   }
 438 }