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.SequenceI;
  24
  25 import java.io.File;
  26 import java.io.IOException;
  27 import java.util.Arrays;
  28 import java.util.Iterator;
  29
  30 import htsjdk.samtools.CigarElement;
  31 import htsjdk.samtools.SAMRecord;
  32 import htsjdk.samtools.SAMRecordIterator;
  33 import htsjdk.samtools.SamReader;
  34 import htsjdk.samtools.SamReaderFactory;
  35 import htsjdk.samtools.ValidationStringency;
  36
  37 public class BamFile extends AlignFile
  38 {
  39
  40   SamReader fileReader;
  41
  42   /**
  43    * Creates a new BamFile object.
  44    */
  45   public BamFile()
  46   {
  47   }
  48
  49   /**
  50    * Creates a new BamFile object.
  51    *
  52    * @param inFile
  53    *          DOCUMENT ME!
  54    * @param sourceType
  55    *          DOCUMENT ME!
  56    *
  57    * @throws IOException
  58    *           DOCUMENT ME!
  59    */
  60   public BamFile(String inFile, DataSourceType sourceType)
  61           throws IOException
  62   {
  63     super(inFile, sourceType);
  64     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  65             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  66                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  67             .validationStringency(ValidationStringency.SILENT);
  68     fileReader = factory.open(new File(inFile));
  69   }
  70
  71   public BamFile(FileParse source) throws IOException
  72   {
  73     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  74             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  75                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  76             .validationStringency(ValidationStringency.SILENT);
  77
  78     // File-based bam
  79     fileReader = factory.open(source.inFile);
  80     parse();
  81   }
  82
  83   @Override
  84   public String print(SequenceI[] seqs, boolean jvsuffix)
  85   {
  86     // TODO Auto-generated method stub
  87     return null;
  88   }
  89
  90   @Override
  91   public void parse() throws IOException
  92   {
  93     SequenceI seq = null;
  94     SAMRecordIterator it = fileReader.iterator();
  95     while (it.hasNext())
  96     {
  97       SAMRecord rec = it.next();
  98       String read = rec.getReadString();
  99       int start = rec.getStart();
 100       int end = rec.getEnd();
 101
 102       Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
 103               .iterator();
 104
 105       seq = parseId(rec.getReadName());
 106       String cigarredRead = parseCigarToSequence(read, cit, '-');
 107       seq.setSequence(cigarredRead);
 108       seq.setStart(start);
 109       seq.setEnd(end);
 110       seqs.addElement(seq);
 111     }
 112
 113   }
 114
 115   /**
 116    * Apply the CIGAR string to a read sequence and return the updated read
 117    *
 118    * @param read
 119    *          the read to update
 120    * @param it
 121    *          iterator over cigar elements
 122    * @param gapChar
 123    *          gap character to use
 124    * @return string representing read with gaps, clipping etc applied
 125    */
 126   private String parseCigarToSequence(String read,
 127           Iterator<CigarElement> it, char gapChar)
 128   {
 129     StringBuilder newRead = new StringBuilder();
 130     int next = 0;
 131
 132     while (it.hasNext())
 133     {
 134       CigarElement el = it.next();
 135       int length = el.getLength();
 136       switch (el.getOperator())
 137       {
 138       case M:
 139         // matched residues
 140         newRead.append(read.substring(next, next + length - 1));
 141         next += length;
 142         break;
 143       case N: // intron in RNA
 144       case D: // deletion
 145         // add gaps
 146         char[] gaps = new char[length];
 147         Arrays.fill(gaps, gapChar);
 148         newRead.append(gaps);
 149         break;
 150       case S:
 151         // soft clipping - just skip this bit of the read
 152         // do nothing
 153         next += length;
 154         break;
 155       case I:
 156         // add gaps to the reference sequence, which we know nothing about just
 157         // now
 158         // TODO
 159         newRead.append(read.substring(next, next + length - 1));
 160         break;
 161       case H:
 162         // hard clipping - this stretch will not appear in the read
 163         break;
 164       default:
 165         // P, X EQ don't know what to do with these
 166         break;
 167       }
 168
 169     }
 170     return newRead.toString();
 171   }
 172
 173 }