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