Java SAMRecord.getMateAlignmentStart方法代码示例

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
public static SAMRecord convertToPairedEndFragment(SAMRecord rec) {
	int insertSize = Math.abs(rec.getInferredInsertSize());
	//TODO: possible to remove (insertSize > 0) clause? What does an insert size of 0 mean?
	if (rec.getReadPairedFlag() && rec.getProperPairFlag() && !rec.getReadUnmappedFlag() 
			&& (insertSize <= MAX_INSERT) && (rec.getAlignmentStart() < rec.getMateAlignmentStart())
			&& rec.getReferenceName() != "chrM") // current paired end representation doesn't do well with circular chromosomes
	{
		//We need to get the full fragment contained by read.getStart to pair.getEnd
		int readEnd=rec.getAlignmentEnd();
		int mateStart=rec.getMateAlignmentStart();
					
		//int extension = insertSize - rec.getReadLength();
		int extension=(mateStart-readEnd)+rec.getReadLength();
		if (extension <= MAX_INSERT) {
			String newRead = rec.getReadString() + StringUtils.repeat("N", extension);
			rec.setReadString(newRead);
			String newQual = StringUtils.repeat("A", newRead.length());
			if(!rec.getBaseQualityString().equals("*")) newQual = rec.getBaseQualityString() + StringUtils.repeat("A", extension);
			rec.setBaseQualityString(newQual);
			String newCigar=newRead.length()+"M";
			rec.setCigarString(newCigar);

			// Change attributes to represent single read
			rec.setMateReferenceName("*");
			rec.setMateAlignmentStart(0);
			rec.setFirstOfPairFlag(false);
			rec.setMateNegativeStrandFlag(false);
			rec.setMateUnmappedFlag(false);
			rec.setProperPairFlag(false);
			rec.setReadPairedFlag(false);
			rec.setSecondOfPairFlag(false);
		}
	} else {
		rec = null;
	}
	return rec;
}
 

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
private boolean isCompatiblePair(SAMRecord r1, SAMRecord r2) {
	if(r1.getReferenceName()==r2.getReferenceName()){
		if((r1.getAlignmentStart()==r2.getMateAlignmentStart())&& (r1.getMateAlignmentStart()==r2.getAlignmentStart())){
			if(r1.getMateNegativeStrandFlag()==r2.getReadNegativeStrandFlag() && r1.getReadNegativeStrandFlag()==r2.getMateNegativeStrandFlag()){
				if(r1.getMateReferenceName().equalsIgnoreCase(r2.getReferenceName()) && r2.getMateReferenceName().equalsIgnoreCase(r1.getReferenceName())){
					return true;
				}
			}
		}
	}
	return false;
}
 

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * Gets a paired end read.  This method assumes that it will only be passed reads which map
 * to the reverse strand since these are the ones which contain enough information to 
 * unambiguously locate both ends of the pair.
 * 
 * @param sections The tab split sections from the SAM file
 * @param flag The binary flag field
 * @return The read which was read
 * @throws SeqMonkException
 */
private SequenceReadWithChromosome getPairedEndRead (SAMRecord samRecord) throws SeqMonkException {
	int strand;
	int start;
	int end;

	if (!samRecord.getReadNegativeStrandFlag()) {
		throw new SeqMonkException("Read passed to parse pair was not on the negative strand");
	}
	
	if (samRecord.getMateNegativeStrandFlag()) {
		throw new SeqMonkException("Ignored discordantly stranded read pair");
	}
	
	end = samRecord.getAlignmentEnd();
	start = samRecord.getMateAlignmentStart();
	
	if (start > end) {
		throw new SeqMonkException("Ignored discordantly stranded read pair");
	}
	
	if (samRecord.getFirstOfPairFlag()) {
		strand = Location.REVERSE;
	}
	else {
		strand = Location.FORWARD;
	}
	
	if ((end - start)+1 > pairedEndDistance) {
		throw new SeqMonkException("Distance between ends "+((end - start)+1)+" was larger than cutoff ("+pairedEndDistance+")");
	}

	ChromosomeWithOffset c;

	try {
		c = dataCollection().genome().getChromosome(samRecord.getReferenceName());
	}
	catch (Exception e) {
		throw new SeqMonkException(e.getLocalizedMessage());
	}
	
	start = c.position(start);
	end = c.position(end);

	// We also don't allow readings which are beyond the end of the chromosome
	if (end > c.chromosome().length()) {
		int overrun = end - c.chromosome().length();
		throw new SeqMonkException("Reading position "+end+" was "+overrun+"bp beyond the end of chr"+c.chromosome().name()+" ("+c.chromosome().length()+")");
	}
	if (start < 1) {
		throw new SeqMonkException("Reading position "+start+" was before the start of chr"+c.chromosome().name()+" ("+c.chromosome().length()+")");			
	}

	// We can now make the new reading
	SequenceReadWithChromosome read = new SequenceReadWithChromosome(c.chromosome(),SequenceRead.packPosition(start,end,strand));

	return read;
}
 

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
public void processSAMRecord(SAMRecord r){
	totalReads++;
	if(!r.getReadUnmappedFlag()){
		totalHits++;
		int count = 1;  //Have to figure out something for BWA when reporting multiple alignments
		if(r.getIntegerAttribute("NH")!=null)
			count = r.getIntegerAttribute("NH");
		if(count==1 && r.getMappingQuality()!=0) //Second clause for BWA
			uniquelyMapped++;
		
		weight += 1/(float)count;

		if(r.getReadPairedFlag()){
			if(r.getMateUnmappedFlag()){
				singleEnd++;
			}else{
				pairMapped++;
				if(r.getMateReferenceName().equals(r.getReferenceName())){
					pairedEndSameChr++;
				}else{
					pairedEndDiffChr++;
				}
			}
		}else{
			singleEnd++;
		}
		
		List<AlignmentBlock> blocks = r.getAlignmentBlocks();
   		if(blocks.size()>=2){
   			junctions+=blocks.size()-1;
   		}
		
		if(!r.getNotPrimaryAlignmentFlag()){
			if(!r.getReadPairedFlag() || r.getFirstOfPairFlag()){
				LHits++;
				if(r.getReadPairedFlag() && r.getProperPairFlag()){
					properPair++;
					if(!r.getReadNegativeStrandFlag() && r.getMateNegativeStrandFlag()){
						double dist = (r.getMateAlignmentStart()+r.getReadLength())-r.getAlignmentStart();
						histo.addValue(dist);
					}
				}
			}else if(r.getSecondOfPairFlag()){
				RHits++;
			}
		}else{
			notPrimary++;
		}
	}
}
 

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
public void processBT1SAMRecord(SAMRecord r){
	totalReads++;
	if(!r.getReadUnmappedFlag())
		if(r.getIntegerAttribute("XM")!=null){
			int xm = r.getIntegerAttribute("XM");
			if(xm!=0)
				weight++;
		}
	else{
		totalHits++;
		int count =1; //TODO: Fix this if using bowtie for multi-mapping reads
		boolean currUnique = true;
    	
		if(count==1 && currUnique){
			uniquelyMapped++;
		}
		
		weight += 1/(float)count;
		
		if(r.getReadPairedFlag()){
			if(r.getMateUnmappedFlag()){
				singleEnd++;
			}else{
				pairMapped++;
				if(r.getMateReferenceName().equals(r.getReferenceName())){
					pairedEndSameChr++;
				}else{
					pairedEndDiffChr++;
				}
			}
		}else{
			singleEnd++;
		}
		
		List<AlignmentBlock> blocks = r.getAlignmentBlocks();
   		if(blocks.size()>=2){
   			junctions+=blocks.size()-1;
   		}
   		
		if(!r.getNotPrimaryAlignmentFlag()){
			if(!r.getReadPairedFlag() || r.getFirstOfPairFlag()){
				LHits++;
				if(r.getReadPairedFlag() && r.getProperPairFlag()){
					properPair++;
					if(!r.getReadNegativeStrandFlag() && r.getMateNegativeStrandFlag()){
						double dist = (r.getMateAlignmentStart()+r.getReadLength())-r.getAlignmentStart();
						histo.addValue(dist);
					}
				}
			}else if(r.getSecondOfPairFlag()){
				RHits++;
			}
		}else{
			notPrimary++;
		}
	}
}
 

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
public void processBT2SAMRecord(SAMRecord r){
	totalReads++;
	if(!r.getReadUnmappedFlag()){
		totalHits++;
		int count =1; //TODO: Fix this if using bowtie2 for multi-mapping reads
		
		int primAScore = r.getIntegerAttribute("AS");
    	int secAScore=-1000000;
    	if(r.getIntegerAttribute("XS")!=null)
    		secAScore = r.getIntegerAttribute("XS");
    	boolean currUnique = primAScore > secAScore ? true : false;
    	
		if(count==1 && currUnique){
			uniquelyMapped++;
		}
		
		weight += 1/(float)count;
		
		if(r.getReadPairedFlag()){
			if(r.getMateUnmappedFlag()){
				singleEnd++;
			}else{
				pairMapped++;
				if(r.getMateReferenceName().equals(r.getReferenceName())){
					pairedEndSameChr++;
				}else{
					pairedEndDiffChr++;
				}
			}
		}else{
			singleEnd++;
		}
		
		List<AlignmentBlock> blocks = r.getAlignmentBlocks();
   		if(blocks.size()>=2){
   			junctions+=blocks.size()-1;
   		}
   		
		if(!r.getNotPrimaryAlignmentFlag()){
			if(!r.getReadPairedFlag() || r.getFirstOfPairFlag()){
				LHits++;
				if(r.getReadPairedFlag() && r.getProperPairFlag()){
					properPair++;
					if(!r.getReadNegativeStrandFlag() && r.getMateNegativeStrandFlag()){
						double dist = (r.getMateAlignmentStart()+r.getReadLength())-r.getAlignmentStart();
						histo.addValue(dist);
					}
				}
			}else if(r.getSecondOfPairFlag()){
				RHits++;
			}
		}else{
			notPrimary++;
		}
	}
}
 

import net.sf.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * Get sam line for record
 * @param rec The record
 * @return Sam formatted line ending in newline
 */
public static String getSamString(SAMRecord rec) {
	
	// Make string representation of record in sam format
	String rtrn = rec.getReadName() + "\t";
	
	// Get sam flags
	int flags = 0;
	// 0x1 template having multiple segments in sequencing
	if(rec.getReadPairedFlag()) flags += 1;
	// 0x4 segment unmapped
	if(rec.getReadUnmappedFlag()) flags += 4;
	// 0x10 SEQ being reverse complemented
	if(rec.getReadNegativeStrandFlag()) flags += 16;
	// 0x100 secondary alignment
	if(rec.getNotPrimaryAlignmentFlag()) flags += 256;
	// 0x200 not passing quality controls
	if(rec.getReadFailsVendorQualityCheckFlag()) flags += 512;
	// 0x400 PCR or optical duplicate
	if(rec.getDuplicateReadFlag()) flags += 1024;
	
	if(rec.getReadPairedFlag()) {
		// 0x2 each segment properly aligned according to the aligner
		if(rec.getProperPairFlag()) flags += 2;
		// 0x8 next segment in the template unmapped
		if(rec.getMateUnmappedFlag()) flags += 8;
		// 0x20 SEQ of the next segment in the template being reversed
		if(rec.getMateNegativeStrandFlag()) flags += 32;
		// 0x40 the first segment in the template
		if(rec.getFirstOfPairFlag()) flags += 64;
		// 0x80 the last segment in the template
		if(rec.getSecondOfPairFlag()) flags += 128;
	}
	
	rtrn += Integer.valueOf(flags).toString() + "\t";
	
	// The rest of the sam fields
	rtrn += rec.getReferenceName() + "\t";
	rtrn += rec.getAlignmentStart() + "\t";
	rtrn += rec.getMappingQuality() + "\t";
	rtrn += rec.getCigarString() + "\t";
	rtrn += rec.getMateReferenceName() + "\t";
	rtrn += rec.getMateAlignmentStart() + "\t";
	rtrn += rec.getInferredInsertSize() + "\t";
	rtrn += rec.getReadString() + "\t";
	rtrn += rec.getBaseQualityString() + "\n";
	
	return rtrn;
}