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

Example 36 with SamReader

use of htsjdk.samtools.SamReader in project polyGembler by c-zhou.

the class GenomeComparison method run.

@Override
public void run() {
    // TODO Auto-generated method stub
    final SamReaderFactory factory = SamReaderFactory.makeDefault().enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS, SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS).validationStringency(ValidationStringency.SILENT);
    final SamReader inputSam1 = factory.open(new File(in_bam1));
    SAMRecordIterator iter1 = inputSam1.iterator();
    final SamReader inputSam2 = factory.open(new File(in_bam2));
    SAMRecordIterator iter2 = inputSam2.iterator();
    this.initial_thread_pool();
    try {
        bw_con = new BufferedWriter(new FileWriter(out_f + ".txt"));
        bw_1 = new BufferedWriter(new FileWriter(out_f + ".1txt"));
        bw_2 = new BufferedWriter(new FileWriter(out_f + ".2txt"));
        SAMRecord record1 = iter1.next(), record2 = iter2.next();
        List<SAMRecord> barcoded_records1 = new ArrayList<SAMRecord>(), barcoded_records2 = new ArrayList<SAMRecord>();
        barcoded_records1.add(record1);
        barcoded_records2.add(record2);
        String barcode1 = record1.getStringAttribute("BX"), barcode2 = record2.getStringAttribute("BX");
        jobs: while (true) {
            if (barcode1.compareTo(barcode2) < 0) {
                // output barcode 1 up to barcode 1
                while (iter1.hasNext()) {
                    while ((record1 = iter1.next()).getStringAttribute("BX").equals(barcode1)) {
                        barcoded_records1.add(record1);
                        if (!iter1.hasNext())
                            break;
                    }
                    // TODO barcoded records 1 output
                    Molecule[] mols = extractMoleculeFromList(barcoded_records1);
                    for (int i = 0; i != mols.length; i++) bw_1.write(mols[i].chr_id + ":" + mols[i].chr_start + "-" + mols[i].chr_end + "\n");
                    barcoded_records1.clear();
                    if (iter1.hasNext()) {
                        barcoded_records1.add(record1);
                        barcode1 = record1.getStringAttribute("BX");
                    } else
                        break jobs;
                    if (barcode1.compareTo(barcode2) >= 0)
                        break;
                }
                if (!barcode1.equals(barcode2))
                    continue;
            } else if (barcode1.compareTo(barcode2) > 0) {
                // output barcode 2 up to barcode 2
                while (iter2.hasNext()) {
                    while ((record2 = iter2.next()).getStringAttribute("BX").equals(barcode2)) {
                        barcoded_records2.add(record2);
                        if (!iter2.hasNext())
                            break;
                    }
                    // TODO barcoded records 2 output
                    Molecule[] mols = extractMoleculeFromList(barcoded_records2);
                    for (int i = 0; i != mols.length; i++) bw_2.write(mols[i].chr_id + ":" + mols[i].chr_start + "-" + mols[i].chr_end + "\n");
                    barcoded_records2.clear();
                    if (iter2.hasNext()) {
                        barcoded_records2.add(record2);
                        barcode2 = record2.getStringAttribute("BX");
                    } else
                        break jobs;
                    if (barcode2.compareTo(barcode1) >= 0)
                        break;
                }
                if (!barcode1.equals(barcode2))
                    continue;
            }
            while ((record1 = iter1.next()).getStringAttribute("BX").equals(barcode1)) {
                barcoded_records1.add(record1);
                if (!iter1.hasNext())
                    break;
            }
            while ((record2 = iter2.next()).getStringAttribute("BX").equals(barcode2)) {
                barcoded_records2.add(record2);
                if (!iter2.hasNext())
                    break;
            }
            executor.submit(new Runnable() {

                List<SAMRecord> list1;

                List<SAMRecord> list2;

                @Override
                public void run() {
                    // TODO Auto-generated method stub
                    try {
                        Molecule[] mols1 = extractMoleculeFromList(list1);
                        Molecule[] mols2 = extractMoleculeFromList(list2);
                        int n1 = mols1.length, n2 = mols2.length;
                        double insec;
                        Set<Integer> m1 = new HashSet<Integer>();
                        Set<Integer> m2 = new HashSet<Integer>();
                        for (int i = 0; i != n1; i++) {
                            for (int j = 0; j != n2; j++) {
                                if ((insec = intersect(mols1[i], mols2[j])) >= overlap_frac) {
                                    bw_con.write(mols1[i].chr_id + ":" + mols1[i].chr_start + "-" + mols1[i].chr_end + "\t" + mols2[j].chr_id + ":" + mols2[j].chr_start + "-" + mols2[j].chr_end + "\t" + insec + "\n");
                                    m1.add(i);
                                    m2.add(j);
                                }
                            }
                        }
                        for (int i = 0; i != n1; i++) {
                            if (!m1.contains(i))
                                bw_1.write(mols1[i].chr_id + ":" + mols1[i].chr_start + "-" + mols1[i].chr_end + "\n");
                        }
                        for (int i = 0; i != n2; i++) {
                            if (!m2.contains(i))
                                bw_2.write(mols2[i].chr_id + ":" + mols2[i].chr_start + "-" + mols2[i].chr_end + "\n");
                        }
                    } catch (Exception e) {
                        Thread t = Thread.currentThread();
                        t.getUncaughtExceptionHandler().uncaughtException(t, e);
                        e.printStackTrace();
                        executor.shutdown();
                        System.exit(1);
                    }
                }

                public Runnable init(List<SAMRecord> list1, List<SAMRecord> list2) {
                    // TODO Auto-generated method stub
                    this.list1 = list1;
                    this.list2 = list2;
                    return this;
                }
            }.init(new ArrayList<SAMRecord>(barcoded_records1), new ArrayList<SAMRecord>(barcoded_records2)));
            if (!iter1.hasNext() || !iter2.hasNext())
                break jobs;
            barcoded_records1.clear();
            barcoded_records2.clear();
            barcoded_records1.add(record1);
            barcoded_records2.add(record2);
            barcode1 = record1.getStringAttribute("BX");
            barcode2 = record2.getStringAttribute("BX");
        }
        if (!iter1.hasNext()) {
            // TODO output remaining iter2
            while (iter2.hasNext()) {
                while ((record2 = iter2.next()).getStringAttribute("BX").equals(barcode2)) {
                    barcoded_records2.add(record2);
                    if (!iter2.hasNext())
                        break;
                }
                Molecule[] mols = extractMoleculeFromList(barcoded_records2);
                for (int i = 0; i != mols.length; i++) bw_2.write(mols[i].chr_id + ":" + mols[i].chr_start + "-" + mols[i].chr_end + "\n");
                barcoded_records2.clear();
                if (iter2.hasNext()) {
                    barcoded_records2.add(record2);
                    barcode2 = record2.getStringAttribute("BX");
                }
            }
        }
        if (!iter2.hasNext()) {
            // TODO output remaining iter1
            while (iter1.hasNext()) {
                while ((record1 = iter1.next()).getStringAttribute("BX").equals(barcode1)) {
                    barcoded_records1.add(record1);
                    if (!iter1.hasNext())
                        break;
                }
                Molecule[] mols = extractMoleculeFromList(barcoded_records1);
                for (int i = 0; i != mols.length; i++) bw_1.write(mols[i].chr_id + ":" + mols[i].chr_start + "-" + mols[i].chr_end + "\n");
                barcoded_records1.clear();
                if (iter1.hasNext()) {
                    barcoded_records1.add(record1);
                    barcode1 = record1.getStringAttribute("BX");
                }
            }
        }
        this.waitFor();
        inputSam1.close();
        inputSam2.close();
        bw_con.close();
        bw_1.close();
        bw_2.close();
    } catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
}

Example 37 with SamReader

use of htsjdk.samtools.SamReader in project polyGembler by c-zhou.

the class TenXMoleculeStats method runSNPCaller.

private void runSNPCaller() {
    // TODO Auto-generated method stub
    final SamReader inputSam = factory.open(new File(this.in_bam));
    if (inputSam.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate)
        throw new RuntimeException("Sort BAM file before SNP calling!!!");
    final SAMSequenceDictionary seqDic = inputSam.getFileHeader().getSequenceDictionary();
    mappedReadsOnChromosome = new long[seqDic.size()];
    oos = Utils.getBufferedWriter(out_stats + ".txt");
    hap = Utils.getBufferedWriter(out_stats + ".haps");
    vcf = Utils.getBufferedWriter(out_stats + ".snps");
    SAMFileHeader header = inputSam.getFileHeader();
    header.setSortOrder(SortOrder.unknown);
    oob = new SAMFileWriterFactory().makeSAMOrBAMWriter(header, true, new File(out_stats + ".bam"));
    try {
        inputSam.close();
    } catch (IOException e1) {
        // TODO Auto-generated catch block
        e1.printStackTrace();
    }
    this.initial_thread_pool();
    for (int i = 0; i < seqDic.size(); i++) {
        executor.submit(new Runnable() {

            private int chr_num;

            @Override
            public void run() {
                // TODO Auto-generated method stub
                try {
                    final SAMSequenceRecord refSeq = seqDic.getSequence(chr_num);
                    final String seqnm = refSeq.getSequenceName();
                    final int seqln = refSeq.getSequenceLength();
                    final Set<Integer> snp_pos = readSNPs(in_vcf, seqnm);
                    final SamReader inputSam = factory.open(new File(in_bam));
                    final SAMRecordIterator iter = inputSam.queryOverlapping(seqnm, 0, seqln);
                    // do NOT call SNPs, read it from file instead
                    final Map<String, List<SAMRecord>> bc_record = new HashMap<String, List<SAMRecord>>();
                    final Map<String, Integer> bc_position = new HashMap<String, Integer>();
                    String bc;
                    SAMRecord tmp_record;
                    int pos;
                    while (iter.hasNext()) {
                        tmp_record = iter.next();
                        if (!processRecord(tmp_record))
                            continue;
                        bc = tmp_record.getStringAttribute("BX");
                        if (bc == null)
                            continue;
                        pos = tmp_record.getAlignmentStart();
                        if (mappedReads % 10000 == 0) {
                            Set<String> bc_removal = new HashSet<String>();
                            for (Map.Entry<String, List<SAMRecord>> entry : bc_record.entrySet()) {
                                String b = entry.getKey();
                                List<SAMRecord> records = entry.getValue();
                                if (bc_position.get(b) + max_gap < pos) {
                                    processMolecule(b, records, snp_pos);
                                    bc_removal.add(b);
                                }
                            }
                            for (String b : bc_removal) {
                                bc_record.remove(b);
                                bc_position.remove(b);
                            }
                        }
                        if (bc_record.containsKey(bc)) {
                            List<SAMRecord> records = bc_record.get(bc);
                            if (bc_position.get(bc) + max_gap < pos) {
                                processMolecule(bc, records, snp_pos);
                                records.clear();
                                bc_position.remove(bc);
                            }
                        } else
                            bc_record.put(bc, new ArrayList<SAMRecord>());
                        bc_record.get(bc).add(tmp_record);
                        bc_position.put(bc, Math.max(tmp_record.getAlignmentEnd(), bc_position.containsKey(bc) ? bc_position.get(bc) : 0));
                    }
                    if (!bc_record.isEmpty())
                        for (Map.Entry<String, List<SAMRecord>> entry : bc_record.entrySet()) processMolecule(entry.getKey(), entry.getValue(), snp_pos);
                    /**
                     *						// do NOT call SNPs, read it from file instead
                     *						int pos;
                     *						SAMRecord buffer = iter.hasNext() ? iter.next() : null;
                     *						int bufferS = buffer==null ? Integer.MAX_VALUE : buffer.getAlignmentStart();
                     *						final Set<SAMRecord> record_pool = new HashSet<SAMRecord>();
                     *						final Set<SAMRecord> tmp_removal = new HashSet<SAMRecord>();
                     *
                     *						final int[] allele_stats = new int[5];
                     *
                     *						final Map<String, TreeSet<SNP>> bc_snp = new HashMap<String, TreeSet<SNP>>();
                     *
                     *						String dna_seq, bc;
                     *						int tmp_int;
                     *						String snp_str;
                     *						char nucl;
                     *						for(int p=0; p!=seqln; p++) {
                     *							if(p%1000000==0) {
                     *								System.err.println(seqnm+":"+p);
                     *								// TODO
                     *								// process bc_snp
                     *								Set<String> bc_removal = new HashSet<String>();
                     *								for(Map.Entry<String, TreeSet<SNP>> entry : bc_snp.entrySet()) {
                     *									if(entry.getValue().last().position+max_gap<p) {
                     *										processMolecule(entry.getKey(), entry.getValue());
                     *										bc_removal.add(entry.getKey());
                     *									}
                     *								}
                     *								for(String b : bc_removal) bc_snp.remove(b);
                     *							}
                     *							pos = p+1; // BAM format is 1-based coordination
                     *
                     *							while(bufferS<=pos) {
                     *								// TODO
                     *								// buffer records upto this position
                     *								while( buffer!=null &&
                     *										buffer.getAlignmentStart()<=pos) {
                     *									if( processRecord(buffer) )
                     *										record_pool.add(buffer);
                     *									buffer=iter.hasNext()?iter.next():null;
                     *									bufferS = buffer==null ? Integer.MAX_VALUE :
                     *										buffer.getAlignmentStart();
                     *								}
                     *							}
                     *							if(record_pool.isEmpty()) continue;
                     *
                     *							Arrays.fill(allele_stats, 0);
                     *							tmp_removal.clear();
                     *							for(SAMRecord record : record_pool) {
                     *
                     *								if(record.getAlignmentEnd()<pos) {
                     *									tmp_removal.add(record);
                     *									continue;
                     *								}
                     *								dna_seq = record.getReadString();
                     *								tmp_int = record.getReadPositionAtReferencePosition(pos);
                     *								nucl = tmp_int==0?'D':Character.toUpperCase(dna_seq.charAt(tmp_int-1));
                     *								switch(nucl) {
                     *								case 'A':
                     *									allele_stats[0]++;
                     *									break;
                     *								case 'C':
                     *									allele_stats[1]++;
                     *									break;
                     *								case 'G':
                     *									allele_stats[2]++;
                     *									break;
                     *								case 'T':
                     *									allele_stats[3]++;
                     *									break;
                     *								case 'D':
                     *									allele_stats[4]++;
                     *									break;
                     *								case 'N':
                     *									// do nothing here
                     *									break;
                     *								default:
                     *									throw new RuntimeException("!!!");
                     *								}
                     *							}
                     *							record_pool.removeAll(tmp_removal);
                     *
                     *							snp_str = feasibleSNP(allele_stats);
                     *
                     *							// not two alleles or minor allele frequency is smaller than threshold
                     *							// don't call indels
                     *							if(snp_str==null) continue;
                     *
                     *							vcf.write(seqnm+"\t"+p+"\t"+snp_str+"\n");
                     *
                     *							for(SAMRecord record : record_pool) {
                     *								bc = record.getStringAttribute("BX");
                     *								if(bc==null) continue;
                     *								dna_seq = record.getReadString();
                     *								tmp_int = record.getReadPositionAtReferencePosition(pos);
                     *								if(tmp_int==0) continue;
                     *								nucl = Character.toUpperCase(dna_seq.charAt(tmp_int-1));
                     *								if(bc_snp.containsKey(bc)) {
                     *									TreeSet<SNP> bc_set = bc_snp.get(bc);
                     *									if(bc_set.last().position+max_gap<p) {
                     *										// TODO
                     *										// gap exceeds max_gap
                     *										// write molecule out
                     *										processMolecule(bc, bc_set);
                     *										bc_set.clear();
                     *									}
                     *								} else {
                     *									bc_snp.put(bc, new TreeSet<SNP>(new SNP.PositionComparator()));
                     *								}
                     *								bc_snp.get(bc).add(new SNP(p, nucl, record));
                     *
                     *							}
                     *						}
                     *
                     *						for(Map.Entry<String, TreeSet<SNP>> entry : bc_snp.entrySet())
                     *							processMolecule(entry.getKey(), entry.getValue());
                     */
                    iter.close();
                    inputSam.close();
                } catch (Exception e) {
                    Thread t = Thread.currentThread();
                    t.getUncaughtExceptionHandler().uncaughtException(t, e);
                    e.printStackTrace();
                    executor.shutdown();
                    System.exit(1);
                }
            }

            private void processMolecule(final String bc, final List<SAMRecord> records, final Set<Integer> snp_pos) {
                // TODO Auto-generated method stub
                MoleculeConstructor molCon = new MoleculeConstructor(records);
                if (molCon.getMolecularLength() < oo_thresh)
                    return;
                molCon.run();
            }

            /**
             *				private void processMolecule(final String bc, final TreeSet<SNP> snp_set)
             *						throws IOException {
             *					// TODO Auto-generated method stub
             *					// write molecule haplotypes
             *					Set<SAMRecord> records = new HashSet<SAMRecord>();
             *					List<String> snp_str = new ArrayList<String>();
             *					SNP snp;
             *					int tmp_position = -1;
             *					while(!snp_set.isEmpty()) {
             *						snp = snp_set.pollFirst();
             *						records.add(snp.alignment);
             *						if(tmp_position==snp.position) {
             *							snp_str.remove(snp_str.size()-1);
             *						} else {
             *							snp_str.add(snp.position+","+
             *									snp.allele);
             *							tmp_position = snp.position;
             *						}
             *					}
             *					// write molecule statistics
             *					MoleculeConstructor molCon  =
             *							new MoleculeConstructor(new ArrayList<SAMRecord>(records));
             *					if(molCon.getMolecularLength()<oo_thresh) return;
             *					molCon.run();
             *
             *					if(snp_str.size()<2) return;
             *					StringBuilder os = new StringBuilder();
             *					os.append(bc);
             *					for(String ss : snp_str) {
             *						os.append(" ");
             *						os.append( ss);
             *					}
             *					os.append("\n");
             *					hap.write(os.toString());
             *				}
             */
            public Runnable init(final int i) {
                // TODO Auto-generated method stub
                this.chr_num = i;
                return this;
            }

            private String feasibleSNP(int[] ints) {
                // TODO Auto-generated method stub
                if (ints[4] > 0)
                    return null;
                List<Integer> obs = new ArrayList<Integer>();
                for (int i = 0; i < 4; i++) if (ints[i] > 0)
                    obs.add(i);
                if (obs.size() != 2)
                    return null;
                int ref = obs.get(0), alt = obs.get(1);
                int d = ints[ref] + ints[alt];
                if (d < min_depth)
                    return null;
                double maf = (double) ints[ref] / d;
                if (maf > 0.5)
                    maf = 1 - maf;
                if (maf < min_maf)
                    return null;
                return Sequence.nucleotide[ref] + "," + Sequence.nucleotide[alt] + "\t" + ints[ref] + "," + ints[alt];
            }
        }.init(i));
    }
    try {
        inputSam.close();
        this.waitFor();
        oos.close();
        hap.close();
        vcf.close();
        oob.close();
        BufferedWriter bw_summary = Utils.getBufferedWriter(this.out_stats + ".summary");
        bw_summary.write("##Reads     mapped: " + mappedReads + "\n");
        bw_summary.write("##Reads   unmapped: " + unmappedReads + "\n");
        bw_summary.write("##Reads duplicated: " + duplicatedReads + "\n");
        bw_summary.write("##Reads  bad pairs: " + badpaired + "\n");
        bw_summary.write("##Reads by pseudo-chromosomes:\n");
        for (int i = 0; i < mappedReadsOnChromosome.length; i++) {
            bw_summary.write(" #" + seqDic.getSequence(i).getSequenceName() + " " + mappedReadsOnChromosome[i] + "\n");
        }
        bw_summary.write("##Reads by barcodes:\n");
        for (Map.Entry<String, Integer> entry : readsOnBarcode.entrySet()) {
            bw_summary.write(" #" + entry.getKey() + " " + entry.getValue() + "\n");
        }
        bw_summary.close();
    } catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
}

Example 38 with SamReader

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

the class SamLocusIteratorTest method testBasicIterator.

@Test
public void testBasicIterator() {
    final String sqHeader = "@HD\tSO:coordinate\tVN:1.0\n@SQ\tSN:chrM\tAS:HG18\tLN:100000\n";
    final String seq1 = "ACCTACGTTCAATATTACAGGCGAACATACTTACTA";
    // phred 10
    final String qual1 = "*";
    final String s1 = "3851612\t16\tchrM\t165\t255\t36M\t*\t0\t0\t" + seq1 + "\t" + qual1 + "\n";
    final String exampleSam = sqHeader + s1 + s1;
    final SamReader samReader = createSamFileReader(exampleSam);
    final SamLocusIterator sli = createSamLocusIterator(samReader);
    // make sure we accumulated depth of 2 for each position
    int pos = 165;
    for (final SamLocusIterator.LocusInfo li : sli) {
        assertEquals(pos++, li.getPosition());
        assertEquals(2, li.getRecordAndPositions().size());
    }
}

Example 39 with SamReader

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

the class SamLocusIteratorTest method testIteratorWithIntron.

@Test
public void testIteratorWithIntron() {
    // See https://github.com/samtools/htsjdk/issues/838
    // Prepare a sam header
    String sqHeader = "@HD\tSO:coordinate\tVN:1.0\n" + "@SQ\tSN:chr1\tAS:HG18\tLN:10000000\n";
    // Prepare one read with a 500,000 bases skipped
    String cigar = "18M500000N18M";
    String s1 = "read1\t0\tchr1\t1\t255\t" + cigar + "\t*\t0\t0\tACCTACGTTCAATATTACAGGCGAACATACTTACTA\t*\n";
    // Prepare sam input and samReader
    String exampleSam = sqHeader + s1 + s1;
    ByteArrayInputStream inputStream = new ByteArrayInputStream(exampleSam.getBytes());
    SamReader samReader = SamReaderFactory.makeDefault().open(SamInputResource.of(inputStream));
    // A small interval to iterate over:
    IntervalList il = new IntervalList(samReader.getFileHeader());
    il.add(new Interval("chr1", 1, 100));
    SamLocusIterator sli = new SamLocusIterator(samReader, il, true);
    // Iterate
    long t0 = System.currentTimeMillis();
    int n = 0;
    for (SamLocusIterator.LocusInfo li : sli) {
        n++;
    }
    long t1 = System.currentTimeMillis();
    System.err.println("Time to iterate " + n + " loci: " + (t1 - t0) / 1000.0 + " sec");
    sli.close();
}

Example 40 with SamReader

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

the class GenomicCoords method setSamSeqDictFromBam.

private boolean setSamSeqDictFromBam(String bamfile) {
    /*  ------------------------------------------------------ */
    /* This chunk prepares SamReader from local bam            */
    SamReaderFactory srf = SamReaderFactory.make();
    srf.validationStringency(ValidationStringency.SILENT);
    SamReader samReader;
    samReader = srf.open(new File(bamfile));
    /*  ------------------------------------------------------ */
    SAMSequenceDictionary seqDict = samReader.getFileHeader().getSequenceDictionary();
    if (seqDict != null && !seqDict.isEmpty()) {
        this.setSamSeqDictSource(new File(bamfile).getAbsolutePath());
        this.setSamSeqDict(seqDict);
        return true;
    }
    return false;
}