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.CigarParser;
  24 import jalview.datamodel.Range;
  25 import jalview.datamodel.Sequence;
  26 import jalview.datamodel.SequenceFeature;
  27 import jalview.datamodel.SequenceI;
  28
  29 import java.io.File;
  30 import java.io.IOException;
  31 import java.net.URL;
  32 import java.util.ArrayList;
  33 import java.util.List;
  34 import java.util.Map;
  35 import java.util.PrimitiveIterator.OfInt;
  36 import java.util.SortedMap;
  37
  38 import htsjdk.samtools.SAMRecord;
  39 import htsjdk.samtools.SAMRecordIterator;
  40 import htsjdk.samtools.SAMSequenceRecord;
  41 import htsjdk.samtools.SamInputResource;
  42 import htsjdk.samtools.SamReader;
  43 import htsjdk.samtools.SamReaderFactory;
  44 import htsjdk.samtools.ValidationStringency;
  45
  46 public class BamFile extends AlignFile
  47 {
  48   // SAM/BAM file reader
  49   private SamReader fileReader;
  50
  51   // start position to read from
  52   private int start = -1;
  53
  54   // end position to read to
  55   private int end = -1;
  56
  57   // chromosome/contig to read
  58   private String chromosome = "";
  59
  60   // first position in alignment
  61   private int alignmentStart = -1;
  62
  63   /**
  64    * Creates a new BamFile object.
  65    */
  66   public BamFile()
  67   {
  68   }
  69
  70   /**
  71    * Creates a new BamFile object.
  72    *
  73    * @param inFile
  74    *          Name of file to read
  75    * @param sourceType
  76    *          Whether data source is file, url or other type of data source
  77    *
  78    * @throws IOException
  79    */
  80   public BamFile(String inFile, DataSourceType sourceType)
  81           throws IOException
  82   {
  83     super(true, inFile, sourceType);
  84     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  85             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  86                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  87             .validationStringency(ValidationStringency.SILENT);
  88     fileReader = factory.open(new File(inFile));
  89   }
  90
  91   /**
  92    * Creates a new BamFile object
  93    *
  94    * @param source
  95    *          wrapper for datasource
  96    * @throws IOException
  97    */
  98   public BamFile(FileParse source) throws IOException
  99   {
 100     super(true, source);
 101     parseSuffix();
 102     final SamReaderFactory factory = SamReaderFactory.makeDefault()
 103             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
 104                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
 105             .validationStringency(ValidationStringency.SILENT);
 106
 107     // File-based bam
 108     if (source.getDataSourceType() == DataSourceType.FILE)
 109     {
 110       fileReader = factory.open(source.inFile);
 111     }
 112     else
 113     {
 114       // locate index ?
 115       String index = source.getDataName() + ".bai";
 116       fileReader = factory.open(SamInputResource.of(source.getDataName())
 117               .index(new URL(index)));
 118     }
 119   }
 120
 121   @Override
 122   public String print(SequenceI[] seqs, boolean jvsuffix)
 123   {
 124     // TODO Auto-generated method stub
 125     return null;
 126   }
 127
 128   @Override
 129   public void parse()
 130   {
 131     // only actually parse if params are set
 132     if (chromosome != null && chromosome != "")
 133     {
 134       SAMRecordIterator it = fileReader.query(chromosome, start, end,
 135               false);
 136       CigarParser parser = new CigarParser('-');
 137       SortedMap<Integer, Integer> insertions = parser.getInsertions(it);
 138       it.close();
 139
 140       it = fileReader.query(chromosome, start, end, false);
 141       while (it.hasNext())
 142       {
 143         SAMRecord rec = it.next();
 144
 145         // set the alignment start to be start of first read (we assume reads
 146         // are sorted)
 147         if (alignmentStart == -1)
 148         {
 149           alignmentStart = rec.getAlignmentStart();
 150         }
 151
 152         // make dataset sequence: start at 1, end at read length
 153         SequenceI seq = new Sequence(rec.getReadName(),
 154                 rec.getReadString().toLowerCase());
 155         seq.setStart(1);
 156         seq.setEnd(rec.getReadLength());
 157         OfInt q = rec.getBaseQualityString().chars()
 158                 .iterator();
 159         int p = seq.getStart();
 160         while (q.hasNext())
 161         {
 162           seq.addSequenceFeature(new SequenceFeature("QUALITY", "", p, p,
 163                   (float) q.next() - ' ', "bamfile"));
 164           p++;
 165         }
 166         String newRead = parser.parseCigarToSequence(rec, insertions,
 167                 alignmentStart, seq);
 168
 169         // make alignment sequences
 170         SequenceI alsq = seq.deriveSequence();
 171         alsq.setSequence(newRead);
 172
 173         // set start relative to soft clip; assume end is set by Sequence code
 174         alsq.setStart(rec.getStart() - rec.getUnclippedStart() + 1);
 175         seqs.add(alsq);
 176       }
 177     }
 178   }
 179
 180   /**
 181    * Get the list of chromosomes or contigs from the file (listed in SQ entries
 182    * in BAM file header)
 183    *
 184    * @return array of chromosome/contig strings
 185    */
 186   @Override
 187   public Object[] preprocess()
 188   {
 189     List<SAMSequenceRecord> refSeqs = fileReader.getFileHeader()
 190             .getSequenceDictionary().getSequences();
 191     List<String> chrs = new ArrayList<>();
 192
 193     for (SAMSequenceRecord ref : refSeqs)
 194     {
 195       chrs.add(ref.getSequenceName());
 196     }
 197     return chrs.toArray();
 198   }
 199
 200   public void setOptions(String chr, int s, int e)
 201   {
 202     chromosome = chr;
 203     start = s;
 204     end = e;
 205     suffix = chromosome + ":" + start + "-" + end;
 206   }
 207
 208   public boolean parseSuffix()
 209   {
 210     if (suffix == null)
 211     {
 212       return false;
 213     }
 214     int csep = suffix.indexOf(":");
 215     int rsep = suffix.indexOf("-", csep);
 216     if (csep < 0 || rsep < 0 || suffix.length() - rsep <= 1)
 217     {
 218       return false;
 219     }
 220     String chr, p1, p2;
 221     chr = suffix.substring(0, csep);
 222     p1 = suffix.substring(csep + 1, rsep);
 223     p2 = suffix.substring(rsep + 1);
 224     int cstart, cend;
 225     try
 226     {
 227       cstart = Integer.parseInt(p1);
 228       cend = Integer.parseInt(p2);
 229     } catch (Exception e)
 230     {
 231       warningMessage = (warningMessage == null ? "" : warningMessage + "\n")
 232               + "Couldn't parse range from " + suffix;
 233       return false;
 234     }
 235     chromosome = chr;
 236     start = cstart;
 237     end = cend;
 238     return true;
 239   }
 240 }