Java SAMRecord.getReadLength方法代码示例

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * HACK TO CREATE GATKSAMRECORD BASED ONLY A SAMRECORD FOR TESTING PURPOSES ONLY
 * @param read
 */
public GATKSAMRecord(final SAMRecord read) {
    super(read.getHeader());
    setReferenceIndex(read.getReferenceIndex());
    setAlignmentStart(read.getAlignmentStart());
    mReadNameLength = (short)read.getReadNameLength();
    setMappingQuality(read.getMappingQuality());
    mCigarLength = read.getCigarLength();
    setFlags(read.getFlags());
    mReadLength = read.getReadLength();
    setMateReferenceIndex(read.getMateReferenceIndex());
    setMateAlignmentStart(read.getMateAlignmentStart());
    setInferredInsertSize(read.getInferredInsertSize());
    mRestOfBinaryData = null;

    super.setReadName(read.getReadName());
    super.setCigarString(read.getCigarString());
    super.setReadBases(read.getReadBases());
    super.setBaseQualities(read.getBaseQualities());

    SAMReadGroupRecord samRG = read.getReadGroup();
    clearAttributes();
    if (samRG != null) {
        GATKSAMReadGroupRecord rg = new GATKSAMReadGroupRecord(samRG);
        setReadGroup(rg);
    }
    List<SAMTagAndValue> attributes = read.getAttributes();
    for(SAMTagAndValue tagAndValue : attributes) {
        setAttribute(tagAndValue.tag, tagAndValue.value);
    }
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
public void loadReads(File[] bams) throws IOException{

int count = 0;
int numOp = 0;

for(File bam : bams){
    HLA.log.appendln("Loading reads from:\t" + bam.getName());
    Object2IntOpenHashMap<String> readLoadingSet = new Object2IntOpenHashMap<String>();
    readLoadingSet.defaultReturnValue(0);
    
    final SamReader reader = SamReaderFactory.makeDefault().open(bam);
    
    //Kourami bam checker added
    if(!checkHeader(reader.getFileHeader())){
	HLA.log.appendln("Unexpected BAM :\t"+ bam.getName() 
			+"\nThe input BAM MUST be aligned to the set of IMGT/HLA alleles in " + HLA.MSAFILELOC + "\n" 
			+ "Please use the recommended preprocessing steps explained on the github page:\n"
			+ "https://github.com/Kingsford-Group/kourami");
	System.err.println("Unexpected BAM :\t"+ bam.getName() 
			   +"\nThe input BAM MUST be aligned to the set of IMGT/HLA alleles in " + HLA.MSAFILELOC + "\n" 
			   + "Please use the recommended preprocessing steps explained on the github page:\n"
			   + "https://github.com/Kingsford-Group/kourami");
	HLA.log.outToFile();
	System.exit(1);
    }

    for(final SAMRecord samRecord : reader){
	if(count == 0){
	    HLA.READ_LENGTH = samRecord.getReadLength();
	    HLA.log.appendln("Setting HLA.READ_LEGNTH = " + HLA.READ_LENGTH);
	}
	//added checking to process reads matching to HLA-type sequences
	//discarding decoy hits (DQB2, DQA2)
	boolean qc = false;
	if( (samRecord.getReferenceName().indexOf("*") > -1) 
	    && !samRecord.getReadUnmappedFlag() 
	    && !samRecord.isSecondaryOrSupplementary() 
	    && !this.startWIns(samRecord)){
	    count++;
	    if(samRecord.getReadPairedFlag())
		numOp += processRecord(samRecord, readLoadingSet);
	    else
		numOp += processRecordUnpaired(samRecord);
	}
	
	if(HLA.DEBUG && count%10000 == 0)
	    HLA.log.appendln("Processed 10000 reads...");
    }
    reader.close();
}
HLA.log.appendln("Loaded a total of " + count + " mapped reads.");
HLA.log.appendln("A total of " + numOp + " bases");
   }
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
public boolean qcCheck(SAMRecord sr){
Cigar cigar = sr.getCigar();
int rLen = sr.getReadLength();
int effectiveLen = 0;
if(cigar==null) 
    return false;
else{
    for(final CigarElement ce : cigar.getCigarElements()){
	CigarOperator op = ce.getOperator();
	int cigarLen = ce.getLength();
	switch(op)
	    {
	    case M:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    case I:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    default: 
		break;
	    }
    }
}
boolean readdebug = false;
if(readdebug){
    HLA.log.appendln(sr.getSAMString());
    HLA.log.appendln("EffectiveLen:\t" + effectiveLen);
    HLA.log.appendln("ReadLen:\t" + rLen);
}
Integer i = sr.getIntegerAttribute("NM");
int nm = 0;
if(i!=null)
    nm = i.intValue();
if(readdebug)
    HLA.log.appendln("NM=\t" + nm);
if(nm < 16){
    if(readdebug)
	HLA.log.appendln("PASSWED QC");
    return true;
}
if(readdebug){
    HLA.log.appendln("FAILED QC");
    HLA.log.appendln(sr.getSAMString());
}
return false;
   }
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
public BAQCalculationResult calcBAQFromHMM(SAMRecord read, byte[] ref, int refOffset) {
    // todo -- need to handle the case where the cigar sum of lengths doesn't cover the whole read
    Pair<Integer, Integer> queryRange = calculateQueryRange(read);
    if ( queryRange == null ) return null; // read has Ns, or is completely clipped away

    int queryStart = queryRange.getFirst();
    int queryEnd = queryRange.getSecond();

    BAQCalculationResult baqResult = calcBAQFromHMM(ref, read.getReadBases(), read.getBaseQualities(), queryStart, queryEnd);

    // cap quals
    int readI = 0, refI = 0;
    for ( CigarElement elt : read.getCigar().getCigarElements() ) {
        int l = elt.getLength();
        switch (elt.getOperator()) {
            case N: // cannot handle these
                return null;
            case H : case P : // ignore pads and hard clips
                break;
            case S : refI += l; // move the reference too, in addition to I
            case I :
                // todo -- is it really the case that we want to treat I and S the same?
                for ( int i = readI; i < readI + l; i++ ) baqResult.bq[i] = baqResult.rawQuals[i];
                readI += l;
                break;
            case D : refI += l; break;
            case M :
                for (int i = readI; i < readI + l; i++) {
                    int expectedPos = refI - refOffset + (i - readI);
                    baqResult.bq[i] = capBaseByBAQ( baqResult.rawQuals[i], baqResult.bq[i], baqResult.state[i], expectedPos );
                }
                readI += l; refI += l;
                break;
            default:
                throw new ReviewedGATKException("BUG: Unexpected CIGAR element " + elt + " in read " + read.getReadName());
        }
    }
    if ( readI != read.getReadLength() ) // odd cigar string
        System.arraycopy(baqResult.rawQuals, 0, baqResult.bq, 0, baqResult.bq.length);

    return baqResult;
}