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.bin.Cache;
  24 import jalview.datamodel.CigarParser;
  25 import jalview.datamodel.Range;
  26 import jalview.datamodel.Sequence;
  27 import jalview.datamodel.SequenceFeature;
  28 import jalview.datamodel.SequenceI;
  29
  30 import java.io.File;
  31 import java.io.IOException;
  32 import java.net.URL;
  33 import java.util.ArrayList;
  34 import java.util.List;
  35 import java.util.Map;
  36 import java.util.PrimitiveIterator.OfInt;
  37 import java.util.SortedMap;
  38
  39 import htsjdk.samtools.SAMFormatException;
  40 import htsjdk.samtools.SAMRecord;
  41 import htsjdk.samtools.SAMRecordIterator;
  42 import htsjdk.samtools.SAMSequenceRecord;
  43 import htsjdk.samtools.SamInputResource;
  44 import htsjdk.samtools.SamReader;
  45 import htsjdk.samtools.SamReaderFactory;
  46 import htsjdk.samtools.ValidationStringency;
  47
  48 public class BamFile extends AlignFile
  49 {
  50   // SAM/BAM file reader
  51   private SamReader fileReader;
  52
  53   // start position to read from
  54   private int start = -1;
  55
  56   // end position to read to
  57   private int end = -1;
  58
  59   // chromosome/contig to read
  60   private String chromosome = "";
  61
  62   // first position in alignment
  63   private int alignmentStart = -1;
  64
  65   private File _bamFile;
  66
  67   /**
  68    * Creates a new BamFile object.
  69    */
  70   public BamFile()
  71   {
  72   }
  73
  74   /**
  75    * Creates a new BamFile object.
  76    *
  77    * @param inFile
  78    *          Name of file to read
  79    * @param sourceType
  80    *          Whether data source is file, url or other type of data source
  81    *
  82    * @throws IOException
  83    */
  84   public BamFile(String inFile, DataSourceType sourceType)
  85           throws IOException
  86   {
  87     super(true, inFile, sourceType);
  88     _bamFile = new File(inFile);
  89     initFileReader();
  90   }
  91   private void initFileReader() throws IOException
  92   {
  93     final SamReaderFactory factory = SamReaderFactory.makeDefault()
  94             .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
  95                     SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
  96             .validationStringency(ValidationStringency.SILENT);
  97     // File-based bam
  98     if (_bamFile!=null)
  99     {
 100       fileReader = factory.open(_bamFile); // will need to be adapted for JalviewJS/etc
 101     }
 102     else
 103     {
 104       // try and locate index ?
 105       String index = getDataName() + ".bai";
 106       fileReader = factory.open(SamInputResource.of(getDataName())
 107               .index(new URL(index)));
 108     }
 109   }
 110   /**
 111    * Creates a new BamFile object
 112    *
 113    * @param source
 114    *          wrapper for datasource
 115    * @throws IOException
 116    */
 117   public BamFile(FileParse source) throws IOException
 118   {
 119     super(false, source);
 120     parseSuffix();
 121     if (source.getDataSourceType() == DataSourceType.FILE)
 122      {
 123        _bamFile = source.inFile;
 124      } else {
 125
 126      }
 127     initFileReader();
 128     doParse();
 129   }
 130
 131   @Override
 132   public String print(SequenceI[] seqs, boolean jvsuffix)
 133   {
 134     // TODO Auto-generated method stub
 135     return null;
 136   }
 137
 138   private StringBuilder insertRefSeq(
 139           SortedMap<Integer, Integer> insertions, Range xtent)
 140   {
 141     StringBuilder refseq = new StringBuilder();
 142     int inserted = 0;
 143     for (int p = xtent.start; p < xtent.end; p++)
 144     {
 145       refseq.append('N');
 146     }
 147     for (Map.Entry<Integer, Integer> insert : insertions.entrySet())
 148     {
 149       int inspos = insert.getKey() - xtent.start + inserted;
 150       if (inspos < 0)
 151       {
 152         System.out.println("Ignoring -ve insert position " + insert.getKey()
 153                 + " of " + insert.getValue() +
 154                 " (alpos: " + inspos + ")");
 155         inspos = 0;
 156         // continue;
 157       }
 158       for (int i = 0, j = insert.getValue(); i < j; i++)
 159       {
 160         inserted++;
 161         refseq.insert(inspos, '-');
 162       }
 163     }
 164     return refseq;
 165   }
 166
 167   private void padRef(SequenceI ref, CigarParser cig)
 168   { // pad to column
 169     int padding = cig.firstAlColumn - ref.findIndex(cig.firstRposCol);
 170     System.out.println("Padding " + padding + " to move refseq position "
 171             + cig.firstRposCol + " to " + cig.firstAlColumn + "col.");
 172
 173     ref.insertCharAt(0, padding, '-');
 174   }
 175
 176   @Override
 177   public void parse()
 178   {
 179     boolean needToReopen=false;
 180     // only actually parse if params are set
 181     if (chromosome != null && chromosome != "")
 182     {
 183       SAMRecordIterator it;
 184       try {
 185         it = fileReader.query(chromosome, start, end,
 186               false);
 187       } catch (UnsupportedOperationException ex)
 188       {
 189         needToReopen=true;
 190         // could be a sam text file, so we just iterate through without query
 191         it = fileReader.iterator();
 192       }
 193       CigarParser parser = new CigarParser('-');
 194       Range[] xtent = new Range[] { new Range(start, end) };
 195       SortedMap<Integer, Integer> insertions[] = parser.getInsertions(it,
 196               xtent);
 197       it.close();
 198
 199       SequenceI refSeq = new Sequence("chr:" + chromosome,
 200               insertRefSeq(insertions[0], xtent[0])
 201                       .toString());
 202
 203       refSeq.setStart(xtent[0].start);
 204       refSeq.setEnd(xtent[0].end);
 205       SequenceI revRefSeq = new Sequence("rev:chr:" + chromosome,
 206               insertRefSeq(insertions[1], xtent[0])
 207                       .toString());
 208       revRefSeq.setStart(xtent[0].start);
 209       revRefSeq.setEnd(xtent[0].end);
 210
 211       // Hack to move the seqs along
 212       padRef(refSeq, parser);
 213       padRef(revRefSeq, parser);
 214
 215       if (needToReopen)
 216       {
 217         try {
 218           initFileReader();
 219         } catch (IOException x)
 220         {
 221           Cache.log.warn("Couldn't reopen S/BAM file",x);
 222         }
 223       }
 224       try {
 225         it = fileReader.query(chromosome, start, end,
 226               false);
 227       } catch (UnsupportedOperationException ex)
 228       {
 229         // could be a sam text file, so we just iterate through without query
 230         it = fileReader.iterator();
 231       }
 232
 233       ArrayList<SequenceI> fwd = new ArrayList(), rev = new ArrayList();
 234       while (it.hasNext())
 235       {
 236         SAMRecord rec=null;
 237         try {
 238           rec = it.next();
 239         } catch (SAMFormatException q)
 240         {
 241           Cache.log.info("Bailing on bad SAM line again",q);
 242           break;
 243         }
 244
 245         // set the alignment start to be start of first read (we assume reads
 246         // are sorted)
 247         if (alignmentStart == -1)
 248         {
 249           alignmentStart = rec.getAlignmentStart();
 250         }
 251
 252         // make dataset sequence: start at 1, end at read length
 253         SequenceI seq = new Sequence(
 254                 "" + (rec.getReadNegativeStrandFlag() ? "rev:" : "")
 255                         + rec.getReadName(),
 256                 rec.getReadString().toLowerCase());
 257         seq.setStart(1);
 258         seq.setEnd(rec.getReadLength());
 259         OfInt q = rec.getBaseQualityString().chars()
 260                 .iterator();
 261         int p = seq.getStart();
 262         while (q.hasNext())
 263         {
 264           seq.addSequenceFeature(new SequenceFeature("QUALITY", "", p, p,
 265                   (float) q.next() - ' ', "bamfile"));
 266           p++;
 267         }
 268         String newRead = parser.parseCigarToSequence(rec,
 269                 insertions[rec.getReadNegativeStrandFlag() ? 1 : 0],
 270                 alignmentStart, seq);
 271
 272         // make alignment sequences
 273         SequenceI alsq = seq.deriveSequence();
 274         alsq.setSequence(newRead);
 275
 276         // set start relative to soft clip; assume end is set by Sequence code
 277         alsq.setStart(rec.getStart() - rec.getUnclippedStart() + 1);
 278         if (rec.getReadNegativeStrandFlag())
 279         {
 280           rev.add(alsq);
 281         }
 282         else
 283         {
 284           fwd.add(alsq);
 285         }
 286       }
 287       // Add forward
 288       if (fwd.size() > 0)
 289       {
 290         seqs.add(refSeq); // FIXME needs to be optional, and properly padded
 291         seqs.addAll(fwd);
 292         fwd.clear();
 293       }
 294       // and reverse strand reads.
 295       if (rev.size() > 0)
 296       {
 297         seqs.add(revRefSeq); // FIXME needs to be optional and properly padded
 298         seqs.addAll(rev);
 299         rev.clear();
 300       }
 301     }
 302   }
 303
 304   /**
 305    * Get the list of chromosomes or contigs from the file (listed in SQ entries
 306    * in BAM file header)
 307    *
 308    * @return array of chromosome/contig strings
 309    */
 310   @Override
 311   public Object[] preprocess()
 312   {
 313     List<SAMSequenceRecord> refSeqs = fileReader.getFileHeader()
 314             .getSequenceDictionary().getSequences();
 315     List<String> chrs = new ArrayList<>();
 316
 317     for (SAMSequenceRecord ref : refSeqs)
 318     {
 319       chrs.add(ref.getSequenceName());
 320     }
 321     return chrs.toArray();
 322   }
 323
 324   public void setOptions(String chr, int s, int e)
 325   {
 326     chromosome = chr;
 327     start = s;
 328     end = e;
 329     suffix = chromosome + ":" + start + "-" + end;
 330   }
 331
 332   public boolean parseSuffix()
 333   {
 334     if (suffix == null)
 335     {
 336       return false;
 337     }
 338     int csep = suffix.indexOf(":");
 339     int rsep = suffix.indexOf("-", csep);
 340     if (csep < 0 || rsep < 0 || suffix.length() - rsep <= 1)
 341     {
 342       return false;
 343     }
 344     String chr, p1, p2;
 345     chr = suffix.substring(0, csep);
 346     p1 = suffix.substring(csep + 1, rsep);
 347     p2 = suffix.substring(rsep + 1);
 348     int cstart, cend;
 349     try
 350     {
 351       cstart = Integer.parseInt(p1);
 352       cend = Integer.parseInt(p2);
 353     } catch (Exception e)
 354     {
 355       warningMessage = (warningMessage == null ? "" : warningMessage + "\n")
 356               + "Couldn't parse range from " + suffix;
 357       return false;
 358     }
 359     chromosome = chr;
 360     start = cstart;
 361     end = cend;
 362     return true;
 363   }
 364 }