src/jalview/io/BamFile.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.io;
  22
  23 import jalview.datamodel.Sequence;
  24 import jalview.datamodel.SequenceI;
  25
  26 import java.io.File;
  27 import java.io.IOException;
  28 import java.util.Arrays;
  29 import java.util.Iterator;
  30 import java.util.Map;
  31 import java.util.Map.Entry;
  32 import java.util.SortedMap;
  33 import java.util.TreeMap;
  34
  35 import htsjdk.samtools.CigarElement;
  36 import htsjdk.samtools.SAMRecord;
  37 import htsjdk.samtools.SAMRecordIterator;
  38 import htsjdk.samtools.SamReader;
  39 import htsjdk.samtools.SamReaderFactory;
  40 import htsjdk.samtools.ValidationStringency;
  41
  42 public class BamFile extends AlignFile
  43 {
  44
  45   SamReader fileReader;
  46
  47   /**
  48    * Creates a new BamFile object.
  49    */
  50   public BamFile()
  51   {
  52   }
  53
  54   /**
  55    * Creates a new BamFile object.
  56    *
  57    * @param inFile
  58    *          DOCUMENT ME!
  59    * @param sourceType
  60    *          DOCUMENT ME!
  61    *
  62    * @throws IOException
  63    *           DOCUMENT ME!
  64    */
  65   public BamFile(String inFile, DataSourceType sourceType)
  66           throws IOException
  67   {
  68     super(inFile, sourceType);
  69     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  70             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  71                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  72             .validationStringency(ValidationStringency.SILENT);
  73     fileReader = factory.open(new File(inFile));
  74   }
  75
  76   public BamFile(FileParse source) throws IOException
  77   {
  78     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  79             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  80                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  81             .validationStringency(ValidationStringency.SILENT);
  82
  83     // File-based bam
  84     fileReader = factory.open(source.inFile);
  85     parse();
  86   }
  87
  88   @Override
  89   public String print(SequenceI[] seqs, boolean jvsuffix)
  90   {
  91     // TODO Auto-generated method stub
  92     return null;
  93   }
  94
  95   @Override
  96   public void parse() throws IOException
  97   {
  98     SAMRecordIterator it = fileReader.iterator();
  99     SortedMap<Integer,Integer> insertions = getInsertions(it);
 100     it.close();
 101
 102     it = fileReader.iterator();
 103     while (it.hasNext())
 104     {
 105       SAMRecord rec = it.next();
 106       String read = rec.getReadString();
 107       int start = rec.getStart();
 108       int end = rec.getEnd();
 109
 110       SequenceI seq = new Sequence(rec.getReadName(), rec.getReadString());
 111
 112       String cigarredRead = parseCigarToSequence(read, rec, insertions,
 113               '-');
 114
 115       SequenceI alsq = seq.deriveSequence();
 116       alsq.setSequence(cigarredRead);
 117       alsq.setStart(start);
 118       alsq.setEnd(end);
 119       seqs.add(alsq);
 120     }
 121
 122     // put insertions back into each sequence
 123     for (SequenceI seq : seqs)
 124     {
 125       int insertCount = 0;
 126       SortedMap<Integer, Integer> seqInserts = insertions
 127               .subMap(seq.getStart(),
 128               seq.getEnd()); // TODO start point should be start of alignment
 129                              // not 0
 130
 131       Iterator<Map.Entry<Integer, Integer>> iit = seqInserts.entrySet()
 132               .iterator();
 133       while (iit.hasNext())
 134       {
 135         // TODO account for sequence already having insertion, and INDELS
 136         Entry<Integer, Integer> entry = iit.next();
 137         int pos = entry.getKey();
 138         int length = entry.getValue();
 139         seq.insertCharAt(pos + insertCount, length, '-');
 140         insertCount++;
 141       }
 142     }
 143   }
 144
 145   /*
 146   1234567
 147   ABCDEFG insert 3,4
 148
 149   R1: insert 3 = gap before 3rd res
 150   AB.CDEFG
 151
 152   R2: insert 4 = gap before 4th res
 153   ABC.DEFG
 154
 155   R3: insert 3,4 = 2 gaps before 3rd res
 156   AB..CDEFG
 157
 158   => AB--C-DEFG
 159
 160   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)
 161   Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added
 162
 163   */
 164   /**
 165    * Apply the CIGAR string to a read sequence and return the updated read.
 166    *
 167    * @param read
 168    *          the read to update
 169    * @param rec
 170    *          the corresponding SAM record
 171    * @param gapChar
 172    *          gap character to use
 173    * @return string representing read with gaps, clipping etc applied
 174    */
 175   private String parseCigarToSequence(String read,
 176           SAMRecord rec, SortedMap<Integer, Integer> insertions,
 177           char gapChar)
 178   {
 179     Iterator<CigarElement> it = rec.getCigar().getCigarElements()
 180             .iterator();
 181
 182     StringBuilder newRead = new StringBuilder();
 183     int next = 0;
 184
 185     // pad with spaces before read
 186     // first count gaps needed to pad to reference position
 187     int gaps = rec.getReferencePositionAtReadPosition(next + 1);
 188     // now count gaps to pad for insertions in other reads
 189     int insertCount = countInsertionsBeforeRead(rec, insertions);
 190     addGaps(newRead, gaps + insertCount, gapChar);
 191
 192     SortedMap<Integer, Integer> seqInserts = insertions
 193             .subMap(rec.getStart(), rec.getEnd() + 1);
 194     // TODO start point should be start of alignment not 0
 195
 196     Iterator<Map.Entry<Integer, Integer>> iit = seqInserts.entrySet()
 197             .iterator();
 198     if (iit.hasNext())
 199     {
 200       // TODO account for sequence already having insertion, and INDELS
 201       Entry<Integer, Integer> entry = iit.next();
 202       int insertPos = entry.getKey();
 203       int insertLen = entry.getValue();
 204       // seq.insertCharAt(pos + insertCount, length, '-');
 205       // insertCount++;
 206     }
 207
 208     while (it.hasNext())
 209     {
 210       CigarElement el = it.next();
 211       int length = el.getLength();
 212
 213
 214       switch (el.getOperator())
 215       {
 216       case M:
 217       case X:
 218       case EQ:
 219         // matched or mismatched residues
 220         newRead.append(read.substring(next, next + length));
 221         next += length;
 222         break;
 223       case N: // intron in RNA
 224       case D: // deletion
 225         addGaps(newRead, length, gapChar);
 226         break;
 227       case S:
 228         // soft clipping - just skip this bit of the read
 229         // do nothing
 230
 231         // work out how many gaps we need before the start of the soft clip -
 232         // don't do this at the end of the read!
 233         if (next == 0) // at start of read
 234         {
 235           addGaps(newRead, rec.getUnclippedStart(), gapChar);
 236         }
 237
 238         newRead.append(read.substring(next, next + length).toLowerCase());
 239         next += length;
 240         break;
 241       case I:
 242         // the reference sequence and other reads should have been gapped for
 243         // this insertion, so just add in the residues
 244         newRead.append(read.substring(next, next + length));
 245         next += length;
 246         break;
 247       case H:
 248         // hard clipping - this stretch will not appear in the read
 249         break;
 250       default:
 251         break;
 252       }
 253
 254
 255     }
 256     return newRead.toString();
 257   }
 258
 259   private String applyCigar(String read, CigarElement el, char gapChar,
 260           int next, int length, int softStart)
 261   {
 262     StringBuilder newRead = new StringBuilder();
 263     switch (el.getOperator())
 264     {
 265     case M:
 266     case X:
 267     case EQ:
 268       // matched or mismatched residues
 269       newRead.append(read.substring(next, next + length));
 270       next += length;
 271       break;
 272     case N: // intron in RNA
 273     case D: // deletion
 274       addGaps(newRead, length, gapChar);
 275       break;
 276     case S:
 277       // soft clipping - just skip this bit of the read
 278       // do nothing
 279
 280       // work out how many gaps we need before the start of the soft clip -
 281       // don't do this at the end of the read!
 282       if (next == 0) // at start of read
 283       {
 284         addGaps(newRead, softStart, gapChar);
 285       }
 286
 287       newRead.append(read.substring(next, next + length).toLowerCase());
 288       next += length;
 289       break;
 290     case I:
 291       // the reference sequence and other reads should have been gapped for
 292       // this insertion, so just add in the residues
 293       newRead.append(read.substring(next, next + length));
 294       next += length;
 295       break;
 296     case H:
 297       // hard clipping - this stretch will not appear in the read
 298       break;
 299     default:
 300       break;
 301     }
 302     return newRead.toString();
 303   }
 304
 305   /**
 306    * Get list of positions inserted to the reference sequence
 307    *
 308    * @param it
 309    * @return
 310    */
 311   private SortedMap<Integer, Integer> getInsertions(Iterator<SAMRecord> it)
 312   {
 313     SortedMap<Integer, Integer> insertions = new TreeMap<>();
 314     while (it.hasNext())
 315     {
 316       // check each record for insertions in the CIGAR string
 317       SAMRecord rec = it.next();
 318       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 319               .iterator();
 320       int next = 0;
 321       while (cit.hasNext())
 322       {
 323         CigarElement el = cit.next();
 324         switch (el.getOperator())
 325         {
 326         case I:
 327           // add to insertions list, and move along the read
 328           int refLocation = rec.getReferencePositionAtReadPosition(next);
 329
 330           // if there's already an insertion at this location, keep the longest
 331           // insertion; if there's no insertion keep this one
 332           if (!insertions.containsKey(refLocation)
 333                   || (insertions.containsKey(refLocation)
 334                           && insertions.get(refLocation) < el.getLength()))
 335           {
 336             insertions.put(refLocation, el.getLength());
 337           }
 338           next += el.getLength();
 339           break;
 340         case M:
 341           // match to reference, move along the read
 342           next += el.getLength();
 343           break;
 344         default:
 345           // deletions, introns etc don't consume any residues from the read
 346           break;
 347         }
 348       }
 349
 350     }
 351     return insertions;
 352   }
 353
 354   /**
 355    * Add sequence of gaps to a read
 356    *
 357    * @param read
 358    *          read to update
 359    * @param length
 360    *          number of gaps
 361    * @param gapChar
 362    *          gap character to use
 363    */
 364   private void addGaps(StringBuilder read, int length, char gapChar)
 365   {
 366     if (length > 0)
 367     {
 368       char[] gaps = new char[length];
 369       Arrays.fill(gaps, gapChar);
 370       read.append(gaps);
 371     }
 372   }
 373
 374   private int countInsertionsBeforeRead(SAMRecord rec,
 375           SortedMap<Integer, Integer> insertions)
 376   {
 377     int gaps = 0;
 378
 379     // add in any insertion gaps before read
 380     // TODO start point should be start of alignment not 0
 381     SortedMap<Integer, Integer> seqInserts = insertions.subMap(0,
 382             rec.getStart());
 383
 384     Iterator<Map.Entry<Integer, Integer>> it = seqInserts.entrySet()
 385             .iterator();
 386     while (it.hasNext())
 387     {
 388       Entry<Integer, Integer> entry = it.next();
 389       gaps += entry.getValue();
 390     }
 391     return gaps;
 392   }
 393
 394 }