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.SAMFormatException;
  13 import htsjdk.samtools.SAMRecord;
  14 import jalview.bin.Cache;
  15
  16 public class CigarParser
  17 {
  18   private final char gapChar;
  19
  20   public CigarParser(char gap)
  21   {
  22     gapChar = gap;
  23   }
  24
  25   /**
  26    * Apply the CIGAR string to a read sequence and return the updated read.
  27    *
  28    * @param rec
  29    *          the SAM record for the read
  30    * @param insertions
  31    *          a map of inserted positions and lengths
  32    * @param alignmentStart
  33    *          start of alignment to be used as offset when padding reads with
  34    *          gaps
  35    * @param seq
  36    *          sequence to be annotated with inserts
  37    * @return string representing read with gaps, clipping etc applied
  38    */
  39   public String parseCigarToSequence(SAMRecord rec,
  40           SortedMap<Integer, Integer> insertions, int alignmentStart,
  41           SequenceI seq)
  42   {
  43     StringBuilder newRead = new StringBuilder();
  44     Iterator<CigarElement> it = rec.getCigar().getCigarElements()
  45             .iterator();
  46
  47     int next = 0;
  48     int refnext = 0;
  49
  50     // pad with spaces before read
  51     // first count gaps needed to pad to reference position
  52     // NB position is 1-based so number of gaps = pos-1
  53     int gaps = rec.getStart() - alignmentStart;
  54
  55     // now count gaps to pad for insertions in other reads
  56     int insertCount = countInsertionsBeforeRead(rec, insertions);
  57     addGaps(newRead, gaps + insertCount);
  58
  59     int lastinserts = 0;
  60     while (it.hasNext())
  61     {
  62       CigarElement el = it.next();
  63       int length = el.getLength();
  64
  65       // which insertions coincide with current CIGAR operation?
  66       SortedMap<Integer, Integer> seqInserts = insertions.subMap(
  67               rec.getStart() + refnext, rec.getStart() + refnext + length);
  68
  69       int[] override = null;
  70       if (lastinserts > 0)
  71       {
  72         if (!seqInserts.containsKey(rec.getStart() + refnext))
  73         {
  74           // ERROR
  75           int pos = rec.getStart() + refnext;
  76           System.out.println(
  77                   "Read " + rec.getReadName() + " has an insertion at "
  78                           + pos + " which is not in seqInserts");
  79         }
  80         else
  81         {
  82           // we just inserted something
  83           // remove these inserts now - should be very first entry
  84           int count = seqInserts.get(rec.getStart() + refnext);
  85           override = new int[] { rec.getStart() + refnext,
  86               count - lastinserts };
  87           lastinserts = 0;
  88         }
  89       }
  90
  91       Iterator<Map.Entry<Integer, Integer>> iit = seqInserts.entrySet()
  92               .iterator();
  93
  94       String nextSegment = applyCigarOp(el, next, rec, iit, override, seq);
  95       newRead.append(nextSegment);
  96
  97       if (el.getOperator().consumesReferenceBases())
  98       {
  99         refnext += length;
 100       }
 101       if (el.getOperator().consumesReadBases())
 102       {
 103         next += length;
 104       }
 105       if (el.getOperator().equals(CigarOperator.I))
 106       {
 107         // count how many insertions we made
 108         lastinserts = length;
 109       }
 110     }
 111     return newRead.toString();
 112   }
 113
 114   /**
 115    * Apply a cigar operation to a segment of a read, accounting for any
 116    * insertions to the reference from other reads in the region
 117    *
 118    * @param el
 119    *          the current CIGAR element
 120    * @param next
 121    *          location in read to start applying CIGAR op
 122    * @param length
 123    *          length of CIGAR op
 124    * @param rec
 125    *          SAMRecord corresponding to read
 126    * @param it
 127    *          iterator over insertions which coincide with rec
 128    * @param override
 129    *          an optional <insertion position, length> which can override a
 130    *          <position,length> in it to change the length. Set to null if not
 131    *          used.
 132    * @param seq
 133    * @return
 134    */
 135   private String applyCigarOp(CigarElement el, int next,
 136           SAMRecord rec, Iterator<Map.Entry<Integer, Integer>> it,
 137           int[] override, SequenceI seq)
 138   {
 139     StringBuilder newRead = new StringBuilder();
 140     String read = rec.getReadString();
 141
 142     int length = el.getLength();
 143     int nextPos = next;
 144
 145     switch (el.getOperator())
 146     {
 147     case M:
 148     case X:
 149     case EQ:
 150       // matched or mismatched residues
 151       newRead = applyCigarConsumeReadAndRef(next, length, rec, it,
 152               override, seq);
 153       break;
 154     case N: // intron in RNA
 155     case D: // deletion
 156       int insertCount = 0;
 157       while (it.hasNext())
 158       {
 159         // just count all the insertions we need to add in the deletion/intron
 160         // region
 161         Entry<Integer, Integer> entry = it.next();
 162         insertCount += entry.getValue();
 163       }
 164       addGaps(newRead, length + insertCount);
 165       break;
 166     case I:
 167       // the reference sequence and other reads should have been gapped for
 168       // this insertion, so just add in the residues
 169       newRead.append(read.substring(nextPos, nextPos + length));
 170       seq.addSequenceFeature(new SequenceFeature("INSERTION", "",
 171               nextPos + 1, nextPos + length, "bamfile"));
 172       break;
 173     case S:
 174       // soft clipping - just skip this bit of the read
 175       seq.addSequenceFeature(new SequenceFeature("SOFTCLIP", "",
 176               nextPos + 1, nextPos + length, "bamfile"));
 177       break;
 178     case H:
 179       // hard clipping - this stretch will not appear in the read
 180       seq.addSequenceFeature(new SequenceFeature("HARDCLIP", "",
 181               nextPos + 1, nextPos + length, "bamfile"));
 182       break;
 183     default:
 184       break;
 185     }
 186
 187     return newRead.toString();
 188   }
 189
 190   /**
 191    * Applies a CIGAR operation which consumes both read and reference i.e. M,=,X
 192    *
 193    * @param next
 194    *          next position in read
 195    * @param length
 196    *          length of cigar operation
 197    * @param rec
 198    *          SAMRecord to apply the CIGAR op to
 199    * @param it
 200    *          iterator over insertions in the CIGAR region
 201    * @param override
 202    *          insertions which have already been accounted for
 203    * @param seq
 204    * @return
 205    */
 206   private StringBuilder applyCigarConsumeReadAndRef(int next, int length,
 207           SAMRecord rec, Iterator<Map.Entry<Integer, Integer>> it,
 208           int[] override, SequenceI seq)
 209   {
 210     StringBuilder newRead = new StringBuilder();
 211     String read = rec.getReadString();
 212
 213     int nextPos = next; // location of next position in read
 214     while (nextPos < next + length)
 215     {
 216       int insertPos = next + length; // location of next insert in read (or end
 217                                      // of read)
 218       int insertLen = 0; // no of gaps to insert after insertPos
 219       if (it.hasNext())
 220       {
 221         // account for sequence already having insertion
 222         // if an override entry exists it gives the number of positions still to
 223         // be inserted at this point, after accounting for insertions driven by
 224         // the current read.
 225         // e.g. if the current read has 2 inserted positions here, and the
 226         // insertions list also has 2 inserted positions, the override will
 227         // cancel the insert for this read. If the insertions list had 3
 228         // insertions, the override would reduce this to 1 insertion.
 229         Entry<Integer, Integer> entry = it.next();
 230         int key = entry.getKey();
 231         insertLen = entry.getValue();
 232         if (override != null && key == override[0])
 233         {
 234           insertLen = override[1];
 235         }
 236         if (insertLen > 0)
 237         {
 238           insertPos = rec.getReadPositionAtReferencePosition(key) - 1;
 239         }
 240         else
 241         {
 242           insertPos = nextPos; // bugfix - trigger by override + insertions in
 243                                // same M stretch
 244         }
 245       }
 246       newRead.append(read.substring(nextPos, insertPos));
 247       nextPos = insertPos;
 248       addGaps(newRead, insertLen); // add insertions
 249     }
 250     return newRead;
 251   }
 252
 253   /**
 254    * Get list of positions inserted to the reference sequence. Note: assumes all
 255    * reads are on the same chromosome.
 256    *
 257    * @param it
 258    *          Iterator over all the SAMRecords
 259    * @return sorted map of insertion positions <position, insertion size>
 260    */
 261   /*
 262   1234567
 263   ABCDEFG insert 3,4
 264
 265   R1: insert 3 = gap before 3rd res
 266   AB.CDEFG
 267
 268   R2: insert 4 = gap before 4th res
 269   ABC.DEFG
 270
 271   R3: insert 3,4 = 2 gaps before 3rd res
 272   AB..CDEFG
 273
 274   => AB--C-DEFG
 275
 276   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)
 277   Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added
 278
 279   */
 280   public SortedMap<Integer, Integer>[] getInsertions(Iterator<SAMRecord> it)
 281   {
 282     return getInsertions(it, null);
 283   }
 284
 285   public int firstAlColumn = -1, firstRposCol = -1;
 286
 287   public SortedMap<Integer, Integer>[] getInsertions(Iterator<SAMRecord> it,
 288           Range[] extent)
 289   {
 290     SortedMap<Integer, Integer> sinserts[] = new SortedMap[] {
 291         new TreeMap<>(), new TreeMap<>() };
 292     Range xten = extent != null && extent[0] != null ? extent[0] : null;
 293     do
 294     {
 295       // check extent of read
 296       SAMRecord rec = null;
 297       try {
 298         rec = it.next();
 299       } catch (SAMFormatException ex)
 300       {
 301         Cache.log.info("Bailing on parsing SAM File - see error below",ex);
 302         break;
 303       }
 304       if (extent != null)
 305       {
 306
 307         int nstart = rec.getAlignmentStart();
 308         if (nstart != 0)
 309         {
 310           nstart = rec.getReferencePositionAtReadPosition(nstart);
 311         }
 312         int nend = rec.getAlignmentEnd();
 313         if (nend != 0)
 314         {
 315           nend = rec.getReferencePositionAtReadPosition(nend);
 316         }
 317
 318         xten = new Range(
 319                 nstart <= 0 ? xten.start : Math.min(nstart, xten.start),
 320                 nend == 0 ? xten.end : Math.max(nend, xten.end));
 321       }
 322
 323       // check each record for insertions in the CIGAR string
 324
 325       // pick the insert map to update
 326       int insmap = rec.getReadNegativeStrandFlag() ? 1 : 0;
 327       SortedMap<Integer, Integer> inserts = sinserts[insmap];
 328
 329       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 330               .iterator();
 331       int refLocation = rec.getAlignmentStart();
 332       int rloc = 0;
 333       int alcol = 0;
 334       int alpos = 0;
 335       while (cit.hasNext())
 336       {
 337         CigarElement el = cit.next();
 338         switch (el.getOperator())
 339         {
 340         case S:
 341         case H:
 342           // ignore soft/hardclipped
 343
 344           break;
 345         default:
 346           if (el.getOperator().consumesReadBases()
 347                   && !el.getOperator().consumesReferenceBases())
 348           {
 349             // add to insertions list, and move along the read
 350             // location is the start of next CIGAR segment
 351             // because getReferencePositionAtReadPosition returns 0 for read
 352             // positions which are not in the reference (like insertions)
 353
 354             // if there's already an insertion at this location, keep the
 355             // longest
 356             // insertion; if there's no insertion keep this one
 357             if (!inserts.containsKey(refLocation)
 358                     || (inserts.containsKey(refLocation)
 359                             && inserts.get(refLocation) < el.getLength()))
 360             {
 361               inserts.put(refLocation, el.getLength());
 362             }
 363             break;
 364           }
 365           if (el.getOperator().consumesReferenceBases())
 366           {
 367             if (el.getOperator().consumesReadBases() && alcol == 0
 368                     && refLocation + 1 > extent[0].start)
 369             {
 370               // first match position : last rloc + first match op
 371               alcol = rloc + 1;
 372               alpos = rec.getReferencePositionAtReadPosition(rloc + 1);
 373             }
 374             refLocation += el.getLength();
 375           }
 376
 377         }
 378         if (el.getOperator().consumesReadBases()
 379                 || el.getOperator().consumesReferenceBases())
 380         {
 381           rloc += el.getLength();
 382         }
 383       }
 384       // Wrong: hack that fails to relocate reference to correct place in
 385       // sequence
 386       if (firstAlColumn < 0)
 387       {
 388         firstAlColumn = alcol;
 389         firstRposCol = alpos;
 390       }
 391     } while (it.hasNext());
 392     if (xten != null)
 393     {
 394       extent[0] = xten;
 395     }
 396     return sinserts;
 397   }
 398
 399   /**
 400    * Add sequence of gaps to a read
 401    *
 402    * @param read
 403    *          read to update
 404    * @param length
 405    *          number of gaps
 406    */
 407   private void addGaps(StringBuilder read, int length)
 408   {
 409     if (length > 0)
 410     {
 411       char[] gaps = new char[length];
 412       Arrays.fill(gaps, gapChar);
 413       read.append(gaps);
 414     }
 415   }
 416
 417   /**
 418    * Get the number of insertions before the start of an aligned read (i.e.
 419    * insertions in soft clipped region are counted too)
 420    *
 421    * @param rec
 422    *          the SAMRecord for the read
 423    * @param inserts
 424    *          the map of insertions
 425    * @return number of insertions before the read starts
 426    */
 427   private int countInsertionsBeforeRead(SAMRecord rec,
 428           SortedMap<Integer, Integer> inserts)
 429   {
 430     int gaps = 0;
 431
 432     // add in any insertion gaps before read
 433     // although we only need to get the submap from alignmentStart, there won't
 434     // be any insertions before that so we can just start at 0
 435     SortedMap<Integer, Integer> seqInserts = inserts.subMap(0,
 436             rec.getStart());
 437
 438     Iterator<Map.Entry<Integer, Integer>> it = seqInserts.entrySet()
 439             .iterator();
 440     while (it.hasNext())
 441     {
 442       Entry<Integer, Integer> entry = it.next();
 443       gaps += entry.getValue();
 444     }
 445     return gaps;
 446   }
 447 }