jpayne@68: package shared;
jpayne@68: 
jpayne@68: import java.io.Serializable;
jpayne@68: import java.util.ArrayList;
jpayne@68: import java.util.Arrays;
jpayne@68: 
jpayne@68: import align2.QualityTools;
jpayne@68: import dna.AminoAcid;
jpayne@68: import stream.Read;
jpayne@68: import stream.SamLine;
jpayne@68: import stream.SiteScore;
jpayne@68: import structures.ByteBuilder;
jpayne@68: 
jpayne@68: /**
jpayne@68:  * Helper class for processes that do inline quality trimming.
jpayne@68:  * @author Brian Bushnell
jpayne@68:  * @date Mar 15, 2013
jpayne@68:  *
jpayne@68:  */
jpayne@68: public final class TrimRead implements Serializable {
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * 
jpayne@68: 	 */
jpayne@68: 	private static final long serialVersionUID = 8791743639124592480L;
jpayne@68: 	
jpayne@68: 	public static void main(String[] args){
jpayne@68: 		byte[] bases=args[0].getBytes();
jpayne@68: 		byte[] quals=(args.length<2 ? null : args[1].getBytes());
jpayne@68: 		if(quals!=null){
jpayne@68: 			for(int i=0; i<quals.length; i++){quals[i]-=32;}
jpayne@68: 		}
jpayne@68: 		byte[] match=(args.length<3 ? null : args[2].getBytes());
jpayne@68: 		float minq=(args.length<4 ? 5 : Float.parseFloat(args[3]));
jpayne@68: 		float minE=(float)QualityTools.phredToProbError(minq);
jpayne@68: 		Read r=new Read(bases, quals, 1);
jpayne@68: 		r.match=match;
jpayne@68: 		System.out.println("Before trim:\n"+r.toFastq()+(r.match==null ? "" : "\n"+new String(r.match)));
jpayne@68: 		System.out.println(Arrays.toString(r.quality));
jpayne@68: 		TrimRead tr=trim(r, true, true, minq, minE, 1);
jpayne@68: 		System.out.println("\nAfter trim:\n"+r.toFastq()+(r.match==null ? "" : "\n"+new String(r.match)));
jpayne@68: 		if(r.match==null){
jpayne@68: 			r.match=new byte[r.length()];
jpayne@68: 			for(int i=0; i<r.length(); i++){r.match[i]='m';}
jpayne@68: 		}
jpayne@68: 		tr.untrim();
jpayne@68: 		System.out.println("\nAfter untrim:\n"+r.toFastq()+(r.match==null ? "" : "\n"+new String(r.match)));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: //	public static TrimRead trim(Read r, boolean trimLeft, boolean trimRight, float trimq, int minlen){
jpayne@68: //		if(r==null || r.bases==null){return null;}
jpayne@68: //		
jpayne@68: //		final int a, b;
jpayne@68: //		if(optimalMode){
jpayne@68: //			long packed=testOptimal(r.bases, r.quality, QualityTools.PROB_ERROR[trimq]);
jpayne@68: //			a=trimLeft ? (int)((packed>>32)&0xFFFFFFFFL) : 0;
jpayne@68: //			b=trimRight ? (int)((packed)&0xFFFFFFFFL) : 0;
jpayne@68: //		}else if(windowMode){
jpayne@68: //			a=0;
jpayne@68: //			b=(trimRight ? testRightWindow(r.bases, r.quality, (byte)trimq, windowLength) : 0);
jpayne@68: //		}else{
jpayne@68: //			a=(trimLeft ? testLeft(r.bases, r.quality, (byte)trimq) : 0);
jpayne@68: //			b=(trimRight ? testRight(r.bases, r.quality, (byte)trimq) : 0);
jpayne@68: //		}
jpayne@68: //		return (a+b==0 ? null : new TrimRead(r, a, b, trimq, minlen));
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	public static TrimRead trim(Read r, boolean trimLeft, boolean trimRight, 
jpayne@68: 			float trimq, float avgErrorRate, int minlen){
jpayne@68: 		if(r==null || r.bases==null){return null;}
jpayne@68: 		
jpayne@68: 		final int a, b;
jpayne@68: 		if(optimalMode){
jpayne@68: 			long packed=testOptimal(r.bases, r.quality, avgErrorRate);
jpayne@68: 			a=trimLeft ? (int)((packed>>32)&0xFFFFFFFFL) : 0;
jpayne@68: 			b=trimRight ? (int)((packed)&0xFFFFFFFFL) : 0;
jpayne@68: 		}else if(windowMode){
jpayne@68: 			a=0;
jpayne@68: 			b=(trimRight ? testRightWindow(r.bases, r.quality, (byte)trimq, windowLength) : 0);
jpayne@68: 		}else{
jpayne@68: 			a=(trimLeft ? testLeft(r.bases, r.quality, (byte)trimq) : 0);
jpayne@68: 			b=(trimRight ? testRight(r.bases, r.quality, (byte)trimq) : 0);
jpayne@68: 		}
jpayne@68: 		return (a+b==0 ? null : new TrimRead(r, a, b, trimq, minlen));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * @param r Read to trim
jpayne@68: 	 * @param trimLeft Trim left side
jpayne@68: 	 * @param trimRight Trim right side
jpayne@68: 	 * @param trimq Maximum quality to trim
jpayne@68: 	 * @param minResult Ensure trimmed read is at least this long
jpayne@68: 	 * @return Number of bases trimmed
jpayne@68: 	 */
jpayne@68: 	public static int trimFast(Read r, boolean trimLeft, boolean trimRight, 
jpayne@68: 			float trimq, float avgErrorRate, int minResult){
jpayne@68: 		return trimFast(r, trimLeft, trimRight, trimq, avgErrorRate, minResult, false);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * @param r Read to trim
jpayne@68: 	 * @param trimLeft Trim left side
jpayne@68: 	 * @param trimRight Trim right side
jpayne@68: 	 * @param trimq Maximum quality to trim
jpayne@68: 	 * @param minResult Ensure trimmed read is at least this long
jpayne@68: 	 * @return Number of bases trimmed
jpayne@68: 	 */
jpayne@68: 	public static int trimFast(Read r, boolean trimLeft, boolean trimRight, 
jpayne@68: 			float trimq, float avgErrorRate, int minResult, boolean trimClip){
jpayne@68: 		return trimFast(r, trimLeft, trimRight, trimq, avgErrorRate, minResult, 0, trimClip);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * This method allows a "discardUnder" parameter, so that qtrim=r will still discard
jpayne@68: 	 * reads that if left-trimmed, would be shorter than the desired minimum length.
jpayne@68: 	 * It is not presently used.
jpayne@68: 	 * @param r Read to trim
jpayne@68: 	 * @param trimLeft Trim left side
jpayne@68: 	 * @param trimRight Trim right side
jpayne@68: 	 * @param trimq Maximum quality to trim
jpayne@68: 	 * @param minResult Ensure trimmed read is at least this long
jpayne@68: 	 * @param discardUnder Resulting reads shorter than this are not wanted
jpayne@68: 	 * @return Number of bases trimmed
jpayne@68: 	 */
jpayne@68: 	public static int trimFast(Read r, boolean trimLeft, boolean trimRight, 
jpayne@68: 			float trimq, float avgErrorRate, int minResult, int discardUnder, boolean trimClip){
jpayne@68: //		assert(avgErrorRate==(float)QualityTools.phredToProbError(trimq)) : trimq+", "+avgErrorRate+", "+(float)QualityTools.phredToProbError(trimq);
jpayne@68: 		final byte[] bases=r.bases, qual=r.quality;
jpayne@68: 		if(bases==null || bases.length<1){return 0;}
jpayne@68: 		
jpayne@68: //		assert(false) : trimLeft+", "+trimRight+", "+trimq+", "+minResult+", "+discardUnder;
jpayne@68: 		
jpayne@68: 		final int len=r.length();
jpayne@68: 		final int a0, b0;
jpayne@68: 		final int a, b;
jpayne@68: 		if(optimalMode){
jpayne@68: 			long packed=testOptimal(bases, qual, avgErrorRate);
jpayne@68: 			a0=(int)((packed>>32)&0xFFFFFFFFL);
jpayne@68: 			b0=(int)((packed)&0xFFFFFFFFL);
jpayne@68: 			if(trimLeft!=trimRight && discardUnder>0 && len-a0-b0<discardUnder){
jpayne@68: 				a=trimLeft ? len : 0;
jpayne@68: 				b=trimRight ? len : 0;
jpayne@68: 			}else{
jpayne@68: 				a=trimLeft ? a0 : 0;
jpayne@68: 				b=trimRight ? b0 : 0;
jpayne@68: 			}
jpayne@68: //			assert(false) : a0+", "+b0+" -> "+a+", "+b;
jpayne@68: 		}else if(windowMode){
jpayne@68: 			a=0;
jpayne@68: 			b=(trimRight ? testRightWindow(bases, qual, (byte)trimq, windowLength) : 0);
jpayne@68: 		}else{
jpayne@68: 			a=(trimLeft ? testLeft(bases, qual, (byte)trimq) : 0);
jpayne@68: 			b=(trimRight ? testRight(bases, qual, (byte)trimq) : 0);
jpayne@68: 		}
jpayne@68: 		return trimByAmount(r, a, b, minResult, trimClip);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public static boolean untrim(Read r){
jpayne@68: 		if(r==null || r.obj==null){return false;}
jpayne@68: 		if(r.obj.getClass()!=TrimRead.class){return false;}
jpayne@68: 		TrimRead tr=(TrimRead)r.obj;
jpayne@68: 		return tr.untrim();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: //	public TrimRead(Read r_, boolean trimLeft, boolean trimRight, float trimq_, int minlen_){
jpayne@68: //		this(r_, (trimLeft ? testLeft(r_.bases, r_.quality, (byte)trimq_) : 0), (trimRight ? testRight(r_.bases, r_.quality, (byte)trimq_) : 0), trimq_, minlen_);
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	public TrimRead(Read r_, int trimLeft, int trimRight, float trimq_, int minlen_){
jpayne@68: 		minlen_=Tools.max(minlen_, 0);
jpayne@68: 		r=r_;
jpayne@68: 		bases1=r.bases;
jpayne@68: 		qual1=r.quality;
jpayne@68: 		trimq=(byte)trimq_;
jpayne@68: 		assert(bases1!=null || qual1==null) : "\n"+new String(bases1)+"\n"+new String(qual1)+"\n";
jpayne@68: 		assert(bases1==null || qual1==null || bases1.length==qual1.length) : "\n"+new String(bases1)+"\n"+new String(qual1)+"\n";
jpayne@68: 		int trimmed=trim(trimLeft, trimRight, minlen_);
jpayne@68: 		if(trimmed>0){
jpayne@68: 			assert(bases2==null || bases2.length>=minlen_ || bases1.length<minlen_) : bases1.length+", "+bases2.length+", "+minlen_+", "+trimLeft+", "+trimRight;
jpayne@68: 			r.bases=bases2;
jpayne@68: 			r.quality=qual2;
jpayne@68: 			r.obj=this;
jpayne@68: 			trimMatch(r);
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: //	/** Trim the left end of the read, from left to right */
jpayne@68: //	private int trim(final boolean trimLeft, final boolean trimRight, final int minlen){
jpayne@68: //		final int a, b;
jpayne@68: //		if(optimalMode){
jpayne@68: //			long packed=testOptimal(bases1, qual1, QualityTools.PROB_ERROR[trimq]);
jpayne@68: //			a=trimLeft ? (int)((packed>>32)&0xFFFFFFFFL) : 0;
jpayne@68: //			b=trimRight ? (int)((packed)&0xFFFFFFFFL) : 0;
jpayne@68: //		}else{
jpayne@68: //			a=(trimLeft ? testLeft(bases1, qual1, (byte)trimq) : 0);
jpayne@68: //			b=(trimRight ? testRight(bases1, qual1, (byte)trimq) : 0);
jpayne@68: //		}
jpayne@68: //		return trim(a, b, minlen);
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	/** Trim the left end of the read, from left to right */
jpayne@68: 	private int trim(int trimLeft, int trimRight, final int minlen){
jpayne@68: 		assert(trimLeft>=0 && trimRight>=0) : "trimLeft="+trimLeft+", trimRight="+trimRight+", minlen="+minlen+", len="+bases1.length;
jpayne@68: 		assert(trimLeft>0 || trimRight>0) : "trimLeft="+trimLeft+", trimRight="+trimRight+", minlen="+minlen+", len="+bases1.length;
jpayne@68: 		final int maxTrim=Tools.min(bases1.length, bases1.length-minlen);
jpayne@68: 		if(trimLeft+trimRight>maxTrim){
jpayne@68: 			int excess=trimLeft+trimRight-maxTrim;
jpayne@68: 			if(trimLeft>0 && excess>0){
jpayne@68: 				trimLeft=Tools.max(0, trimLeft-excess);
jpayne@68: 				excess=trimLeft+trimRight-maxTrim;
jpayne@68: 			}
jpayne@68: 			if(trimRight>0 && excess>0){
jpayne@68: 				trimRight=Tools.max(0, trimRight-excess);
jpayne@68: 				excess=trimLeft+trimRight-maxTrim;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		leftTrimmed=trimLeft;
jpayne@68: 		rightTrimmed=trimRight;
jpayne@68: 		final int sum=leftTrimmed+rightTrimmed;
jpayne@68: 		
jpayne@68: 		if(verbose){
jpayne@68: 			System.err.println("leftTrimmed="+leftTrimmed+", rightTrimmed="+rightTrimmed+", sum="+sum);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(sum==0){
jpayne@68: 			bases2=bases1;
jpayne@68: 			qual2=qual1;
jpayne@68: 		}else{
jpayne@68: 			bases2=KillSwitch.copyOfRange(bases1, trimLeft, bases1.length-trimRight);
jpayne@68: 			qual2=((qual1==null || (trimLeft+trimRight>=qual1.length)) ? null : KillSwitch.copyOfRange(qual1, trimLeft, qual1.length-trimRight));
jpayne@68: 		}
jpayne@68: 		return sum;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Trim bases outside of leftLoc and rightLoc, excluding leftLoc and rightLoc */
jpayne@68: 	public static int trimToPosition(Read r, int leftLoc, int rightLoc, int minResultingLength){
jpayne@68: 		final int len=r.length();
jpayne@68: 		return trimByAmount(r, leftLoc, len-rightLoc-1, minResultingLength, false);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Remove non-genetic-code from reads */
jpayne@68: 	public static int trimBadSequence(Read r){
jpayne@68: 		final byte[] bases=r.bases, quals=r.quality;
jpayne@68: 		if(bases==null){return 0;}
jpayne@68: 		final int minGenetic=20;
jpayne@68: 		int lastNon=-1;
jpayne@68: 		for(int i=0; i<bases.length; i++){
jpayne@68: 			byte b=bases[i];
jpayne@68: 			if(!AminoAcid.isACGTN(b)){lastNon=i;}
jpayne@68: 			if(i-lastNon>minGenetic){break;}
jpayne@68: 		}
jpayne@68: 		if(lastNon>=0){
jpayne@68: 			r.bases=KillSwitch.copyOfRange(bases, lastNon+1, bases.length);
jpayne@68: 			if(quals!=null){
jpayne@68: 				r.quality=KillSwitch.copyOfRange(quals, lastNon+1, quals.length);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return lastNon+1;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Trim this many bases from each end */
jpayne@68: 	public static int trimByAmount(Read r, int leftTrimAmount, int rightTrimAmount, int minResultingLength){
jpayne@68: 		return trimByAmount(r, leftTrimAmount, rightTrimAmount, minResultingLength, false);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Trim this many bases from each end */
jpayne@68: 	public static int trimByAmount(Read r, int leftTrimAmount, int rightTrimAmount, int minResultingLength, boolean trimClip){
jpayne@68: 
jpayne@68: 		leftTrimAmount=Tools.max(leftTrimAmount, 0);
jpayne@68: 		rightTrimAmount=Tools.max(rightTrimAmount, 0);
jpayne@68: 		
jpayne@68: 		//These assertions are unnecessary if the mapping information will never be used or output.
jpayne@68: //		assert(r.match==null) : "TODO: Handle trimming of reads with match strings.";
jpayne@68: 		assert(r.sites==null) : "TODO: Handle trimming of reads with SiteScores.";
jpayne@68: 		
jpayne@68: 		if(r.match!=null){
jpayne@68: 			return trimReadWithMatch(r, r.samline, leftTrimAmount, rightTrimAmount, minResultingLength, Integer.MAX_VALUE, trimClip);
jpayne@68: 		}else if(r.samline!=null && r.samline.cigar!=null && r.samline.cigarContainsOnlyME()){
jpayne@68: 			return trimReadWithMatchFast(r, r.samline, leftTrimAmount, rightTrimAmount, minResultingLength);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		final byte[] bases=r.bases, qual=r.quality;
jpayne@68: 		final int len=(bases==null ? 0 : bases.length), qlen=(qual==null ? 0 : qual.length);
jpayne@68: 		if(len<1){return 0;}
jpayne@68: 		minResultingLength=Tools.min(len, Tools.max(minResultingLength, 0));
jpayne@68: 		if(leftTrimAmount+rightTrimAmount+minResultingLength>len){
jpayne@68: 			rightTrimAmount=Tools.max(1, len-minResultingLength);
jpayne@68: 			leftTrimAmount=0;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		final int total=leftTrimAmount+rightTrimAmount;
jpayne@68: //		System.err.println("D: L="+leftTrimAmount+", R="+rightTrimAmount+", len="+r.length()+", tot="+total);
jpayne@68: 		if(total>0){
jpayne@68: 			r.bases=KillSwitch.copyOfRange(bases, leftTrimAmount, len-rightTrimAmount);
jpayne@68: 			r.quality=(leftTrimAmount+rightTrimAmount>=qlen ? null : KillSwitch.copyOfRange(qual, leftTrimAmount, qlen-rightTrimAmount));
jpayne@68: 			trimMatch(r);
jpayne@68: 			if(r.stop>r.start){ //TODO:  Fixing mapped coordinates needs more work.
jpayne@68: 				r.start+=leftTrimAmount;
jpayne@68: 				r.stop-=rightTrimAmount;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(verbose){
jpayne@68: 			System.err.println("leftTrimmed="+leftTrimAmount+", rightTrimmed="+rightTrimAmount+
jpayne@68: 					", sum="+total+", final length="+r.length());
jpayne@68: 		}
jpayne@68: 		return total;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Count number of bases that need trimming on each side, and pack into a long */
jpayne@68: 	public static long testOptimal(byte[] bases, byte[] qual, float avgErrorRate){
jpayne@68: 		if(optimalBias>=0){avgErrorRate=optimalBias;}//Override
jpayne@68: 		assert(avgErrorRate>0 && avgErrorRate<=1) : "Average error rate ("+avgErrorRate+") must be between 0 (exclusive) and 1 (inclusive)";
jpayne@68: 		if(bases==null || bases.length==0){return 0;}
jpayne@68: 		if(qual==null){return avgErrorRate>=1 ? 0 : ((((long)testLeftN(bases))<<32) | (((long)testRightN(bases))&0xFFFFFFFFL));}
jpayne@68: 		
jpayne@68: 		float maxScore=0;
jpayne@68: 		float score=0;
jpayne@68: 		int maxLoc=-1;
jpayne@68: 		int maxCount=-1;
jpayne@68: 		int count=0;
jpayne@68: 		
jpayne@68: 		final float nprob=Tools.max(Tools.min(avgErrorRate*1.1f, 1), NPROB);
jpayne@68: 		
jpayne@68: 		for(int i=0; i<bases.length; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			byte q=qual[i];
jpayne@68: //			float probError=(b=='N' ? nprob : ADJUST_QUALITY ? CalcTrueQuality.estimateErrorProbAvg(qual, bases, i) : QualityTools.PROB_ERROR[q]);
jpayne@68: //			float probError=(b=='N' ? nprob : ADJUST_QUALITY ? CalcTrueQuality.estimateErrorProbGeoAvg(qual, bases, i) : QualityTools.PROB_ERROR[q]);
jpayne@68: //			float probError=(b=='N' ? nprob : ADJUST_QUALITY ? CalcTrueQuality.estimateErrorProb2(qual, bases, i) : QualityTools.PROB_ERROR[q]);
jpayne@68: 
jpayne@68: //			float probError=(b=='N' ? nprob : q==1 ? PROB1 : QualityTools.PROB_ERROR[q]);
jpayne@68: //			float probError=(b=='N' ? nprob : QualityTools.PROB_ERROR[q]);
jpayne@68: 			float probError=((b=='N' || q<1) ? nprob : QualityTools.PROB_ERROR[q]);
jpayne@68: 			
jpayne@68: //			assert(q>0 || b=='N') : "index "+i+": q="+q+", b="+(char)b+"\n"+new String(bases)+"\n"+Arrays.toString(qual)+"\n";
jpayne@68: 			
jpayne@68: 			float delta=avgErrorRate-probError;
jpayne@68: 			score=score+delta;
jpayne@68: 			if(score>0){
jpayne@68: 				count++;
jpayne@68: 				if(score>maxScore || (score==maxScore && count>maxCount)){
jpayne@68: 					maxScore=score;
jpayne@68: 					maxCount=count;
jpayne@68: 					maxLoc=i;
jpayne@68: 				}
jpayne@68: 			}else{
jpayne@68: 				score=0;
jpayne@68: 				count=0;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		final int left, right;
jpayne@68: 		if(maxScore>0){
jpayne@68: 			assert(maxLoc>=0);
jpayne@68: 			assert(maxCount>0);
jpayne@68: 			left=maxLoc-maxCount+1;
jpayne@68: 			assert(left>=0 && left<=bases.length);
jpayne@68: 			right=bases.length-maxLoc-1;
jpayne@68: 		}else{
jpayne@68: 			left=0;
jpayne@68: 			right=bases.length;
jpayne@68: 		}
jpayne@68: 		final long packed=((((long)left)<<32) | (((long)right)&0xFFFFFFFFL));
jpayne@68: 		
jpayne@68: 		if(verbose){
jpayne@68: 			System.err.println(Arrays.toString(qual));
jpayne@68: 			System.err.println("After testLocal: maxScore="+maxScore+", maxLoc="+maxLoc+", maxCount="+maxCount+
jpayne@68: 					", left="+left+", right="+right+", returning "+Long.toHexString(packed));
jpayne@68: 		}
jpayne@68: 		return packed;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Count number of bases that need trimming on left side */
jpayne@68: 	private static int testLeft(byte[] bases, byte[] qual, final byte trimq){
jpayne@68: 		if(bases==null || bases.length==0){return 0;}
jpayne@68: 		if(qual==null){return trimq<0 ? 0 : testLeftN(bases);}
jpayne@68: 		int good=0;
jpayne@68: 		int lastBad=-1;
jpayne@68: 		int i=0;
jpayne@68: 		for(; i<bases.length && good<minGoodInterval; i++){
jpayne@68: 			final byte q=qual[i];
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			assert(q>0 || b=='N') : "index "+i+": q="+q+", b="+(char)b+"\n"+new String(bases)+"\n"+Arrays.toString(qual)+"\n";
jpayne@68: 			if(q>trimq){good++;}
jpayne@68: 			else{good=0; lastBad=i;}
jpayne@68: 		}
jpayne@68: 		if(verbose){
jpayne@68: //			System.err.println(Arrays.toString(qual));
jpayne@68: 			System.err.println("After testLeft: good="+good+", lastBad="+lastBad+", i="+i+", returning "+(lastBad+1));
jpayne@68: //			assert(false);
jpayne@68: 		}
jpayne@68: 		return lastBad+1;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Count number of bases that need trimming on right side using a sliding window */
jpayne@68: 	private static int testRightWindow(byte[] bases, byte[] qual, final byte trimq, final int window){
jpayne@68: 		if(bases==null || bases.length==0){return 0;}
jpayne@68: 		if(qual==null || qual.length<window){return trimq>0 ? 0 : testRightN(bases);}
jpayne@68: 		final int thresh=Tools.max(window*trimq, 1);
jpayne@68: 		int sum=0;
jpayne@68: 		for(int i=0, j=-window; i<qual.length; i++, j++){
jpayne@68: 			final byte q=qual[i];
jpayne@68: 			sum+=q;
jpayne@68: 			if(j>=-1){
jpayne@68: 				if(j>=0){sum-=qual[j];}
jpayne@68: 				if(sum<thresh){
jpayne@68: 					return qual.length-j-1;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return 0;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Count number of bases that need trimming on right side */
jpayne@68: 	private static int testRight(byte[] bases, byte[] qual, final byte trimq){
jpayne@68: 		if(bases==null || bases.length==0){return 0;}
jpayne@68: 		if(qual==null){return trimq<0 ? 0 : testRightN(bases);}
jpayne@68: 		int good=0;
jpayne@68: 		int lastBad=bases.length;
jpayne@68: 		int i=bases.length-1;
jpayne@68: 		for(; i>=0 && good<minGoodInterval; i--){
jpayne@68: 			final byte q=qual[i];
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			assert(q>0 || b=='N') : "index "+i+": q="+q+", b="+(char)b+"\n"+new String(bases)+"\n"+Arrays.toString(qual)+"\n";
jpayne@68: 			if(q>trimq){good++;}
jpayne@68: 			else{good=0; lastBad=i;}
jpayne@68: 		}
jpayne@68: 		if(verbose){
jpayne@68: 			System.err.println("After trimLeft: good="+good+", lastBad="+lastBad+", i="+i+", returning "+(bases.length-lastBad));
jpayne@68: 		}
jpayne@68: 		return bases.length-lastBad;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Count number of bases that need trimming on left side, considering only N as bad */
jpayne@68: 	public static int testLeftN(byte[] bases){
jpayne@68: 		if(bases==null || bases.length==0){return 0;}
jpayne@68: 		int good=0;
jpayne@68: 		int lastBad=-1;
jpayne@68: 		for(int i=0; i<bases.length && good<minGoodInterval; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			//if(dna.AminoAcid.isFullyDefined(b)){good++;}
jpayne@68: 			if(b!=((byte)'N')){good++;}
jpayne@68: 			else{good=0; lastBad=i;}
jpayne@68: 		}
jpayne@68: 		return lastBad+1;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Count number of bases that need trimming on right side, considering only N as bad */
jpayne@68: 	public static int testRightN(byte[] bases){
jpayne@68: 		if(bases==null || bases.length==0){return 0;}
jpayne@68: 		int good=0;
jpayne@68: 		int lastBad=bases.length;
jpayne@68: 		for(int i=bases.length-1; i>=0 && good<minGoodInterval; i--){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			//if(dna.AminoAcid.isFullyDefined(b)){good++;}
jpayne@68: 			if(b!=((byte)'N')){good++;}
jpayne@68: 			else{good=0; lastBad=i;}
jpayne@68: 		}
jpayne@68: 		return bases.length-lastBad;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean untrim(){
jpayne@68: 		if(leftTrimmed==0 && rightTrimmed==0){return false;}
jpayne@68: 		r.setPerfect(false);
jpayne@68: 		
jpayne@68: 		final int lt, rt;
jpayne@68: 		if(r.strand()==Shared.PLUS){
jpayne@68: 			lt=leftTrimmed;
jpayne@68: 			rt=rightTrimmed;
jpayne@68: 		}else{
jpayne@68: 			lt=rightTrimmed;
jpayne@68: 			rt=leftTrimmed;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		boolean returnToShort=false;
jpayne@68: 		if(verbose){System.err.println("Untrimming");}
jpayne@68: 		if(r.match!=null){
jpayne@68: 			if(r.shortmatch()){
jpayne@68: 				r.toLongMatchString(false);
jpayne@68: 				returnToShort=true;
jpayne@68: 			}
jpayne@68: 			byte[] match2=new byte[r.match.length+lt+rt];
jpayne@68: 			int i=0;
jpayne@68: 			for(; i<lt; i++){
jpayne@68: 				match2[i]='C';
jpayne@68: 			}
jpayne@68: 			for(int j=0; j<r.match.length; i++, j++){
jpayne@68: 				match2[i]=r.match[j];
jpayne@68: 			}
jpayne@68: 			for(; i<match2.length; i++){
jpayne@68: 				match2[i]='C';
jpayne@68: 			}
jpayne@68: 			r.match=match2;
jpayne@68: 		}
jpayne@68: 		r.bases=bases1;
jpayne@68: 		r.quality=qual1;
jpayne@68: 		r.start-=lt;
jpayne@68: 		r.stop+=rt;
jpayne@68: 		if(returnToShort){
jpayne@68: 			r.toShortMatchString(true);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(r.sites!=null){
jpayne@68: 			for(SiteScore ss : r.sites){untrim(ss);}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		return true;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private boolean untrim(SiteScore ss){
jpayne@68: 		if(ss==null){return false;}
jpayne@68: 		if(leftTrimmed==0 && rightTrimmed==0){return false;}
jpayne@68: 		ss.perfect=ss.semiperfect=false;
jpayne@68: 		
jpayne@68: 		final int lt, rt;
jpayne@68: 		if(ss.strand==Shared.PLUS){
jpayne@68: 			lt=leftTrimmed;
jpayne@68: 			rt=rightTrimmed;
jpayne@68: 		}else{
jpayne@68: 			lt=rightTrimmed;
jpayne@68: 			rt=leftTrimmed;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		boolean returnToShort=false;
jpayne@68: 		if(verbose){System.err.println("Untrimming ss "+ss);}
jpayne@68: 		if(ss.match!=null){
jpayne@68: 			
jpayne@68: 			boolean shortmatch=false;
jpayne@68: 			for(byte b : ss.match){
jpayne@68: 				if(Tools.isDigit(b)){shortmatch=true; break;}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(shortmatch){
jpayne@68: 				ss.match=Read.toLongMatchString(ss.match);
jpayne@68: 				returnToShort=true;
jpayne@68: 			}
jpayne@68: 			byte[] match2=new byte[ss.match.length+lt+rt];
jpayne@68: 			int i=0;
jpayne@68: 			for(; i<lt; i++){
jpayne@68: 				match2[i]='C';
jpayne@68: 			}
jpayne@68: 			for(int j=0; j<ss.match.length; i++, j++){
jpayne@68: 				match2[i]=ss.match[j];
jpayne@68: 			}
jpayne@68: 			for(; i<match2.length; i++){
jpayne@68: 				match2[i]='C';
jpayne@68: 			}
jpayne@68: 			ss.match=match2;
jpayne@68: 		}
jpayne@68: 		ss.setLimits(ss.start-lt, ss.stop+rt);
jpayne@68: 		if(returnToShort){ss.match=Read.toShortMatchString(ss.match);}
jpayne@68: 		return true;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private static boolean trimMatch(Read r){
jpayne@68: 		if(r.match==null && r.sites==null){return false;}
jpayne@68: 		
jpayne@68: 		//Easy mode!
jpayne@68: 		r.match=null;
jpayne@68: 		if(r.sites!=null){
jpayne@68: 			for(SiteScore ss : r.sites){
jpayne@68: 				if(ss!=null){ss.match=null;}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return true;
jpayne@68: 		
jpayne@68: 		//TODO - need to adjust read start and stop based on this information.  Also check strand!
jpayne@68: //		byte[] match=r.match;
jpayne@68: //		if(r.shortmatch()){match=Read.toLongMatchString(match);}
jpayne@68: //		byte[] match2=new byte[match.length-leftTrimmed-rightTrimmed];
jpayne@68: //		for(int mpos=0, bpos=0; bpos<leftTrimmed; mpos++){
jpayne@68: //			byte m=match[mpos];
jpayne@68: //			if(m=='D'){
jpayne@68: //				//do nothing
jpayne@68: //			}else{
jpayne@68: //				bpos++;
jpayne@68: //			}
jpayne@68: //		}
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	/** Special case of 100% match, or no match string */
jpayne@68: 	public static int trimReadWithMatchFast(final Read r, final SamLine sl, final int leftTrimAmount, final int rightTrimAmount, final int minFinalLength){
jpayne@68: 		assert(r.match!=null || sl!=null);
jpayne@68: 		if(r.bases==null){return 0;}
jpayne@68: 		if(leftTrimAmount<1 && rightTrimAmount<1){return 0;}
jpayne@68: 		if(leftTrimAmount+rightTrimAmount>=r.length()){return -leftTrimAmount-rightTrimAmount;}
jpayne@68: 		
jpayne@68: 		final boolean shortmatch=r.shortmatch();
jpayne@68: 		final byte[] oldMatch=r.match;
jpayne@68: 		r.match=null;
jpayne@68: 		r.samline=null;
jpayne@68: 		final int trimmed;
jpayne@68: //		System.err.println(rightTrimAmount+", "+leftTrimAmount);
jpayne@68: 		if(sl!=null && sl.strand()==Shared.MINUS){
jpayne@68: 			trimmed=trimByAmount(r, rightTrimAmount, leftTrimAmount, minFinalLength, false);
jpayne@68: 		}else{
jpayne@68: 			trimmed=trimByAmount(r, leftTrimAmount, rightTrimAmount, minFinalLength, false);
jpayne@68: 		}
jpayne@68: 		r.samline=sl;
jpayne@68: 		if(trimmed<1){
jpayne@68: 			r.match=oldMatch;
jpayne@68: 			return 0;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		final int len=r.length();
jpayne@68: 		ByteBuilder bb=new ByteBuilder();
jpayne@68: 		if(oldMatch!=null){
jpayne@68: 			if(shortmatch){
jpayne@68: 				bb.append((byte)'m');
jpayne@68: 				if(len>1){bb.append(len);}
jpayne@68: 			}else{
jpayne@68: 				for(int i=0; i<len; i++){bb.append((byte)'m');}
jpayne@68: 			}
jpayne@68: 			r.match=bb.toBytes();
jpayne@68: 			bb.clear();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(sl!=null){
jpayne@68: 			sl.pos+=leftTrimAmount;
jpayne@68: 			if(sl.cigar!=null && sl.cigar.length()>0){
jpayne@68: 				char c=sl.cigar.charAt(sl.cigar.length()-1);
jpayne@68: 				assert(c=='M' || c=='=') : c+"; "+sl.cigar+"\n"+sl;
jpayne@68: 				bb.append(len);
jpayne@68: 				bb.append(SamLine.VERSION>1.3 ? '=' : 'M');
jpayne@68: 				sl.cigar=bb.toString();
jpayne@68: 			}
jpayne@68: 			sl.seq=r.bases;
jpayne@68: 			sl.qual=r.quality;
jpayne@68: 			if(trimmed>0 && sl.optional!=null && sl.optional.size()>0){
jpayne@68: 				ArrayList<String> list=new ArrayList<String>(2);
jpayne@68: 				for(int i=0; i<sl.optional.size(); i++){
jpayne@68: 					String s=sl.optional.get(i);
jpayne@68: 					if(s.startsWith("PG:") || s.startsWith("RG:") || s.startsWith("X") || s.startsWith("Y") || s.startsWith("Z")){list.add(s);} //Only keep safe flags.
jpayne@68: 				}
jpayne@68: 				sl.optional.clear();
jpayne@68: 				sl.optional.addAll(list);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return trimmed;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: //	//Should be unneeded, just use the above function
jpayne@68: //	public static int trimReadWithoutMatch(final Read r, final SamLine sl, final int leftTrimAmount, final int rightTrimAmount, final int minFinalLength){
jpayne@68: //		if(r.bases==null){return 0;}
jpayne@68: //		if(leftTrimAmount<1 && rightTrimAmount<1){return 0;}
jpayne@68: //		if(leftTrimAmount+rightTrimAmount>=r.length()){return -leftTrimAmount-rightTrimAmount;}
jpayne@68: //		
jpayne@68: //		assert(r.match==null);
jpayne@68: //		final int trimmed;
jpayne@68: //		if(sl!=null && sl.strand()==Shared.MINUS){
jpayne@68: //			trimmed=trimByAmount(r, rightTrimAmount, leftTrimAmount, minFinalLength, false);
jpayne@68: //		}else{
jpayne@68: //			trimmed=trimByAmount(r, leftTrimAmount, rightTrimAmount, minFinalLength, false);
jpayne@68: //		}
jpayne@68: //		if(trimmed<1){return 0;}
jpayne@68: //		
jpayne@68: //		if(sl!=null){
jpayne@68: //			sl.pos+=leftTrimAmount;
jpayne@68: //			assert(sl.cigar==null);
jpayne@68: //			sl.seq=r.bases;
jpayne@68: //			sl.qual=r.quality;
jpayne@68: //			if(trimmed>0 && sl.optional!=null && sl.optional.size()>0){
jpayne@68: //				ArrayList<String> list=new ArrayList<String>(2);
jpayne@68: //				for(int i=0; i<sl.optional.size(); i++){
jpayne@68: //					String s=sl.optional.get(i);
jpayne@68: //					if(s.startsWith("PG:") || s.startsWith("RG:") || s.startsWith("X") || s.startsWith("Y") || s.startsWith("Z")){list.add(s);} //Only keep safe flags.
jpayne@68: //				}
jpayne@68: //				sl.optional.clear();
jpayne@68: //				sl.optional.addAll(list);
jpayne@68: //			}
jpayne@68: //		}
jpayne@68: //		return trimmed;
jpayne@68: //	}
jpayne@68: 
jpayne@68: 	//TODO: This is slow
jpayne@68: 	//TODO: Note, returns a negative number if the whole read is supposed to be trimmed
jpayne@68: 	public static int trimReadWithMatch(final Read r, final SamLine sl, 
jpayne@68: 			int leftTrimAmount, int rightTrimAmount, int minFinalLength, long scafLen, boolean trimClip){
jpayne@68: 		if(r.bases==null || (leftTrimAmount<1 && rightTrimAmount<1 && !trimClip)){return 0;}
jpayne@68: 		if(!r.containsNonM() && !trimClip){//TODO: H not handled
jpayne@68: 			return trimReadWithMatchFast(r, sl, leftTrimAmount, rightTrimAmount, minFinalLength);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: //		if(r.match==null){
jpayne@68: //			assert(sl.cigar==null);
jpayne@68: //			int trimmed=trimByAmount(r, leftTrimAmount, rightTrimAmount, minFinalLength, false);
jpayne@68: //			if(sl!=null){
jpayne@68: //				assert(sl.cigar==null);
jpayne@68: //				int start=sl.pos-1;
jpayne@68: //				int stop=start+sl.seq.length;
jpayne@68: //				sl.seq=r.bases;
jpayne@68: //				sl.qual=r.quality;
jpayne@68: //				if(trimmed>0 && sl.optional!=null && sl.optional.size()>0){
jpayne@68: //					ArrayList<String> list=new ArrayList<String>(2);
jpayne@68: //					for(int i=0; i<sl.optional.size(); i++){
jpayne@68: //						String s=sl.optional.get(i);
jpayne@68: //						if(s.startsWith("PG:") || s.startsWith("RG:") || s.startsWith("X") || s.startsWith("Y") || s.startsWith("Z")){list.add(s);} //Only keep safe flags.
jpayne@68: //					}
jpayne@68: //					sl.optional.clear();
jpayne@68: //					sl.optional.addAll(list);
jpayne@68: //				}
jpayne@68: //			}
jpayne@68: //		}
jpayne@68: 		
jpayne@68: 		if(trimClip){
jpayne@68: 			r.toLongMatchString(false);
jpayne@68: 			byte[] match=r.match;
jpayne@68: 			int leftClip=0, rightClip=0;
jpayne@68: 			for(int i=0; i<match.length; i++){
jpayne@68: 				if(match[i]=='C'){leftClip++;}
jpayne@68: 				else{break;}
jpayne@68: 			}
jpayne@68: 			for(int i=match.length-1; i>=0; i--){
jpayne@68: 				if(match[i]=='C'){rightClip++;}
jpayne@68: 				else{break;}
jpayne@68: 			}
jpayne@68: 			leftTrimAmount=Tools.max(leftTrimAmount, leftClip);
jpayne@68: 			rightTrimAmount=Tools.max(rightTrimAmount, rightClip);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(leftTrimAmount<1 && rightTrimAmount<1){return 0;}
jpayne@68: 		if(leftTrimAmount+rightTrimAmount>=r.length()){return -leftTrimAmount-rightTrimAmount;}
jpayne@68: 		
jpayne@68: 		final int oldPos=(sl==null ? 1 : sl.pos);
jpayne@68: 		
jpayne@68: 		r.toLongMatchString(false);
jpayne@68: 		byte[] match=r.match;
jpayne@68: 		
jpayne@68: 		assert(!r.shortmatch());
jpayne@68: 
jpayne@68: //		System.err.println("Q: "+sl);
jpayne@68: //		System.err.println(new String(match));
jpayne@68: 		
jpayne@68: 		ByteBuilder bb=new ByteBuilder(match.length);
jpayne@68: //		System.err.println("leftTrimAmount="+leftTrimAmount+", rightTrimAmount="+rightTrimAmount);
jpayne@68: 		
jpayne@68: 		int leftTrim=leftTrimAmount>0 ? r.length() : 0, rightTrim=rightTrimAmount>0 ? r.length() : 0;
jpayne@68: //		System.err.println("leftTrim="+leftTrim+", rightTrim="+rightTrim);
jpayne@68: 		{
jpayne@68: 			final int mlen=match.length;
jpayne@68: 			boolean keep=leftTrimAmount<1;
jpayne@68: 			int rpos=(sl==null ? 1 : sl.pos);
jpayne@68: 			for(int mpos=0, cpos=0; mpos<mlen; mpos++){
jpayne@68: 				final byte m=match[mpos];
jpayne@68: 				if(m=='m' || m=='S' || m=='V' || m=='N'){
jpayne@68: 					cpos++;
jpayne@68: 					rpos++;
jpayne@68: 				}else if(m=='D'){
jpayne@68: 					rpos++;
jpayne@68: 				}else if(m=='I' || m=='C'){
jpayne@68: 					cpos++;
jpayne@68: 				}else if(m=='d'){
jpayne@68: 					rpos++;
jpayne@68: 				}else if(m=='i'){
jpayne@68: 					cpos++;
jpayne@68: 				}else{
jpayne@68: 					assert(false) : "Unknown symbol "+(char)m;
jpayne@68: 				}
jpayne@68: 
jpayne@68: 				if(keep){
jpayne@68: 					bb.append(m);
jpayne@68: 				}else if(cpos>=leftTrimAmount){
jpayne@68: 					byte next=(mpos<mlen-1 ? match[mpos+1] : (byte)'m');
jpayne@68: 					if(m=='m' || m=='S' ||  m=='V' ||  m=='N' || next=='m' || next=='S' || next=='V' || next=='N'){
jpayne@68: 						keep=true;
jpayne@68: 						if(sl!=null){sl.pos=rpos;}
jpayne@68: 						leftTrim=cpos;
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: //		System.err.println("R: "+sl);
jpayne@68: 		match=bb.toBytes();
jpayne@68: 		Tools.reverseInPlace(match);
jpayne@68: 		bb.clear();
jpayne@68: 		{
jpayne@68: 			final int mlen=match.length;
jpayne@68: 			boolean keep=rightTrimAmount<1;
jpayne@68: 			for(int mpos=0, cpos=0; mpos<mlen; mpos++){
jpayne@68: 				final byte m=match[mpos];
jpayne@68: 				if(m=='m' || m=='S' || m=='V' || m=='N' || m=='I' || m=='C'){
jpayne@68: 					cpos++;
jpayne@68: 				}else if(m=='D'){
jpayne@68: 				}else if(m=='d'){
jpayne@68: 				}else if(m=='i'){
jpayne@68: 					cpos++;
jpayne@68: 				}else{
jpayne@68: 					assert(false) : "Unknown symbol "+(char)m;
jpayne@68: 				}
jpayne@68: 
jpayne@68: 				if(keep){
jpayne@68: 					bb.append(m);
jpayne@68: 				}else if(cpos>=rightTrimAmount){
jpayne@68: 					byte next=(mpos<mlen-1 ? match[mpos+1] : (byte)'m');
jpayne@68: //					if(m=='m' || m=='S' ||  m=='V' ||  m=='N' || next=='m' || next=='S' || next=='V' || next=='N'){
jpayne@68: 					if(next!='I' && next!='D' && next!='i' && next!='d'){
jpayne@68: 						keep=true;
jpayne@68: 						rightTrim=cpos;
jpayne@68: 					}
jpayne@68: //					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: //		System.err.println("S: "+sl);
jpayne@68: 		match=bb.toBytes();
jpayne@68: 		Tools.reverseInPlace(match);
jpayne@68: 		
jpayne@68: //		System.err.println("leftTrim="+leftTrim+", rightTrim="+rightTrim);
jpayne@68: 		
jpayne@68: 		if(leftTrim+rightTrim>=r.length()){
jpayne@68: 			if(sl!=null){sl.pos=oldPos;}
jpayne@68: 			return -leftTrim-rightTrim;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: //		System.err.println("T: "+sl);
jpayne@68: 		r.match=null;
jpayne@68: 		if(sl!=null && sl.strand()==Shared.MINUS){
jpayne@68: 			int temp=leftTrim;
jpayne@68: 			leftTrim=rightTrim;
jpayne@68: 			rightTrim=temp;
jpayne@68: 		}
jpayne@68: 		final int trimmed=trimByAmount(r, leftTrim, rightTrim, minFinalLength, false);
jpayne@68: //		System.err.println("tba: "+leftTrim+", "+rightTrim+", "+minFinalLength);
jpayne@68: 		r.match=match;
jpayne@68: 		
jpayne@68: 		if(sl!=null){
jpayne@68: 			int start=sl.pos-1;
jpayne@68: 			int stop=start+Read.calcMatchLength(match)-1;
jpayne@68: 			if(SamLine.VERSION>1.3){
jpayne@68: 				sl.cigar=SamLine.toCigar14(match, start, stop, scafLen, r.bases);
jpayne@68: 			}else{
jpayne@68: 				sl.cigar=SamLine.toCigar13(match, start, stop, scafLen, r.bases);
jpayne@68: 			}
jpayne@68: 			sl.seq=r.bases;
jpayne@68: 			sl.qual=r.quality;
jpayne@68: 			if(trimmed>0 && sl.optional!=null && sl.optional.size()>0){
jpayne@68: 				ArrayList<String> list=new ArrayList<String>(2);
jpayne@68: 				for(int i=0; i<sl.optional.size(); i++){
jpayne@68: 					String s=sl.optional.get(i);
jpayne@68: 					if(s.startsWith("PG:") || s.startsWith("RG:") || s.startsWith("X") || s.startsWith("Y") || s.startsWith("Z")){list.add(s);} //Only keep safe flags.
jpayne@68: 				}
jpayne@68: 				sl.optional.clear();
jpayne@68: 				sl.optional.addAll(list);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return trimmed;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public int trimmed() {
jpayne@68: 		return leftTrimmed+rightTrimmed;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public final Read r;
jpayne@68: 	
jpayne@68: 	/** untrimmed bases */
jpayne@68: 	public final byte[] bases1;
jpayne@68: 	/** untrimmed qualities */
jpayne@68: 	public final byte[] qual1;
jpayne@68: 	/** trimmed bases */
jpayne@68: 	public byte[] bases2;
jpayne@68: 	/** trimmed qualities */
jpayne@68: 	public byte[] qual2;
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	public final float trimq;
jpayne@68: 	public int leftTrimmed;
jpayne@68: 	public int rightTrimmed;
jpayne@68: 	
jpayne@68: 	/** Require this many consecutive good bases to stop trimming.  Minimum is 1.
jpayne@68: 	 * This is for the old trimming mode and not really used anymore */
jpayne@68: 	public static int minGoodInterval=2;
jpayne@68: 	
jpayne@68: 	public static boolean verbose=false;
jpayne@68: 	public static boolean optimalMode=true;
jpayne@68: 	public static boolean windowMode=false;
jpayne@68: 	public static float optimalBias=-1f;
jpayne@68: 	
jpayne@68: 	public static int windowLength=4;
jpayne@68: 	
jpayne@68: 	private static final float NPROB=0.75f;
jpayne@68: //	public static float PROB1=QualityTools.PROB_ERROR[1];
jpayne@68: 	
jpayne@68: 	
jpayne@68: }