Examples with SamReader htsjdk.samtools.SamReader used on opensource projects

Example 41 with SamReader

use of htsjdk.samtools.SamReader in project ASCIIGenome by dariober.

the class Utils method countReadsInWindow.

/**
 * Count reads in interval using the given filters.
 * @param bam
 * @param gc
 * @param filters List of filters to apply
 * @return
 * @throws MalformedURLException
 */
public static long countReadsInWindow(String bam, GenomicCoords gc, List<SamRecordFilter> filters) throws MalformedURLException {
    /*  ------------------------------------------------------ */
    /* This chunk prepares SamReader from local bam or URL bam */
    UrlValidator urlValidator = new UrlValidator();
    SamReaderFactory srf = SamReaderFactory.make();
    srf.validationStringency(ValidationStringency.SILENT);
    SamReader samReader;
    if (urlValidator.isValid(bam)) {
        samReader = srf.open(SamInputResource.of(new URL(bam)).index(new URL(bam + ".bai")));
    } else {
        samReader = srf.open(new File(bam));
    }
    /*  ------------------------------------------------------ */
    long cnt = 0;
    Iterator<SAMRecord> sam = samReader.query(gc.getChrom(), gc.getFrom(), gc.getTo(), false);
    AggregateFilter aggregateFilter = new AggregateFilter(filters);
    while (sam.hasNext()) {
        SAMRecord rec = sam.next();
        if (!aggregateFilter.filterOut(rec)) {
            cnt++;
        }
    }
    try {
        samReader.close();
    } catch (IOException e) {
        e.printStackTrace();
    }
    return cnt;
}

Example 42 with SamReader

use of htsjdk.samtools.SamReader in project ASCIIGenome by dariober.

the class Utils method sortAndIndexSamOrBam.

/**
 * Sort and index input sam or bam.
 * @throws IOException
 */
public static void sortAndIndexSamOrBam(String inSamOrBam, String sortedBam, boolean deleteOnExit) throws IOException {
    /*  ------------------------------------------------------ */
    /* This chunk prepares SamReader from local bam or URL bam */
    UrlValidator urlValidator = new UrlValidator();
    SamReaderFactory srf = SamReaderFactory.make();
    srf.validationStringency(ValidationStringency.SILENT);
    SamReader samReader;
    if (urlValidator.isValid(inSamOrBam)) {
        samReader = SamReaderFactory.makeDefault().open(SamInputResource.of(new URL(inSamOrBam)));
    } else {
        samReader = srf.open(new File(inSamOrBam));
    }
    /*  ------------------------------------------------------ */
    samReader.getFileHeader().setSortOrder(SortOrder.coordinate);
    File out = new File(sortedBam);
    if (deleteOnExit) {
        out.deleteOnExit();
        File idx = new File(out.getAbsolutePath().replaceAll("\\.bam$", "") + ".bai");
        idx.deleteOnExit();
    }
    SAMFileWriter outputSam = new SAMFileWriterFactory().setCreateIndex(true).makeSAMOrBAMWriter(samReader.getFileHeader(), false, out);
    for (final SAMRecord samRecord : samReader) {
        outputSam.addAlignment(samRecord);
    }
    samReader.close();
    outputSam.close();
}

Example 43 with SamReader

use of htsjdk.samtools.SamReader in project ASCIIGenome by dariober.

the class Utils method getSamReader.

/**
 * Prepare SamReader from local bam or URL bam
 */
public static SamReader getSamReader(String workFilename) throws MalformedURLException {
    UrlValidator urlValidator = new UrlValidator();
    SamReaderFactory srf = SamReaderFactory.make();
    srf.validationStringency(ValidationStringency.SILENT);
    SamReader samReader;
    if (urlValidator.isValid(workFilename)) {
        samReader = srf.open(SamInputResource.of(new URL(workFilename)).index(new URL(workFilename + ".bai")));
    } else {
        samReader = srf.open(new File(workFilename));
    }
    return samReader;
}

Example 44 with SamReader

use of htsjdk.samtools.SamReader in project ASCIIGenome by dariober.

the class TrackReads method update.

/* M e t h o d s */
public void update() throws InvalidGenomicCoordsException, IOException {
    if (this.getyMaxLines() == 0) {
        return;
    }
    this.userWindowSize = this.getGc().getUserWindowSize();
    this.readStack = new ArrayList<List<SamSequenceFragment>>();
    if (this.getGc().getGenomicWindowSize() < this.MAX_REGION_SIZE) {
        SamReader samReader = Utils.getSamReader(this.getWorkFilename());
        List<Boolean> passFilter = this.filterReads(samReader, this.getGc().getChrom(), this.getGc().getFrom(), this.getGc().getTo());
        this.nRecsInWindow = 0;
        for (boolean x : passFilter) {
            // The count of reads in window is the count of reads passing filters
            if (x) {
                this.nRecsInWindow++;
            }
        }
        samReader = Utils.getSamReader(this.getWorkFilename());
        Iterator<SAMRecord> sam = samReader.query(this.getGc().getChrom(), this.getGc().getFrom(), this.getGc().getTo(), false);
        float max_reads = Float.parseFloat(Config.get(ConfigKey.max_reads_in_stack));
        float probSample = max_reads / this.nRecsInWindow;
        // Add this random String to the read name so different screenshot will generate
        // different samples.
        String rndOffset = Integer.toString(new Random().nextInt());
        List<TextRead> textReads = new ArrayList<TextRead>();
        ListIterator<Boolean> pass = passFilter.listIterator();
        while (sam.hasNext() && textReads.size() < max_reads) {
            SAMRecord rec = sam.next();
            if (pass.next()) {
                String templ_name = Utils.templateNameFromSamReadName(rec.getReadName());
                // Hashing.md5().hashBytes((templ_name + rndOffset).getBytes()).asLong();
                long v = (templ_name + rndOffset).hashCode();
                Random rand = new Random(v);
                if (rand.nextFloat() < probSample) {
                    // Downsampler
                    TextRead tr = new TextRead(rec, this.getGc(), Utils.asBoolean(Config.get(ConfigKey.show_soft_clip)));
                    textReads.add(tr);
                }
            }
        }
        this.readStack = stackReads(textReads);
    } else {
        this.nRecsInWindow = -1;
    }
}

Example 45 with SamReader

use of htsjdk.samtools.SamReader in project ASCIIGenome by dariober.

the class FilterTest method STUBcanCallBS.

// @Test
public void STUBcanCallBS() {
    System.out.println("\u203E");
    int f_incl = 0;
    int F_excl = 0;
    String chrom = "chrY";
    int from = 1;
    int to = 1;
    SamReaderFactory srf = SamReaderFactory.make();
    SamReader samReader = srf.open(new File("test_data/mjb050_oxBS.bam"));
    SAMFileHeader fh = samReader.getFileHeader();
    IntervalList il = new IntervalList(fh);
    Interval interval = new Interval(chrom, from, to);
    il.add(interval);
    List<SamRecordFilter> filters = FlagToFilter.flagToFilterList(f_incl, F_excl);
    SamLocusIterator samLocIter = new SamLocusIterator(samReader, il, true);
    samLocIter.setSamFilters(filters);
    while (samLocIter.hasNext()) {
        LocusInfo locus = samLocIter.next();
        int M = 0;
        int U = 0;
        int mism = 0;
        for (RecordAndOffset recOff : locus.getRecordAndPositions()) {
            int pos = locus.getPosition();
            // Code to get ref sequence at pos
            // If ref sequence is C, count read bases if:
            char refbase = Character.toUpperCase('\0');
            char rb = Character.toUpperCase((char) recOff.getReadBase());
            boolean isTopStrand = (new ReadFromTopStrandFilter(true)).filterOut(recOff.getRecord());
            if (refbase == 'C') {
                if (isTopStrand) {
                    // -ve 2nd pair
                    if (rb == 'C') {
                        M++;
                    } else if (rb == 'T') {
                        U++;
                    } else {
                        mism++;
                    }
                }
            } else if (refbase == 'G') {
                if (!isTopStrand) {
                    if (rb == 'G') {
                        M++;
                    } else if (rb == 'A') {
                        U++;
                    } else {
                        mism++;
                    }
                // System.out.println(locus.getPosition() + " ");
                }
            } else {
            // Not a C on ref
            }
        }
    // if(locus.getRecordAndPositions().size() > 0){
    // System.out.println(locus.getPosition() + " " + locus.getRecordAndPositions().size() + " " +
    // locus.getRecordAndPositions().get(0).getRecord().getFlags());
    // }
    }
}