X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fio%2FBamFile.java;h=10e2d0ec7db91f30360e5d5309d81b05e075c896;hb=846629a1a7ab7715d14a18296bc3024c5de4ac5c;hp=e9919ff11d0098b5c6160953a5e14d57f6b2a642;hpb=932cc89cff4bb8fb8d5885d1c4bdc110f34486d9;p=jalview.git diff --git a/src/jalview/io/BamFile.java b/src/jalview/io/BamFile.java index e9919ff..10e2d0e 100644 --- a/src/jalview/io/BamFile.java +++ b/src/jalview/io/BamFile.java @@ -20,27 +20,39 @@ */ 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 java.util.TreeMap; -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. @@ -53,17 +65,16 @@ public class BamFile extends AlignFile * 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) @@ -71,8 +82,16 @@ public class BamFile extends AlignFile 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) @@ -80,7 +99,6 @@ public class BamFile extends AlignFile // File-based bam fileReader = factory.open(source.inFile); - parse(); } @Override @@ -91,191 +109,73 @@ public class BamFile extends AlignFile } @Override - public void parse() throws IOException + public void parse() { - SAMRecordIterator it = fileReader.iterator(); - SortedMap insertions = getInsertions(it); - it.close(); - - 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(); + SAMRecordIterator it = fileReader.query(chromosome, start, end, + false); + CigarParser parser = new CigarParser('-'); + SortedMap insertions = parser.getInsertions(it); + it.close(); + + it = fileReader.query(chromosome, start, end, false); + while (it.hasNext()) + { + SAMRecord rec = it.next(); - SequenceI seq = new Sequence(rec.getReadName(), rec.getReadString()); + // set the alignment start to be start of first read (we assume reads + // are sorted) + if (alignmentStart == -1) + { + alignmentStart = rec.getAlignmentStart(); + } - String cigarredRead = parseCigarToSequence(read, rec, '-'); + // 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()); - SequenceI alsq = seq.deriveSequence(); - alsq.setSequence(cigarredRead); - alsq.setStart(start); - alsq.setEnd(end); - seqs.add(alsq); - } + String newRead = parser.parseCigarToSequence(rec, insertions, + alignmentStart, seq); - for (SequenceI seq : seqs) - { - int insertCount = 0; - SortedMap seqInserts = insertions.subMap(0, - seq.getEnd()); // TODO start point should be start of alignment - // not 0 - seqInserts.forEach(action); - while ((nextInsertion != -1) && (nextInsertion < seq.getEnd())) - { - seq.insertCharsAt(nextInsertion + insertCount, '-'); - // nextInsertion = insertions.nextSetBit(nextInsertion + 1); - insertCount++; + // 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); } } } - /* - 1234567 - ABCDEFG insert 3,4 - - R1: insert 3 = gap before 3rd res - AB.CDEFG - - R2: insert 4 = gap before 4th res - ABC.DEFG - - R3: insert 3,4 = 2 gaps before 3rd res - AB..CDEFG - - => AB--C-DEFG - - So need to store insertions as (pos, length) giving here (3,1),(4,1),(3,2) - then can sub longest length at position so (3,2), (4,1) - Then when inserting need to include previously inserted count AND determine if sequence already has an insertion or needs gap(s) added - - */ /** - * 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 rec - * the corresponding SAM record - * @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, - SAMRecord rec, char gapChar) + @Override + public Object[] preprocess() { - Iterator it = rec.getCigar().getCigarElements() - .iterator(); - - StringBuilder newRead = new StringBuilder(); - int next = 0; - int ap = 0; // position of query region qi.start; + List refSeqs = fileReader.getFileHeader() + .getSequenceDictionary().getSequences(); + List chrs = new ArrayList<>(); - // pad with spaces before read - for (int alp = rec - .getReferencePositionAtReadPosition(next + 1); ap < alp;) + for (SAMSequenceRecord ref : refSeqs) { - newRead.append(" "); - ap++; + chrs.add(ref.getSequenceName()); } - - while (it.hasNext()) - { - CigarElement el = it.next(); - int length = el.getLength(); - char[] gaps; - - switch (el.getOperator()) - { - case M: - // matched residues - newRead.append( - read.substring(next, next + length)); - next += length; - break; - case N: // intron in RNA - case D: // deletion - // add gaps - 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 - - // work out how many gaps we need before the start of the soft clip - - // don't do this at the end of the read! - if (next == 0) // at start of read - { - int numgaps = rec.getUnclippedStart(); - gaps = new char[numgaps]; - Arrays.fill(gaps, ' '); - newRead.append(gaps); - } - - newRead.append( - read.substring(next, next + length).toLowerCase()); - next += length; - break; - case I: - // the reference sequence and other reads should have been gapped for - // this insertion, so just add in the residues - newRead.append(read.substring(next, next + length)); - next += length; - 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; - } - - } - return newRead.toString(); + return chrs.toArray(); } - /** - * Get list of positions inserted to the reference sequence - * - * @param it - * @return - */ - private SortedMap getInsertions(Iterator it) + public void setOptions(String chr, int s, int e) { - SortedMap insertions = new TreeMap<>(); - while (it.hasNext()) - { - // check each record for insertions in the CIGAR string - SAMRecord rec = it.next(); - Iterator cit = rec.getCigar().getCigarElements() - .iterator(); - int next = 0; - while (cit.hasNext()) - { - CigarElement el = cit.next(); - switch (el.getOperator()) - { - case I: - // add to insertions list, and move along the read - int refLocation = rec.getReferencePositionAtReadPosition(next); - insertions.put(refLocation, el.getLength()); - next += el.getLength(); - break; - case M: - // match to reference, move along the read - next += el.getLength(); - break; - default: - // deletions, introns etc don't consume any residues from the read - break; - } - } - - } - return insertions; + chromosome = chr; + start = s; + end = e; } - }