*/
package jalview.io;
+import jalview.datamodel.CigarParser;
+import jalview.datamodel.Sequence;
import jalview.datamodel.SequenceI;
import java.io.File;
import java.io.IOException;
-import java.util.Arrays;
-import java.util.Iterator;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.SortedMap;
-import htsjdk.samtools.CigarElement;
import htsjdk.samtools.SAMRecord;
import htsjdk.samtools.SAMRecordIterator;
+import htsjdk.samtools.SAMSequenceRecord;
import htsjdk.samtools.SamReader;
import htsjdk.samtools.SamReaderFactory;
import htsjdk.samtools.ValidationStringency;
public class BamFile extends AlignFile
{
+ // SAM/BAM file reader
+ private SamReader fileReader;
- SamReader fileReader;
+ // start position to read from
+ private int start = -1;
+
+ // end position to read to
+ private int end = -1;
+
+ // chromosome/contig to read
+ private String chromosome = "";
+
+ // first position in alignment
+ private int alignmentStart = -1;
/**
* Creates a new BamFile object.
* Creates a new BamFile object.
*
* @param inFile
- * DOCUMENT ME!
+ * Name of file to read
* @param sourceType
- * DOCUMENT ME!
+ * Whether data source is file, url or other type of data source
*
* @throws IOException
- * DOCUMENT ME!
*/
public BamFile(String inFile, DataSourceType sourceType)
throws IOException
{
- super(inFile, sourceType);
+ super(true, inFile, sourceType);
final SamReaderFactory factory = SamReaderFactory.makeDefault()
.enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
fileReader = factory.open(new File(inFile));
}
+ /**
+ * Creates a new BamFile object
+ *
+ * @param source
+ * wrapper for datasource
+ * @throws IOException
+ */
public BamFile(FileParse source) throws IOException
{
+ super(true, source);
final SamReaderFactory factory = SamReaderFactory.makeDefault()
.enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS,
SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
// File-based bam
fileReader = factory.open(source.inFile);
- parse();
}
@Override
}
@Override
- public void parse() throws IOException
+ public void parse()
{
- SequenceI seq = null;
- SAMRecordIterator it = fileReader.iterator();
- while (it.hasNext())
+ // only actually parse if params are set
+ if (chromosome != null && chromosome != "")
{
- SAMRecord rec = it.next();
- String read = rec.getReadString();
- int start = rec.getStart();
- int end = rec.getEnd();
-
- Iterator<CigarElement> cit = rec.getCigar().getCigarElements()
- .iterator();
-
- seq = parseId(rec.getReadName());
- String cigarredRead = parseCigarToSequence(read, cit, '-');
- seq.setSequence(cigarredRead);
- seq.setStart(start);
- seq.setEnd(end);
- seqs.addElement(seq);
+ SAMRecordIterator it = fileReader.query(chromosome, start, end,
+ false);
+ CigarParser parser = new CigarParser('-');
+ SortedMap<Integer, Integer> insertions = parser.getInsertions(it);
+ it.close();
+
+ it = fileReader.query(chromosome, start, end, false);
+ while (it.hasNext())
+ {
+ SAMRecord rec = it.next();
+
+ // set the alignment start to be start of first read (we assume reads
+ // are sorted)
+ if (alignmentStart == -1)
+ {
+ alignmentStart = rec.getAlignmentStart();
+ }
+
+ // make dataset sequence: start at 1, end at read length
+ SequenceI seq = new Sequence(rec.getReadName(),
+ rec.getReadString().toLowerCase());
+ seq.setStart(1);
+ seq.setEnd(rec.getReadLength());
+
+ String newRead = parser.parseCigarToSequence(rec, insertions,
+ alignmentStart);
+
+ // make alignment sequences
+ SequenceI alsq = seq.deriveSequence();
+ alsq.setSequence(newRead);
+
+ // set start relative to soft clip; assume end is set by Sequence code
+ alsq.setStart(rec.getStart() - rec.getUnclippedStart() + 1);
+ seqs.add(alsq);
+ }
}
-
}
/**
- * Apply the CIGAR string to a read sequence and return the updated read
+ * Get the list of chromosomes or contigs from the file (listed in SQ entries
+ * in BAM file header)
*
- * @param read
- * the read to update
- * @param it
- * iterator over cigar elements
- * @param gapChar
- * gap character to use
- * @return string representing read with gaps, clipping etc applied
+ * @return array of chromosome/contig strings
*/
- private String parseCigarToSequence(String read,
- Iterator<CigarElement> it, char gapChar)
+ @Override
+ public Object[] preprocess()
{
- StringBuilder newRead = new StringBuilder();
- int next = 0;
+ List<SAMSequenceRecord> refSeqs = fileReader.getFileHeader()
+ .getSequenceDictionary().getSequences();
+ List<String> chrs = new ArrayList<>();
- while (it.hasNext())
+ for (SAMSequenceRecord ref : refSeqs)
{
- CigarElement el = it.next();
- int length = el.getLength();
- switch (el.getOperator())
- {
- case M:
- // matched residues
- newRead.append(read.substring(next, next + length - 1));
- next += length;
- break;
- case N: // intron in RNA
- case D: // deletion
- // add gaps
- char[] gaps = new char[length];
- Arrays.fill(gaps, gapChar);
- newRead.append(gaps);
- break;
- case S:
- // soft clipping - just skip this bit of the read
- // do nothing
- next += length;
- break;
- case I:
- // add gaps to the reference sequence, which we know nothing about just
- // now
- // TODO
- newRead.append(read.substring(next, next + length - 1));
- break;
- case H:
- // hard clipping - this stretch will not appear in the read
- break;
- default:
- // P, X EQ don't know what to do with these
- break;
- }
-
+ chrs.add(ref.getSequenceName());
}
- return newRead.toString();
+ return chrs.toArray();
}
+ public void setOptions(String chr, int s, int e)
+ {
+ chromosome = chr;
+ start = s;
+ end = e;
+ }
}
git://source.jalview.org / jalview.git / blobdiff