zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Minor cleanup of BAQ calculation (final variables, etc)

This commit is contained in:
Mark DePristo 2012-12-17 16:37:20 -05:00
parent 67fe81391c
commit 66d32f646b
1 changed files with 16 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
public/java/src/org/broadinstitute/sting/utils/baq/BAQ.java
View File

@ -6,6 +6,7 @@ import net.sf.samtools.CigarElement;
import net.sf.samtools.CigarOperator;
import net.sf.samtools.SAMRecord;
import net.sf.samtools.SAMUtils;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
@ -37,6 +38,7 @@ import org.broadinstitute.sting.utils.sam.ReadUtils;
state[i] being wrong.
*/
public class BAQ {
private final static Logger logger = Logger.getLogger(BAQ.class);
private final static boolean DEBUG = false;
public enum CalculationMode {
@ -179,8 +181,7 @@ public class BAQ {
/*** initialization ***/
// change coordinates
int l_ref = ref.length;
final int l_ref = ref.length;
// set band width
int bw2, bw = l_ref > l_query? l_ref : l_query;
@ -266,26 +267,6 @@ public class BAQ {
s[l_query+1] = sum; // the last scaling factor
}
//gdbebug+
/*
double cac=0.;
// undo scaling of forward probabilities to obtain plain probability of observation given model
double[] su = new double[f[l_query].length];
{
double sum = 0.;
double[] logs = new double[s.length];
for (k=0; k < logs.length; k++) {
logs[k] = Math.log10(s[k]);
sum += logs[k];
}
for (k=0; k < f[l_query].length; k++)
su[k]= Math.log10(f[l_query][k])+ sum;
cac = MathUtils.softMax(su);
}
System.out.format("s:%f\n",cac);
// gdebug-
*/
/*** backward ***/
// b[l_query] (b[l_query+1][0]=1 and thus \tilde{b}[][]=1/s[l_query+1]; this is where s[l_query+1] comes from)
for (k = 1; k <= l_ref; ++k) {
@ -305,8 +286,8 @@ public class BAQ {
for (k = end; k >= beg; --k) {
int u, v11, v01, v10;
u = set_u(bw, i, k); v11 = set_u(bw, i+1, k+1); v10 = set_u(bw, i+1, k); v01 = set_u(bw, i, k+1);
double e = (k >= l_ref? 0 : calcEpsilon(ref[k], qyi1, _iqual[qstart+i])) * bi1[v11];
bi[u+0] = e * m[0] + EI * m[1] * bi1[v10+1] + m[2] * bi[v01+2]; // bi1[v11] has been foled into e.
final double e = (k >= l_ref? 0 : calcEpsilon(ref[k], qyi1, _iqual[qstart+i])) * bi1[v11];
bi[u+0] = e * m[0] + EI * m[1] * bi1[v10+1] + m[2] * bi[v01+2]; // bi1[v11] has been folded into e.
bi[u+1] = e * m[3] + EI * m[4] * bi1[v10+1];
bi[u+2] = (e * m[6] + m[8] * bi[v01+2]) * y;
}
@ -332,12 +313,12 @@ public class BAQ {
/*** MAP ***/
for (i = 1; i <= l_query; ++i) {
double sum = 0., max = 0.;
double[] fi = f[i], bi = b[i];
final double[] fi = f[i], bi = b[i];
int beg = 1, end = l_ref, x, max_k = -1;
x = i - bw; beg = beg > x? beg : x;
x = i + bw; end = end < x? end : x;
for (k = beg; k <= end; ++k) {
int u = set_u(bw, i, k);
final int u = set_u(bw, i, k);
double z;
sum += (z = fi[u+0] * bi[u+0]); if (z > max) { max = z; max_k = (k-1)<<2 | 0; }
sum += (z = fi[u+1] * bi[u+1]); if (z > max) { max = z; max_k = (k-1)<<2 | 1; }
@ -531,7 +512,11 @@ public class BAQ {
}
}
// final SimpleTimer total = new SimpleTimer();
// final SimpleTimer local = new SimpleTimer();
// int n = 0;
public BAQCalculationResult calcBAQFromHMM(byte[] ref, byte[] query, byte[] quals, int queryStart, int queryEnd ) {
// total.restart();
if ( queryStart < 0 ) throw new ReviewedStingException("BUG: queryStart < 0: " + queryStart);
if ( queryEnd < 0 ) throw new ReviewedStingException("BUG: queryEnd < 0: " + queryEnd);
if ( queryEnd < queryStart ) throw new ReviewedStingException("BUG: queryStart < queryEnd : " + queryStart + " end =" + queryEnd);
@ -539,7 +524,12 @@ public class BAQ {
// note -- assumes ref is offset from the *CLIPPED* start
BAQCalculationResult baqResult = new BAQCalculationResult(query, quals, ref);
int queryLen = queryEnd - queryStart;
// local.restart();
hmm_glocal(baqResult.refBases, baqResult.readBases, queryStart, queryLen, baqResult.rawQuals, baqResult.state, baqResult.bq);
// local.stop();
// total.stop();
// if ( n++ % 100000 == 0 )
// logger.info("n = " + n + ": Total " + total.getElapsedTimeNano() + " local " + local.getElapsedTimeNano());
return baqResult;
}