jpayne@68: package prok;
jpayne@68: 
jpayne@68: import dna.AminoAcid;
jpayne@68: import shared.Shared;
jpayne@68: import shared.Tools;
jpayne@68: import structures.ByteBuilder;
jpayne@68: 
jpayne@68: /**
jpayne@68:  * ORF means Open Reading Frame.
jpayne@68:  * It starts at the first base of a start codon and ends at the last base of a stop codon.
jpayne@68:  * The length is divisible by 3. 
jpayne@68:  * @author Brian Bushnell
jpayne@68:  * @date Sep 20, 2018
jpayne@68:  *
jpayne@68:  */
jpayne@68: public class Orf extends Feature {
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Initialization        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/** 
jpayne@68: 	 * Bases and coordinates are assumed to be the correct strand.
jpayne@68: 	 * Minus-strand ORFs can be flipped at the end of the constructor.
jpayne@68: 	 * @param scafName_
jpayne@68: 	 * @param start_
jpayne@68: 	 * @param stop_
jpayne@68: 	 * @param strand_
jpayne@68: 	 * @param frame_
jpayne@68: 	 * @param bases
jpayne@68: 	 * @param flip
jpayne@68: 	 */
jpayne@68: 	public Orf(String scafName_, int start_, int stop_, int strand_, int frame_, byte[] bases, boolean flip, int type_) {
jpayne@68: 		super(scafName_, start_, stop_, strand_, bases.length);
jpayne@68: 		frame=frame_;
jpayne@68: 		startCodon=getCodon(start, bases);
jpayne@68: 		stopCodon=getCodon(stop-2, bases);
jpayne@68: 		type=type_;
jpayne@68: 		
jpayne@68: 		if(flip && strand==Shared.MINUS){flip();}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Init Helpers         ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Grab the codon starting at from.
jpayne@68: 	 * Assumes bases are in the correct strand
jpayne@68: 	 * @param from
jpayne@68: 	 * @param bases
jpayne@68: 	 * @return
jpayne@68: 	 */
jpayne@68: 	private static int getCodon(int from, byte[] bases){
jpayne@68: 		int codon=0;
jpayne@68: 		for(int i=0; i<3; i++){
jpayne@68: //			assert(i+from<bases.length) : i+", "+from+", "+bases.length;
jpayne@68: 			byte b=bases[i+from];
jpayne@68: 			int x=AminoAcid.baseToNumber[b];
jpayne@68: 			codon=(codon<<2)|x;
jpayne@68: 		}
jpayne@68: 		return codon;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public float calcOrfScore(){
jpayne@68: 		return calcOrfScore(0);
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	/**
jpayne@68: 	 * The score of an ORF alone is a factor of the length, start score, stop score, and kmer score.
jpayne@68: 	 * The score of an ORF in the context of an overlapping gene also includes a penalty for the overlap length.
jpayne@68: 	 * @param overlap
jpayne@68: 	 * @return Calculated score
jpayne@68: 	 */
jpayne@68: 	public float calcOrfScore(int overlap){
jpayne@68: 		double a=Math.sqrt(Tools.max(f1, e1+startScore));
jpayne@68: //		double b=Math.sqrt(f2/*Tools.max(f2, e2+stopScore)*/);//This is better, ignoring stopscore completely
jpayne@68: 		double b=Math.sqrt(Tools.max(f2, e2+0.35f*stopScore));
jpayne@68: 		double c=Tools.max(f3, e3+averageKmerScore());
jpayne@68: 		assert(a!=Double.NaN);
jpayne@68: 		assert(b!=Double.NaN);
jpayne@68: 		assert(c!=Double.NaN);
jpayne@68: 		c=4*Math.pow(c, 2.2);
jpayne@68: 		double d=(0.1*a*b*c*(Math.pow(length()-overlap, 2.5)-(overlap<1 ? 0 : Math.pow(overlap+50, 2))));//TODO: Adjust these constants
jpayne@68: 		if(d>0){d=Math.sqrt(d);}
jpayne@68: 		assert(d!=Double.NaN);
jpayne@68: 		return (float)d;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float averageKmerScore(){
jpayne@68: 		return kmerScore/(length()-GeneModel.kInnerCDS-2); //This slightly affects score if kInnerCDS is changed
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Public Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public boolean isValidPrev(Orf prev, int maxOverlap){
jpayne@68: 		if(prev.stop>=stop || prev.stop>=start+maxOverlap || prev.start>=start){return false;}
jpayne@68: 		if(prev.frame==frame && prev.strand==strand && prev.stop>=start){return false;}
jpayne@68: 		return true;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public float pathScore() {return Tools.max(pathScorePlus, pathScoreMinus);}
jpayne@68: 	public float pathScore(int prevStrand) {return prevStrand==0 ? pathScorePlus : pathScoreMinus;}
jpayne@68: 
jpayne@68: 	public Orf prev(){return pathScorePlus>=pathScoreMinus ? prevPlus : prevMinus;}
jpayne@68: 	public Orf prev(int prevStrand){return prevStrand==0 ? prevPlus : prevMinus;}
jpayne@68: 
jpayne@68: 	public int pathLength(int prevStrand){return prevStrand==0 ? pathLengthPlus : pathLengthMinus;}
jpayne@68: 	public int pathLength(){return pathScorePlus>=pathScoreMinus ? pathLengthPlus : pathLengthMinus;}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------           ToString           ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public String toStringFlipped(){
jpayne@68: 		if(strand==flipped()){
jpayne@68: 			return toString();
jpayne@68: 		}
jpayne@68: 		flip();
jpayne@68: 		String s=toString();
jpayne@68: 		flip();
jpayne@68: 		return s;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public String toString(){
jpayne@68: 		return toGff();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public String toGff(){
jpayne@68: 		ByteBuilder bb=new ByteBuilder();
jpayne@68: 		appendGff(bb);
jpayne@68: 		return bb.toString();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public ByteBuilder appendGff(ByteBuilder bb){
jpayne@68: 		if(scafName==null){
jpayne@68: 			bb.append('.').tab();
jpayne@68: 		}else{
jpayne@68: 			for(int i=0, max=scafName.length(); i<max; i++){
jpayne@68: 				char c=scafName.charAt(i);
jpayne@68: 				if(c==' ' || c=='\t'){break;}
jpayne@68: 				bb.append(c);
jpayne@68: 			}
jpayne@68: 			bb.tab();
jpayne@68: 		}
jpayne@68: 		bb.append("BBTools").append('\t');
jpayne@68: 		bb.append(typeStrings2[type]).append('\t');
jpayne@68: 		bb.append(start+1).append('\t');
jpayne@68: 		bb.append(stop+1).append('\t');
jpayne@68: 		
jpayne@68: 		if(orfScore<0){bb.append('.').append('\t');}
jpayne@68: 		else{bb.append(orfScore, 2).append('\t');}
jpayne@68: 		
jpayne@68: 		bb.append(strand<0 ? '.' : Shared.strandCodes2[strand]).append('\t');
jpayne@68: 		
jpayne@68: 		bb.append('0').append('\t');
jpayne@68: 
jpayne@68: 		//bb.append('.');
jpayne@68: 		bb.append(typeStrings[type]).append(',');
jpayne@68: 		if(type==0){
jpayne@68: 			bb.append("fr").append(frame).append(',');
jpayne@68: 		}
jpayne@68: //		bb.append(startCodon).append(',');
jpayne@68: //		bb.append(stopCodon).append(',');
jpayne@68: 		bb.append("startScr:").append(startScore, 3).append(',');
jpayne@68: 		bb.append("stopScr:").append(stopScore, 3).append(',');
jpayne@68: 		bb.append("innerScr:").append(averageKmerScore(), 3).append(',');
jpayne@68: 		bb.append("len:").append(length());
jpayne@68: 		if(type==0){
jpayne@68: 			bb.append(',');
jpayne@68: 			bb.append("start:").append(AminoAcid.codonToString(startCodon)).append(',');
jpayne@68: 			bb.append("stop:").append(AminoAcid.codonToString(stopCodon));
jpayne@68: 		}
jpayne@68: 		return bb;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean isSSU(){
jpayne@68: 		return type==r16S || type==r18S;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean is16S(){
jpayne@68: 		return type==r16S;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean is18S(){
jpayne@68: 		return type==r18S;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------          Overrides           ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public float score() {
jpayne@68: 		return orfScore;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Fields            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public final int frame;
jpayne@68: 
jpayne@68: 	//These are not needed but nice for printing
jpayne@68: 	public final int startCodon;
jpayne@68: 	public final int stopCodon;
jpayne@68: 	
jpayne@68: 	public float startScore;
jpayne@68: 	public float stopScore;
jpayne@68: 	public float kmerScore;
jpayne@68: 	
jpayne@68: 	public float orfScore;
jpayne@68: 
jpayne@68: 	//Path scores are for pathfinding phase
jpayne@68: 	
jpayne@68: 	public float pathScorePlus;
jpayne@68: 	public int pathLengthPlus=1;
jpayne@68: 	public Orf prevPlus;
jpayne@68: 	
jpayne@68: 	public float pathScoreMinus;
jpayne@68: 	public int pathLengthMinus=1;
jpayne@68: 	public Orf prevMinus;
jpayne@68: 	
jpayne@68: 	public final int type;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Static Fields        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/* for kinnercds=6 */ 
jpayne@68: //	static float e1=0.1f; 
jpayne@68: //	static float e2=-0.04f; 
jpayne@68: //	static float e3=0.01f;//Decreasing this decreases TP but increases SNR
jpayne@68: //	
jpayne@68: //	static float f1=0.08f; 
jpayne@68: //	static float f2=0.06f; 
jpayne@68: //	static float f3=0.09f;
jpayne@68: 
jpayne@68: 	/* for kinnercds=7 */ 
jpayne@68: 	static float e1=0.35f; 
jpayne@68: 	static float e2=-0.1f; 
jpayne@68: 	static float e3=-0.01f;//Decreasing this decreases TP but increases SNR
jpayne@68: 	
jpayne@68: 	static float f1=0.08f; 
jpayne@68: 	static float f2=0.02f; 
jpayne@68: 	static float f3=0.09f;
jpayne@68: 	
jpayne@68: }