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.SortedMap;
  31 import java.util.TreeMap;
  32
  33 import htsjdk.samtools.CigarElement;
  34 import htsjdk.samtools.SAMRecord;
  35 import htsjdk.samtools.SAMRecordIterator;
  36 import htsjdk.samtools.SamReader;
  37 import htsjdk.samtools.SamReaderFactory;
  38 import htsjdk.samtools.ValidationStringency;
  39
  40 public class BamFile extends AlignFile
  41 {
  42
  43   SamReader fileReader;
  44
  45   /**
  46    * Creates a new BamFile object.
  47    */
  48   public BamFile()
  49   {
  50   }
  51
  52   /**
  53    * Creates a new BamFile object.
  54    *
  55    * @param inFile
  56    *          DOCUMENT ME!
  57    * @param sourceType
  58    *          DOCUMENT ME!
  59    *
  60    * @throws IOException
  61    *           DOCUMENT ME!
  62    */
  63   public BamFile(String inFile, DataSourceType sourceType)
  64           throws IOException
  65   {
  66     super(inFile, sourceType);
  67     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  68             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  69                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  70             .validationStringency(ValidationStringency.SILENT);
  71     fileReader = factory.open(new File(inFile));
  72   }
  73
  74   public BamFile(FileParse source) throws IOException
  75   {
  76     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  77             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  78                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  79             .validationStringency(ValidationStringency.SILENT);
  80
  81     // File-based bam
  82     fileReader = factory.open(source.inFile);
  83     parse();
  84   }
  85
  86   @Override
  87   public String print(SequenceI[] seqs, boolean jvsuffix)
  88   {
  89     // TODO Auto-generated method stub
  90     return null;
  91   }
  92
  93   @Override
  94   public void parse() throws IOException
  95   {
  96     SAMRecordIterator it = fileReader.iterator();
  97     SortedMap<Integer,Integer> insertions = getInsertions(it);
  98     it.close();
  99
 100     it = fileReader.iterator();
 101     while (it.hasNext())
 102     {
 103       SAMRecord rec = it.next();
 104       String read = rec.getReadString();
 105       int start = rec.getStart();
 106       int end = rec.getEnd();
 107
 108       SequenceI seq = new Sequence(rec.getReadName(), rec.getReadString());
 109
 110       String cigarredRead = parseCigarToSequence(read, rec, '-');
 111
 112       SequenceI alsq = seq.deriveSequence();
 113       alsq.setSequence(cigarredRead);
 114       alsq.setStart(start);
 115       alsq.setEnd(end);
 116       seqs.add(alsq);
 117     }
 118
 119     for (SequenceI seq : seqs)
 120     {
 121       int insertCount = 0;
 122       SortedMap<Integer, Integer> seqInserts = insertions.subMap(0,
 123               seq.getEnd()); // TODO start point should be start of alignment
 124                              // not 0
 125       seqInserts.forEach(action);
 126       while ((nextInsertion != -1) && (nextInsertion < seq.getEnd()))
 127       {
 128         seq.insertCharsAt(nextInsertion + insertCount, '-');
 129         // nextInsertion = insertions.nextSetBit(nextInsertion + 1);
 130         insertCount++;
 131       }
 132     }
 133   }
 134
 135   /*
 136   1234567
 137   ABCDEFG insert 3,4
 138
 139   R1: insert 3 = gap before 3rd res
 140   AB.CDEFG
 141
 142   R2: insert 4 = gap before 4th res
 143   ABC.DEFG
 144
 145   R3: insert 3,4 = 2 gaps before 3rd res
 146   AB..CDEFG
 147
 148   => AB--C-DEFG
 149
 150   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)
 151   Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added
 152
 153   */
 154   /**
 155    * Apply the CIGAR string to a read sequence and return the updated read.
 156    *
 157    * @param read
 158    *          the read to update
 159    * @param rec
 160    *          the corresponding SAM record
 161    * @param gapChar
 162    *          gap character to use
 163    * @return string representing read with gaps, clipping etc applied
 164    */
 165   private String parseCigarToSequence(String read,
 166           SAMRecord rec, char gapChar)
 167   {
 168     Iterator<CigarElement> it = rec.getCigar().getCigarElements()
 169             .iterator();
 170
 171     StringBuilder newRead = new StringBuilder();
 172     int next = 0;
 173     int ap = 0; // position of query region qi.start;
 174
 175     // pad with spaces before read
 176     for (int alp = rec
 177             .getReferencePositionAtReadPosition(next + 1); ap < alp;)
 178     {
 179       newRead.append(" ");
 180       ap++;
 181     }
 182
 183     while (it.hasNext())
 184     {
 185       CigarElement el = it.next();
 186       int length = el.getLength();
 187       char[] gaps;
 188
 189       switch (el.getOperator())
 190       {
 191       case M:
 192         // matched residues
 193         newRead.append(
 194                 read.substring(next, next + length));
 195         next += length;
 196         break;
 197       case N: // intron in RNA
 198       case D: // deletion
 199         // add gaps
 200         gaps = new char[length];
 201         Arrays.fill(gaps, gapChar);
 202         newRead.append(gaps);
 203         break;
 204       case S:
 205         // soft clipping - just skip this bit of the read
 206         // do nothing
 207
 208         // work out how many gaps we need before the start of the soft clip -
 209         // don't do this at the end of the read!
 210         if (next == 0) // at start of read
 211         {
 212           int numgaps = rec.getUnclippedStart();
 213           gaps = new char[numgaps];
 214           Arrays.fill(gaps, ' ');
 215           newRead.append(gaps);
 216         }
 217
 218         newRead.append(
 219                 read.substring(next, next + length).toLowerCase());
 220         next += length;
 221         break;
 222       case I:
 223         // the reference sequence and other reads should have been gapped for
 224         // this insertion, so just add in the residues
 225         newRead.append(read.substring(next, next + length));
 226         next += length;
 227         break;
 228       case H:
 229         // hard clipping - this stretch will not appear in the read
 230         break;
 231       default:
 232         // P, X EQ don't know what to do with these
 233         break;
 234       }
 235
 236     }
 237     return newRead.toString();
 238   }
 239
 240   /**
 241    * Get list of positions inserted to the reference sequence
 242    *
 243    * @param it
 244    * @return
 245    */
 246   private SortedMap<Integer, Integer> getInsertions(Iterator<SAMRecord> it)
 247   {
 248     SortedMap<Integer, Integer> insertions = new TreeMap<>();
 249     while (it.hasNext())
 250     {
 251       // check each record for insertions in the CIGAR string
 252       SAMRecord rec = it.next();
 253       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 254               .iterator();
 255       int next = 0;
 256       while (cit.hasNext())
 257       {
 258         CigarElement el = cit.next();
 259         switch (el.getOperator())
 260         {
 261         case I:
 262           // add to insertions list, and move along the read
 263           int refLocation = rec.getReferencePositionAtReadPosition(next);
 264           insertions.put(refLocation, el.getLength());
 265           next += el.getLength();
 266           break;
 267         case M:
 268           // match to reference, move along the read
 269           next += el.getLength();
 270           break;
 271         default:
 272           // deletions, introns etc don't consume any residues from the read
 273           break;
 274         }
 275       }
 276
 277     }
 278     return insertions;
 279   }
 280
 281 }