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   private StringBuilder insertRefSeq(
 129           SortedMap<Integer, Integer> insertions, Range xtent)
 130   {
 131     StringBuilder refseq = new StringBuilder();
 132     int inserted = 0;
 133     for (int p = xtent.start; p < xtent.end; p++)
 134     {
 135       refseq.append('N');
 136     }
 137     for (Map.Entry<Integer, Integer> insert : insertions.entrySet())
 138     {
 139       int inspos = insert.getKey() - xtent.start + inserted;
 140       if (inspos < 0)
 141       {
 142         System.out.println("Ignoring -ve insert position " + insert.getKey()
 143                 + " of " + insert.getValue() +
 144                 " (alpos: " + inspos + ")");
 145         inspos = 0;
 146         // continue;
 147       }
 148       for (int i = 0, j = insert.getValue(); i < j; i++)
 149       {
 150         inserted++;
 151         refseq.insert(inspos, '-');
 152       }
 153     }
 154     return refseq;
 155   }
 156
 157   private void padRef(SequenceI ref, CigarParser cig)
 158   { // pad to column
 159     int padding = cig.firstAlColumn - ref.findIndex(cig.firstRposCol);
 160     System.out.println("Padding " + padding + " to move refseq position "
 161             + cig.firstRposCol + " to " + cig.firstAlColumn + "col.");
 162
 163     ref.insertCharAt(0, padding, '-');
 164   }
 165
 166   @Override
 167   public void parse()
 168   {
 169     // only actually parse if params are set
 170     if (chromosome != null && chromosome != "")
 171     {
 172       SAMRecordIterator it = fileReader.query(chromosome, start, end,
 173               false);
 174       CigarParser parser = new CigarParser('-');
 175       Range[] xtent = new Range[] { new Range(start, end) };
 176       SortedMap<Integer, Integer> insertions[] = parser.getInsertions(it,
 177               xtent);
 178       it.close();
 179
 180       SequenceI refSeq = new Sequence("chr:" + chromosome,
 181               insertRefSeq(insertions[0], xtent[0])
 182                       .toString());
 183
 184       refSeq.setStart(xtent[0].start);
 185       refSeq.setEnd(xtent[0].end);
 186       SequenceI revRefSeq = new Sequence("rev:chr:" + chromosome,
 187               insertRefSeq(insertions[1], xtent[0])
 188                       .toString());
 189       revRefSeq.setStart(xtent[0].start);
 190       revRefSeq.setEnd(xtent[0].end);
 191
 192       // Hack to move the seqs along
 193       padRef(refSeq, parser);
 194       padRef(revRefSeq, parser);
 195
 196       it = fileReader.query(chromosome, start, end, false);
 197
 198       ArrayList<SequenceI> fwd = new ArrayList(), rev = new ArrayList();
 199       while (it.hasNext())
 200       {
 201         SAMRecord rec = it.next();
 202
 203         // set the alignment start to be start of first read (we assume reads
 204         // are sorted)
 205         if (alignmentStart == -1)
 206         {
 207           alignmentStart = rec.getAlignmentStart();
 208         }
 209
 210         // make dataset sequence: start at 1, end at read length
 211         SequenceI seq = new Sequence(
 212                 "" + (rec.getReadNegativeStrandFlag() ? "rev:" : "")
 213                         + rec.getReadName(),
 214                 rec.getReadString().toLowerCase());
 215         seq.setStart(1);
 216         seq.setEnd(rec.getReadLength());
 217         OfInt q = rec.getBaseQualityString().chars()
 218                 .iterator();
 219         int p = seq.getStart();
 220         while (q.hasNext())
 221         {
 222           seq.addSequenceFeature(new SequenceFeature("QUALITY", "", p, p,
 223                   (float) q.next() - ' ', "bamfile"));
 224           p++;
 225         }
 226         String newRead = parser.parseCigarToSequence(rec,
 227                 insertions[rec.getReadNegativeStrandFlag() ? 1 : 0],
 228                 alignmentStart, seq);
 229
 230         // make alignment sequences
 231         SequenceI alsq = seq.deriveSequence();
 232         alsq.setSequence(newRead);
 233
 234         // set start relative to soft clip; assume end is set by Sequence code
 235         alsq.setStart(rec.getStart() - rec.getUnclippedStart() + 1);
 236         if (rec.getReadNegativeStrandFlag())
 237         {
 238           rev.add(alsq);
 239         }
 240         else
 241         {
 242           fwd.add(alsq);
 243         }
 244       }
 245       // Add forward
 246       if (fwd.size() > 0)
 247       {
 248         seqs.add(refSeq); // FIXME needs to be optional, and properly padded
 249         seqs.addAll(fwd);
 250         fwd.clear();
 251       }
 252       // and reverse strand reads.
 253       if (rev.size() > 0)
 254       {
 255         seqs.add(revRefSeq); // FIXME needs to be optional and properly padded
 256         seqs.addAll(rev);
 257         rev.clear();
 258       }
 259     }
 260   }
 261
 262   /**
 263    * Get the list of chromosomes or contigs from the file (listed in SQ entries
 264    * in BAM file header)
 265    *
 266    * @return array of chromosome/contig strings
 267    */
 268   @Override
 269   public Object[] preprocess()
 270   {
 271     List<SAMSequenceRecord> refSeqs = fileReader.getFileHeader()
 272             .getSequenceDictionary().getSequences();
 273     List<String> chrs = new ArrayList<>();
 274
 275     for (SAMSequenceRecord ref : refSeqs)
 276     {
 277       chrs.add(ref.getSequenceName());
 278     }
 279     return chrs.toArray();
 280   }
 281
 282   public void setOptions(String chr, int s, int e)
 283   {
 284     chromosome = chr;
 285     start = s;
 286     end = e;
 287     suffix = chromosome + ":" + start + "-" + end;
 288   }
 289
 290   public boolean parseSuffix()
 291   {
 292     if (suffix == null)
 293     {
 294       return false;
 295     }
 296     int csep = suffix.indexOf(":");
 297     int rsep = suffix.indexOf("-", csep);
 298     if (csep < 0 || rsep < 0 || suffix.length() - rsep <= 1)
 299     {
 300       return false;
 301     }
 302     String chr, p1, p2;
 303     chr = suffix.substring(0, csep);
 304     p1 = suffix.substring(csep + 1, rsep);
 305     p2 = suffix.substring(rsep + 1);
 306     int cstart, cend;
 307     try
 308     {
 309       cstart = Integer.parseInt(p1);
 310       cend = Integer.parseInt(p2);
 311     } catch (Exception e)
 312     {
 313       warningMessage = (warningMessage == null ? "" : warningMessage + "\n")
 314               + "Couldn't parse range from " + suffix;
 315       return false;
 316     }
 317     chromosome = chr;
 318     start = cstart;
 319     end = cend;
 320     return true;
 321   }
 322 }