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