Java SAMRecord类代码示例

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
private SAMRecord createSAMRecord(String read, String cigar, String md, int NM, boolean reverse) {
	SAMFileHeader header = new SAMFileHeader();
	List<SAMProgramRecord> PRs = new ArrayList<>();
	PRs.add(new SAMProgramRecord("Bowtie"));
	header.setProgramRecords(PRs);

	SAMRecord record = new SAMRecord(header);
	record.setReadString(read);
	record.setAttribute("NM", NM); //number of mismatches
	record.setCigarString(cigar);
	record.setAttribute("MD", md);

	if (reverse) {
		record.setReadNegativeStrandFlag(true);
	}

	return record;
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Test
public void testAlignmentOfRecord() throws Exception {
    SAMRecord record = new SAMRecord(this.header);
    record.setReadBases("TCCGTTAGTC".getBytes());
    record.setReadName(queryName);
    record.setCigarString("1S3M1I4M1S");
    record.setReferenceIndex(0);
    record.setReferenceName(this.referenceName);
    record.setAlignmentStart(2);

    assertEquals(referenceName, record.getReferenceName());

    Alignment result = SAMBEASTUtils.alignmentOfRecord(record, this.reference);
    assertEquals(2, result.getTaxonCount());
    assertEquals(new String(reference), result.getAlignedSequenceString(0));
    assertEquals("-CCGTAGT-", result.getAlignedSequenceString(1));
    assertEquals(referenceName, result.getTaxonId(0));
    assertEquals(queryName, result.getTaxonId(1));
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Override
public void run() {

    // count the number of poison pills received in queue
    int pills = 0;

    try {
        OneReadWithAlignment one;
        while (true) {
            one = queue.take();
            if (one.isOk()) {
                SAMRecord rec = one.getSAMRecord(samheader);
                alignSAM.addAlignment(rec);
            } else {
                pills++;
                if (pills >= poisonpills) {
                    break;
                }
            }
        }
    } catch (Exception ex) {
        System.out.println("Exception in readoutput: " + ex.getMessage());
    }

}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 * Collects information about SNV positions in a given alignment file.
 *
 *
 * @param snvpositions
 *
 * a set of genomic position with base changes
 *
 * @param bamf
 *
 * file to scan to fill in coverage information
 *
 * @param ginfo
 * @return
 *
 * a list of same length as snvposition in the input. This new list includes
 * counters for each base in the alphabet. The function scans the alignment
 * to fill in the counters.
 *
 * @throws IOException
 */
private SNVPositionDetailsList collectBAFDetailsFromBam(ArrayList<SNVPosition> snvpositions,
        File bamf, GenomeInfo ginfo) throws IOException {

    // from a list of positions, create a list that includes base/coverage counters
    SNVPositionDetailsList alldetails = makeDetailsList(snvpositions);

    // in the constructor for bamregions, I had the minimum mapping quality set at "minmapqual"
    // but this only worked well for minmapqual<=1. So now I set this to zero regardless of 
    // the mapping quality used in variant calling.        
    BamRegionsMap bamregions = new BamRegionsMap(bamf, ginfo, "SILENT", 0);

    int numloci = alldetails.size();
    for (int i = 0; i < numloci; i++) {
        SNVPositionDetails entry = alldetails.get(i);
        int nowchr = entry.getChrIndex();
        // look up the position in the bam
        SAMRecord[] varinbam = bamregions.lookup(ginfo.getChrName(nowchr), entry.getPosition());
        // collect information about the locus from each bam record
        countDetailsFromRecords(entry, varinbam);
    }
    bamregions.close();

    return (alldetails);
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 * Serves two purposes. 1 - to wrap the SAMRecord objects into
 * ThesaurusSAMRecords 2 - to clip the records if they have mismatches near
 * the ends. (But if the clipped records do not support the variant, the
 * clipping is not carried out)
 *
 *
 * @param varinbam
 * @param genomereader
 * @param cliplength
 * @return
 */
private ThesaurusSAMRecord[] makeClippedRecords(SAMRecord[] varinbam, FastaReader genomereader,
        int cliplength, int variantposition) {

    ThesaurusSAMRecord[] tbamrecords = new ThesaurusSAMRecord[varinbam.length];
    for (int k = 0; k < varinbam.length; k++) {
        SAMRecord record = varinbam[k];

        byte[] refsequence = genomereader.getSequenceBase1(record.getAlignmentStart(), record.getAlignmentEnd());

        ThesaurusSAMRecord trecord = new ThesaurusSAMRecord(record, refsequence, cliplength);
        // if the clipping cancel the variant position, revert the clipping                        
        if (trecord.getAlignmentStart() > variantposition || trecord.getAlignmentEnd() < variantposition) {
            trecord = new ThesaurusSAMRecord(record, refsequence, 0);
        }

        tbamrecords[k] = trecord;
    }
    return tbamrecords;
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 * make a copy of the current state of the regions object, into a primitive
 * array
 *
 * @return
 */
public SAMRecord[] getRegions() {

    // simple exit if no regions available
    if (isEmpty()) {
        return null;
    }

    // figure out how many regions are current stored in memory
    int count = getRegionsCount();

    // make an array of ThesaurusEntries here
    SAMRecord[] ans = new SAMRecord[count];
    int index = 0;
    for (Map.Entry<Integer, ArrayList<SAMRecord>> entry : regions.entrySet()) {
        ArrayList<SAMRecord> map1 = entry.getValue();
        for (int i = 0; i < map1.size(); i++) {
            ans[index] = map1.get(i);
            index++;
        }
    }        
    return ans;
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 * Converts a SAM record into a bed entry and outputs it straight away
 *
 * @param record
 */
private void processOneUnalignedRecord(SAMRecord record) {

    // output the alignment positions into a bed file
    int start = record.getAlignmentStart() - 1;
    int end = start + record.getReadLength();
    String out = record.getReferenceName() + "\t" + start + "\t" + end + "\n";

    synchronized (outunaligned) {
        try {
            outunaligned.write(out.getBytes());
        } catch (Exception ex) {
            System.out.println("Exception writing unaligned bed: " + ex.getMessage());
        }
    }

    // create a dummy thesaurus entry
    ThesaurusEntry entry = new ThesaurusEntry(record, ginfo);
    // complete it (need to fix alignment end manually, then use function to get anchors set right)
    entry.alignEnd = entry.originEnd;
    completeThesaurusEntry(entry, record);
    // set the origin chromosome to -1 to signal unaligned setup
    entry.originChrIndex = -1;
    outputThesaurusEntry(entry, outstreams.get(entry.getAlignChr()));

}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 * converts a SAM record into a thesaurus entry
 *
 * @param record
 * @param seqmap
 * @param outsream
 */
private void processOneRecord(SAMRecord record) {
    // create a new entry from the record, but give up if correct alignment
    ThesaurusEntry entry = new ThesaurusEntry(record, ginfo);
    if (entry.isTrivial()) {
        return;
    }

    if (isAlignmentClean(record)) {
        // fill in the missing information about mismatches in the entry
        completeThesaurusEntry(entry, record);
        outputThesaurusEntry(entry, outstreams.get(entry.getAlignChr()));
        // symmetrize the thesaurus by also outputing an entry with align<->origin
        ThesaurusEntry reverseentry = entry.getReverseEntry();
        outputThesaurusEntry(reverseentry, outstreams.get(reverseentry.getAlignChr()));
    }

}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 *
 * @return
 *
 * true if the cigar contains only one element that describes a sequence
 * match, e.g. 100M
 *
 *
 */
private boolean isAlignmentClean(SAMRecord record) {

    Cigar cigar = record.getCigar();
    int numcigel = cigar.numCigarElements();
    if (numcigel > 1) {
        return false;
    }
    CigarOperator co = cigar.getCigarElement(0).getOperator();
    if (co != CigarOperator.M) {
        return false;
    }

    return true;
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
/**
 * copy reads from one bam file and put them into files by chromosome.
 *
 * @param nowbam
 */
private void splitIntoChroms(File nowbam) {
    SAMFileReader inputSam = new SAMFileReader(nowbam);

    // write each aligned record into its separate file, or into "unmapped"
    for (final SAMRecord record : inputSam) {
        String nowchrom = record.getReferenceName();
        SAMFileWriter nowwriter = chrombams.get(nowchrom);
        if (nowwriter == null) {
            nowwriter = chrombams.get("unmapped");
        }
        nowwriter.addAlignment(record);
    }

    inputSam.close();
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Test
public void testProcessSequenceRaise() {
	System.out.println("Running test InsertLengthDistributionTest.testProcessSequenceRaise");
	log.info("Running test InsertLengthDistributionTest.testProcessSequenceRaise");

	testObjectFactory = new TestObjectFactory();
	samRecords = testObjectFactory.getSamRecords();
	
	samRecords.get(0).setInferredInsertSize(7);
	samRecords.get(1).setInferredInsertSize(3);
	samRecords.get(2).setInferredInsertSize(1);

	for (SAMRecord samRecord : samRecords) {
		insertLengthDistribution.processSequence(samRecord);
	}
	assertTrue(insertLengthDistribution.raisesError());
	assertTrue(insertLengthDistribution.raisesWarning());
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Test
public void testProcessSequenceRaiseCalculation() {

	System.out.println("Running test InsertLengthDistributionTest.testProcessSequenceRaiseCalculation");
	log.info("Running test InsertLengthDistributionTest.testProcessSequenceRaiseCalculation");
		
	List<Long> insertSizes = UtilityTest.readInsertSizesLong();
	int index = 0;
	SAMFileHeader samFileHeader = TestObjectFactory.getInstance();
	SAMRecord samRecord = new SAMRecord(samFileHeader);

	samRecord.setProperPairFlag(true);
	samRecord.setReadPairedFlag(true);

	for (long count : insertSizes) {
		for (int i = 0; i < count; i++) {
			samRecord.setInferredInsertSize(index);
			
			insertLengthDistribution.processSequence(samRecord);
		}
		index++;
	}
	assertFalse(insertLengthDistribution.raisesError());
	assertFalse(insertLengthDistribution.raisesWarning());
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Test
public void testBooleans() {
	System.out.println("Running test InsertLengthDistributionTest.testBooleans");	
	log.info("Running test InsertLengthDistributionTest.testBooleans");
	
	testObjectFactory = new TestObjectFactory();
	samRecords = testObjectFactory.getSamRecords();
	samRecords.get(0).setInferredInsertSize(7);
	samRecords.get(1).setInferredInsertSize(3);
	samRecords.get(2).setInferredInsertSize(1);
	for (SAMRecord samRecord : samRecords) {
		insertLengthDistribution.processSequence(samRecord);
	}
	
	assertFalse(insertLengthDistribution.ignoreInReport());
	assertFalse(insertLengthDistribution.needsToSeeAnnotation());
	assertTrue(insertLengthDistribution.raisesError());
	assertTrue(insertLengthDistribution.raisesWarning());
	assertTrue(insertLengthDistribution.needsToSeeSequences());
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Before
public void setUp() throws Exception {	
	variantCallDetection = new VariantCallDetection();
	basicStatistics = new BasicStatistics(variantCallDetection);
	
	String filename = new String(new File("").getAbsolutePath() + "/test/resources/test_header.sam");
	samRecords = SAMRecordLoader.loadSAMFile(filename);		
	if(samRecords.isEmpty()) { 
		log.warn("Impossible to run the test as " + filename + " seems empty");
		return; 
	}

	for(SAMRecord read : samRecords) {
		variantCallDetection.processSequence(read);
		basicStatistics.processSequence(read);
	}
}
 

import net.sf.samtools.SAMRecord; //导入依赖的package包/类
@Test
public void testProcessSequence() {
	System.out.println("Running test SequenceQualityDistributionTest.testProcessSequence");
	log.info("Running test SequenceQualityDistributionTest.testProcessSequence");
	
	for (SAMRecord samRecord : samRecords) {
		sequenceQualityDistribution.processSequence(samRecord);
	}
	List<Integer> distribution = sequenceQualityDistribution.getDistribution();
	
	assertEquals(0, (int) distribution.get(0));
	assertEquals(0, (int) distribution.get(1));
	assertEquals(0, (int) distribution.get(2));
	assertEquals(2, (int) distribution.get(3));
	assertEquals(0, (int) distribution.get(4));
	assertEquals(0, (int) distribution.get(5));
	assertEquals(1, (int) distribution.get(6));
}