package gem;
import java.util.*;
import gem.util.Progress;
import gem.util.ErrorFunction;
/**
* This class is for calculation difference of propostions and their p-value.
*/
public class Difference implements Constants
{
/**
* Calculates and assigns gamma p-values to the given triplets.
*
* @param trips
*/
public static void assignGammaPval(Collection<Triplet> trips)
{
System.out.print("Assigining significance ...");
Progress p = new Progress(trips.size());
for (Triplet t : trips)
{
t.pvalGamma = calcGammaPval(t.cnt);
// t.pvalGamma = calcModulationPval(t.cnt);
// int[][][] c = CaseCounter.countComplete(t);
// t.dcmi = MINDy.calcDMI(c);
// t.pvalDcmi = MINDy.calcPval(c);
p.tick();
}
}
public static List<Triplet> assignBestCountsConsiderAbsent(Collection<Triplet> trips, Set<Gene> absent)
{
List<Triplet> list = new ArrayList<Triplet>(trips);
for (Triplet t : trips)
{
if (absent.contains(t.M))
{
list.remove(t);
continue;
}
double mod_orig = Difference.calcModulation(t);
// Gene best = null;
// double bestMod = 0;
for (Gene abs : absent)
{
CaseCounter.count(t, abs, ABSENT);
double moda = Difference.calcModulation(t);
CaseCounter.count(t, abs, PRESENT);
double modp = Difference.calcModulation(t);
if (modp * mod_orig > 0 && Math.abs(modp) > Math.abs(mod_orig) && (moda * modp < 0 || Math.abs(moda) < Math.abs(mod_orig) / 2))
{
list.remove(t);
break;
}
else
{
CaseCounter.count(t);
}
// if (mod_orig * mod < 0 && Math.abs(mod) > Math.abs(bestMod))
// {
// bestMod = mod;
// best = abs;
// }
}
// if (best != null)
// CaseCounter.count(t, best, ABSENT);
// else
// list.remove(t);
}
return list;
}
/**
* Calculates and assigns betaM p-values.
*
* @param trips
*/
public static void assignBetaMPval(Collection<Triplet> trips)
{
for (Triplet t : trips)
{
t.pvalBetaM = calcBetaMpval(t.cnt, calcTotals(t.cnt));
}
}
public static double calcGammaPval(Triplet t)
{
return calcGammaPval(t.cnt);
}
public static double calcGammaPval(int[] cnt)
{
int[] n = calcTotals(cnt);
double[] p = calcProportions(cnt, n);
double d = calcGamma(p);
double stdev = calcGammaStdev(cnt);
return 1 - ErrorFunction.getSignif(Math.abs(d) / (stdev * SQRT2));
}
public static double calcGammaStdev(int[] cnt)
{
int[] n = calcTotals(cnt);
double p0123 = (cnt[4] + cnt[5] + cnt[6] + cnt[7]) / (double) (n[0] + n[1] + n[2] + n[3]);
double q0123 = 1 - p0123;
double ovn = (1D / n[0]) + (1D / n[1]) + (1D / n[2]) + (1D / n[3]);
double stdev = Math.sqrt(p0123 * q0123 * ovn);
return stdev;
}
public static double[] calcProportions(int[] cnt, int[] n)
{
double[] p = new double[4];
if (n == null) n = calcTotals(cnt);
for (int i = 0; i < 4; i++)
{
p[i] = cnt[i + 4] / (double) n[i];
}
return p;
}
public static double[] calcProportions(Triplet t)
{
double[] p = new double[4];
int[] n = calcTotals(t.cnt);
for (int i = 0; i < 4; i++)
{
p[i] = t.cnt[i + 4] / (double) n[i];
}
return p;
}
public static double calcGamma(Triplet t)
{
return calcGamma(calcProportions(t.cnt, calcTotals(t.cnt)));
}
public static double calcGamma(int[] cnt)
{
return calcGamma(calcProportions(cnt, calcTotals(cnt)));
}
public static double calcGamma(double[] p)
{
return p[3] + p[0] - p[1] - p[2];
}
public static double calcAlphaF(Triplet t)
{
double[] p = calcProportions(t);
return calcAlphaF(p);
}
public static double calcAlphaF(double[] p)
{
return p[1] - p[0];
}
public static double calcAlphaM(double[] p)
{
return p[2] - p[0];
}
public static double calcAlphaM(Triplet t)
{
double[] p = calcProportions(t);
return p[2] - p[0];
}
public static double calcBetaF(double[] p)
{
return p[3] - p[2];
}
public static double calcBetaF(Triplet t)
{
double[] p = calcProportions(t);
return calcBetaF(p);
}
public static double calcBetaM(double[] p)
{
return p[3] - p[1];
}
public static double calcBetaM(Triplet t)
{
double[] p = calcProportions(t);
return p[3] - p[1];
}
public static double calcAlphaFplusBetaM(double[] p)
{
return p[3] - p[0];
}
public static int[] calcTotals(int[] cnt)
{
int[] n = new int[4];
for (int i = 0; i < 4; i++)
{
n[i] = cnt[i] + cnt[i + 4];
}
return n;
}
public static double calcAlphaFpval(Triplet t)
{
return calcAlphaFpval(t.cnt, calcTotals(t.cnt));
}
public static double calcAlphaFpval(int[] cnt, int[] n)
{
if (n == null) n = calcTotals(cnt);
return calcPairwisePval(cnt[1], n[1], cnt[0], n[0]);
}
public static double calcBetaFpval(int[] cnt, int[] n)
{
if (n == null) n = calcTotals(cnt);
return calcPairwisePval(cnt[3], n[3], cnt[2], n[2]);
}
public static double calcBetaFpval(Triplet t)
{
return calcBetaFpval(t.cnt, calcTotals(t.cnt));
}
public static double calcAlphaMpval(int[] cnt, int[] n)
{
if (n == null) n = calcTotals(cnt);
return calcPairwisePval(cnt[2], n[2], cnt[0], n[0]);
}
public static double calcAlphaMpval(Triplet t)
{
return calcAlphaMpval(t.cnt, calcTotals(t.cnt));
}
public static double calcBetaMpval(Triplet t)
{
return calcBetaMpval(t.cnt, calcTotals(t.cnt));
}
public static double calcBetaMpval(int[] cnt, int[] n)
{
if (n == null) n = calcTotals(cnt);
return calcPairwisePval(cnt[3], n[3], cnt[1], n[1]);
}
public static double calcAlphaFplusBetaMpval(int[] cnt, int[] n)
{
if (n == null) n = calcTotals(cnt);
return calcPairwisePval(cnt[3], n[3], cnt[0], n[0]);
}
public static double calcModulation(Triplet t)
{
return calcModulation(t.cnt);
}
public static double calcModulation(int[] cnt)
{
int[] n = calcTotals(cnt);
double[] p = calcProportions(cnt, n);
double gamma = calcGamma(p);
double alfaM = calcAlphaM(p);
double betaM = calcBetaM(p);
if (betaM * alfaM < 0) return betaM;
if (Math.abs(betaM) > Math.abs(alfaM)) return gamma;
return 0;
}
public static double calcReverseModulation(int[] cnt)
{
int[] n = calcTotals(cnt);
double[] p = calcProportions(cnt, n);
double gamma = calcGamma(p);
double alfaF = calcAlphaF(p);
double betaF = calcBetaF(p);
if (betaF * alfaF < 0) return betaF;
if (Math.abs(betaF) > Math.abs(alfaF)) return gamma;
return 0;
}
public static double calcModulationPval(Triplet t)
{
return calcModulationPval(t.cnt);
}
public static double calcModulationPval(int[] cnt)
{
double mod = calcModulation(cnt);
double stdev = calcGammaStdev(cnt);
return calcPval(mod, stdev);
}
public static double calcOR(int[] cnt)
{
int[] n = calcTotals(cnt);
double[] p = calcProportions(cnt, n);
double alfaF = calcAlphaF(p);
double alfaM = calcAlphaM(p);
double gamma = calcGamma(p);
if (alfaF * alfaM > 0 && gamma * alfaF < 0)
{
if (alfaF < 0) return Math.max(alfaF, alfaM);
else return Math.min(alfaF, alfaM);
}
else
{
return 0;
}
}
public static double calcORPval(int[] cnt)
{
double or = calcOR(cnt);
double stdev = calcGammaStdev(cnt);
return calcPval(or, stdev);
}
public static double calcFMod(int[] cnt)
{
return calcGamma(cnt) - calcModulation(cnt) + calcOR(cnt);
}
public static double calcFModPval(int[] cnt)
{
double fmod = calcFMod(cnt);
double stdev = calcGammaStdev(cnt);
return calcPval(fmod, stdev);
}
/**
* Calculates the p-val of any two proportion differences.
*
* @param c1
* @param n1
* @param c2
* @param n2
* @return p-val
*/
public static double calcPairwisePval(int c1, int n1, int c2, int n2)
{
double d = (c1 / (double) n1) - (c2 / (double) n2);
double stdev = calcStdev(new int[]{c1, c2}, new int[]{n1, n2});
return calcPval(d, stdev);
}
public static double calcStdev(int[] c, int[] n)
{
int tot = 0;
double p = 0;
double ovn = 0;
for (int i = 0; i < c.length; i++)
{
tot += n[i];
p += c[i];
ovn += 1D / n[i];
}
p /= tot;
return Math.sqrt(p * (1 - p) * ovn);
}
public static double calcPval(double dif, double stdev)
{
return 1 - ErrorFunction.getSignif(Math.abs(dif) / (stdev * SQRT2));
}
public static double calcPairwiseCoefficient(Gene x, Gene y)
{
int c[][] = CaseCounter.countTwo(x, y);
double p1 = c[1][1] / (double) (c[1][1] + c[1][0]);
double p0 = c[0][1] / (double) (c[0][1] + c[0][0]);
return p1 - p0;
}
public static double calcPairwiseCoefficientConditionalThird(Gene x, Gene y, Gene z,
boolean positiveCorr)
{
int c[][] = CaseCounter.countTwoWhenThird(x, y, z, !positiveCorr);
double p1 = c[1][1] / (double) (c[1][1] + c[1][0]);
double p0 = c[0][1] / (double) (c[0][1] + c[0][0]);
return p1 - p0;
}
/**
* Checks if A depends B
*/
public static double calcPairwiseDependency(Gene A, Gene B, Gene T)
{
int ca[][] = CaseCounter.countTwo(A, T);
double p1 = ca[1][1] / (double) (ca[1][1] + ca[1][0]);
double p0 = ca[0][1] / (double) (ca[0][1] + ca[0][0]);
double coefA = p1 - p0;
// System.out.println("coefA = " + coefA);
int cb[][] = CaseCounter.countTwo(B, T);
// p1 = cb[1][1] / (double) (cb[1][1] + cb[1][0]);
// p0 = cb[0][1] / (double) (cb[0][1] + cb[0][0]);
//
// double coefB = p1 - p0;
double p = ca[0][1] + ca[1][1] + cb[0][1] + cb[1][1];
double n = p + ca[0][0] + ca[1][0] + cb[0][0] + cb[1][0];
p /= n;
double ovn = (1D / (cb[1][1] + cb[1][0])) + (1D / (cb[0][1] + cb[0][0])) +
(1D / (ca[1][1] + ca[1][0])) + (1D / (ca[0][1] + ca[0][0]));
double std = Math.sqrt(p * (1 - p) * ovn);
boolean corrsign = calcPairwiseCoefficient(A, T) * calcPairwiseCoefficient(B, T) > 0;
double coefAp = calcPairwiseCoefficientConditionalThird(A, T, B, corrsign);
// double coefBp = calcPairwiseCoefficientConditionalThird(B, T, A, corrsign);
// System.out.println("coefAp = " + coefAp);
double difA = coefA - coefAp;
// double difB = coefB - coefBp;
if (Math.abs(coefA) < Math.abs(coefAp)) return 1;
double pvalA = 1 - ErrorFunction.getSignif(Math.abs(difA) / (std * SQRT2));
// double pvalB = 1 - ErrorFunction.getSignif(Math.abs(difB) / (std * SQRT2));
return pvalA;
}
public static double calcPairwiseCoefficientPval(Gene x, Gene y)
{
int c[][] = CaseCounter.countTwo(x, y);
int n0 = c[0][1] + c[0][0];
int n1 = c[1][1] + c[1][0];
double p = (c[1][1] + c[0][1]) / (double) (n0 + n1);
double ovn = (1D / n0) + (1D / n1);
double p1 = c[1][1] / (double) (c[1][1] + c[1][0]);
double p0 = c[0][1] / (double) (c[0][1] + c[0][0]);
double stdev = Math.sqrt(p * (1 - p) * ovn);
return 1 - ErrorFunction.getSignif(Math.abs(p1 - p0) / (stdev * SQRT2));
}
public static boolean complexBetaM(Triplet t, double thr)
{
int[][][] inds = new int[2][][];
double p1 = t.cnt[4] / (double) (t.cnt[0] + t.cnt[4]);
double p2 = t.cnt[5] / (double) (t.cnt[1] + t.cnt[5]);
double p3 = t.cnt[6] / (double) (t.cnt[2] + t.cnt[6]);
double p4 = t.cnt[7] / (double) (t.cnt[3] + t.cnt[7]);
inds[0] = new int[][]{{0, 2},{1, 2},{2, 2}};
if (Math.abs(p2 - p1) > Math.abs(p4 - p3))
{
inds[1] = new int[][]{{0, 0},{0, 1},{0, 2}};
}
else
{
inds[1] = new int[][]{{2, 0},{2, 1},{2, 2}};
}
for (int[][] ind : inds)
{
if (isComplex(t, ind, thr))
{
return true;
}
}
return false;
}
public static boolean isComplex(Triplet t, int[][] inds, double thr)
{
int[][][] cnt = CaseCounter.countComplete(t);
int x0 = cnt[inds[0][0]][inds[0][1]][2];
int x1 = cnt[inds[1][0]][inds[1][1]][2];
int x2 = cnt[inds[2][0]][inds[2][1]][2];
int n0 = cnt[inds[0][0]][inds[0][1]][0] + x0;
int n1 = cnt[inds[1][0]][inds[1][1]][0] + x1;
int n2 = cnt[inds[2][0]][inds[2][1]][0] + x2;
double p0 = x0 / (double) n0;
double p1 = x1 / (double) n1;
double p2 = x2 / (double) n2;
if ((p1 < p0 && p1 < p2) || (p1 > p0 && p1 > p2))
{
if (calcPairwisePval(x1, n1, x0, n0) < thr &&
calcPairwisePval(x1, n1, x2, n2) < thr)
{
return true;
}
}
return false;
}
//----------------------------------------------------------------------------------------------
// Section: Support related
//----------------------------------------------------------------------------------------------
public static double getGammaSupport(int[] cnt, boolean positive)
{
double support = getPositiveGammaSupport(cnt);
if (!positive) support = 1 - support;
return support;
}
public static double getPositiveGammaSupport(int[] cnt)
{
double n = 0;
for (int x : cnt) n += x;
return (cnt[1] + cnt[2] + cnt[4] + cnt[7]) / n;
}
public static boolean tissueSupportIsGreater(Triplet t)
{
double gamma = calcGamma(t);
double support = getGammaSupport(t.cnt, gamma > 0);
double tiss_sup = getGammaSupport(t.cnt_tiss, gamma > 0);
return tiss_sup > support;
}
}