R/packageAW.R

In metaOmics/MetaPath: Perform the Meta-Analysis for Pathway Enrichment

####### August 1, 2015 #######
#######   Added aw weights ###

####################################
### Main Function ##################

###### COLUMN --> samples
###### ROW --> multiple features
###### Vector is acceptable and treated as one feature

aw.fisher.pvalue <- function(p.values, method, log=FALSE, weight.matrix=FALSE) {
    if(NCOL(p.values) == 1) p.values= t(p.values)

    n = NCOL(p.values)

    pCumLog= -log(p.values)

    if(weight.matrix) {
        orderMatrix = t(apply(p.values, 1, order))
        for(i in 1:NROW(p.values)) {
            pCumLog[i,] = cumsum(pCumLog[i,orderMatrix[i,]])
        }
    } else {
        for(i in 1:NROW(p.values)) {
            pCumLog[i,] = cumsum(pCumLog[i,order(-pCumLog[i,])])
        }
    }

    sumWeight = rep(1, NROW(p.values))
    if(method == "original") {
        bestStat = -pchisq(2*pCumLog[,1], 2, lower.tail=F, log=T)
        for(i in 2:n) {
            statNew = -pchisq(2*pCumLog[,i], 2*i, lower.tail=F, log=T)
            sumWeight[statNew > bestStat] = i
            bestStat = pmax(bestStat, statNew)
        }
    } else if(method == "uncond") {
        bestStat = -pbeta(exp(-pCumLog[,1]), 1, n, log=T)
        statNew =  -pchisq(2*pCumLog[,n], 2*n, lower.tail=F, log=T)
        sumWeight[statNew > bestStat] = n
        bestStat = pmax(bestStat, statNew)
        
        if(n>2) for(i in 2:(n-1)) { 
            statNew = uncondTopDist(2*pCumLog[,i],n,i)
            sumWeight[statNew > bestStat] = i
            bestStat = pmax(bestStat, statNew)
        }

    } else if (method == "cond") {
        bestStat = -pchisq(2*pCumLog[,n], 2*n, lower.tail=F, log=T)
        for(i in 1:(n-1)) {
            statNew = -pchisq(2*(pCumLog[,i]-i*(pCumLog[,i+1]-pCumLog[,i])), 2*i, lower.tail=F, log=T)
            sumWeight[statNew > bestStat] = i
            bestStat = pmax(bestStat, statNew)
        }
    } else {
          stop("Incorrect method!")
    }

    if(weight.matrix) {
        for(i in 1:NROW(orderMatrix)) {
            orderMatrix[i,orderMatrix[i,1:sumWeight[i]]] = 0
        }
        orderMatrix[orderMatrix != 0] = 1
        list(pvalues = aw.fisher.stat(bestStat, n, method, log), weights = 1-orderMatrix, sum.weight = sumWeight)
    } else {
        list(pvalues = aw.fisher.stat(bestStat, n, method, log), sum.weight = sumWeight)
    }
}

### Other Internal Functions #############
aw.fisher.stat <- function(pstat, n, method, log=FALSE) {
    index = match(n, awFisherData[["nList"]])
#    if(!is.na(index)) {
#        quant = awFisherData[[method]][index,]
#    } else {
        ### smooth estimation over n
        numN = NCOL(awFisherData[[method]])
        quant = rep(0, numN)
        for(i in 1:numN) {
            f = splinefun(c(1, awFisherData[["nList"]]), c(awFisherData[["logPTarget"]][i], awFisherData[[method]][,i]))
            quant[i] = f(n)
        }
#    }

    ##### Estimating ###########
    f = splinefun(c(0,quant), c(0,awFisherData[["logPTarget"]]), method="monoH.FC")
    if(log)  -f(pstat)  else exp(-f(pstat))
}


preDefList = exp(log(501)*(1:50)/50)-1
uncondTopDist <- function(stat, n, which) {
    estimates = pvalueTop(A=preDefList, n=n, m=which)
    f = splinefun(c(0,preDefList), c(0,estimates), method="monoH.FC")
    f(stat)
}

intFunc <- function(p,A,n,m) {
    pchisq(A+2*m*log(p), 2*m, lower.tail=F)*dbeta(p, m+1, n-m)
}
pvalueTop <- Vectorize(function(A,n,m) {
                pLim = exp(-A/(2*m))
                -log(integrate(intFunc, lower=pLim,upper=1,A = A, n=n,m=m, stop.on.error=F,rel.tol=1e-3)$value 
                        + pbeta(pLim, m+1, n-m))}
             )

######################
####  Data 
######################

awFisherData = list(logPTarget=-log(c(0.99, 0.95,0.9,0.8,0.7,0.5,0.3,1e-1,1e-2,1e-3,1e-5,1e-7,1e-10, 1e-15, 1e-25, 1e-80,1e-200)),
          original=matrix(c(0.0991299,0.2307506,0.3616958,0.4908891,0.6012938,0.7007693,0.8010211,0.8852742,0.9674844,
                              1.113856,1.304995,1.612286,2.274926,3.816328,6.379592,10.03372,15.8871,28.39492,50.09994,107.3815,
                            0.2532749,0.4589056,0.6426983,0.7975563,0.9350853,1.056027,1.165338,1.267054,1.36635,1.545324,
                              1.810291,2.243907,3.162265,5.079755,8.110184,12.3064,18.90731,32.46575,55.64708,115.5134,
                            0.3785702,0.6240129,0.8244336,0.9922922,1.142877,1.27297,1.397758,1.508944,1.621461,1.83459,
                              2.140719,2.646127,3.666674,5.782382,9.040087,13.49952,20.3788,34.47629,58.35517,119.4172,
                            0.5911318,0.8747565,1.10027,1.286399,1.455918,1.606457,1.752032,1.885532,2.017578,2.267131,
                              2.638745,3.232473,4.391949,6.732669,10.29312,15.05824,22.34432,37.068,61.78113,124.354,
                            0.7887108,1.100477,1.346518,1.550477,1.736986,1.905264,2.067322,2.21877,2.36439,2.644688,
                              3.053511,3.714803,4.978453,7.494532,11.26784,16.27226,23.84434,38.99281,64.31994,128.0403,
                            1.21693,1.582346,1.865415,2.105853,2.32561,2.519655,2.706602,2.884535,3.060112,3.404064,3.876032,
                              4.631569,6.076903,8.897185,13.01535,18.40227,26.42127,42.35352,68.65255,134.2008,
                            1.812172,2.224871,2.54456,2.815906,3.066485,3.299503,3.521331,3.727825,3.931819,4.326293,
                              4.869171,5.727013,7.353057,10.47594,14.96737,20.72703,29.19919,45.86845,73.2044,140.5048,
                            3.006446,3.498343,3.882972,4.202215,4.504013,4.763724,5.026574,5.270521,5.514426,5.981952,6.619728,
                              7.632539,9.526118,13.06491,18.04471,24.38708,33.53347,51.32331,80.00293,149.9905,
                            5.412572,5.978331,6.421937,6.818137,7.178373,7.489069,7.806217,8.101734,8.40125,8.966433,
                              9.724746,10.9367,13.20626,17.32554,22.98528,30.12292,40.1709,59.49315,90.07466,163.6995,
                            7.771121,8.360937,8.862098,9.283587,9.705136,10.04135,10.41136,10.72468,11.07541,11.6783,
                              12.54049,13.87575,16.40208,20.96055,27.11554,34.78738,45.51271,65.95685,97.89573,174.1557,
                            12.44521,13.09599,13.6526,14.10011,14.54778,14.986,15.30801,15.77355,16.13772,16.79953,17.82213,
                              19.38547,22.27196,27.42194,34.3746,42.84719,54.5786,76.72491,110.749,191.039,
                            17.09525,17.75282,18.3716,18.88985,19.37994,19.78108,20.14133,20.75287,21.01495,21.71198,
                              22.83626,24.6096,27.71989,33.2676,40.9529,49.94472,62.55643,85.9504,121.8311,205.3097,
                            24.03921,24.68517,25.40527,25.95085,26.45911,26.87557,27.39865,27.67894,28.19278,29.0289,
                              30.23832,32.32301,35.47315,41.61328,49.97596,59.59598,73.60026,98.42317,136.9592,224.2329,
                            35.56544,36.24341,36.8793,37.71307,38.04985,39.3861,38.91839,39.42456,39.97605,40.79525,42.19192,
                              44.34679,47.92825,54.43455,64.27763,74.33595,90.52811,116.6218,161.774,248.6012,
                            58.59329,59.29901,59.84986,60.87042,61.81504,61.85231,61.88179,62.51328,63.15607,64.2352,65.50203,
                              68.01901,72.74197,79.66,89.91935,101.7426,118.9691,150.7868,192.5185,292.6912,
                            185.1815,186.014,186.6892,187.0414,187.5587,187.9156,189.0134,189.5565,190.0352,192.4073,
                              192.6559,194.5913,201.7345,209.2415,223.7879,239.5384,258.8478,297.156,352.8421,474.5017,
                            461.2872,462.3472,462.4968,462.8381,463.8254,463.8606,464.7668,464.4796,465.2994,465.9452,
                              467.9054,471.0704,474.8733,484.9073,498.1799,512.5592,541.8607,585.8783,651.6489,798.9819), ncol=17),
            cond=matrix(c(0.03810959,0.0621685,0.08131234,0.09920804,0.111354,0.122629,0.1349937,0.1432366,0.1520866,
                              0.1679167,0.184224,0.2087301,0.2434008,0.2871737,0.331107,0.3677778,0.4102335,0.4527058,0.5024291,0.5630171,
                           0.127912,0.1818843,0.2207695,0.2516294,0.2775123,0.2983333,0.3172555,0.337089,0.3547296,0.3786171,
                              0.4092174,0.4455226,0.4985854,0.5779685,0.6441868,0.69751,0.7509269,0.827722,0.8875773,0.9727183,
                           0.2186853,0.2913149,0.3406531,0.3804417,0.4133604,0.4433051,0.4685546,0.4876554,0.5078706,0.541914,0.579527,
                              0.6256311,0.6896221,0.7820427,0.8602653,0.9251285,0.9848993,1.076998,1.144793,1.237292,
                           0.3893207,0.4826316,0.5499885,0.6051856,0.6505681,0.6844553,0.7160721,0.7424228,0.7674958,0.8070018,
                              0.8569188,0.9147883,0.9964052,1.10293,1.196161,1.270607,1.347474,1.446223,1.527901,1.636886,
                           0.5585807,0.6742425,0.7560492,0.8207702,0.8732926,0.9129082,0.9493195,0.9797121,1.008026,1.05531,
                              1.111715,1.180343,1.269962,1.391778,1.492729,1.57213,1.662787,1.768504,1.850884,1.979766,
                           0.95566,1.110985,1.216068,1.29076,1.355766,1.405023,1.445523,1.484516,1.518811,1.580653,1.645448,
                              1.726037,1.839056,1.978376,2.099252,2.193948,2.290192,2.411119,2.51721,2.651896,
                           1.52536,1.713372,1.832813,1.918937,1.993254,2.056027,2.107726,2.152054,2.201163,2.271805,2.348409,
                              2.440764,2.568491,2.726969,2.870715,2.971539,3.078847,3.214829,3.323123,3.47035,
                           2.70562,2.929678,3.080541,3.193354,3.277625,3.338019,3.4043,3.450513,3.511548,3.599451,3.674784,
                              3.804746,3.947966,4.123648,4.289046,4.406467,4.52903,4.700287,4.825234,4.960964,
                           5.077873,5.34332,5.526001,5.657974,5.772755,5.848972,5.931133,6.003743,6.089314,6.176835,
                              6.254821,6.350939,6.530929,6.826697,6.97927,7.100861,7.447727,7.758332,7.515395,7.543793,
                           7.415679,7.711341,7.889183,8.047575,8.198342,8.252384,8.337884,8.467926,8.543158,8.638408,
                              8.671341,8.793482,9.047245,9.584596,9.452446,9.45232,9.892516,10.00523,9.830353,10.04764,
                           12.0838,12.39359,12.6006,12.75797,12.90931,13.02841,13.0343,13.156,13.22006,13.26682,13.39121,
                              13.35111,13.73286,14.03898,14.02277,14.26102,14.27311,14.44403,14.33591,14.14482,
                           16.70767,17.00769,17.31219,17.4414,17.57655,17.72012,17.67593,17.99483,17.70582,17.76762,
                              18.03914,18.0793,18.0617,18.34354,19.13029,20.76382,18.30004,18.7526,18.83723,18.87949,
                           23.64351,23.91838,24.24045,24.47294,24.67672,24.56047,24.69437,24.52458,24.53108,24.68174,
                              24.84996,25.85694,25.28258,25.2019,25.52906,25.05026,25.52343,26.48087,25.91837,25.62153,
                           35.17247,35.44539,35.75883,36.02557,35.93328,36.91222,35.91838,35.94448,36.09843,35.89931,
                              37.04063,36.88587,36.29341,36.53636,42.25099,36.58584,36.77084,36.50841,41.56533,37.13722,
                           58.20533,58.44961,58.59536,59.20801,59.48873,59.02293,58.63219,58.86448,58.81296,58.91426,
                              59.27569,59.34844,60.65453,59.62144,59.50907,59.61764,59.74752,59.77,59.80256,66.04942,
                           184.7801,185.2772,185.0311,185.1535,185.1581,185.0703,185.2399,185.5662,185.5465,186.0258,
                              185.8041,185.8345,185.7155,186.1403,186.0538,186.0937,186.0333,186.2792,186.202,186.6681,
                           460.8942,461.4731,461.2796,461.2338,461.6603,461.3963,461.3641,461.362,461.3927,462.4369,462.3063,
                              461.5865,461.6614,461.646,461.8454,462.9842,462.1843,462.2884,462.3873,463.0256), ncol=17),
           uncond=matrix(c(0.01167915,0.01380246,0.01517486,0.01701253,0.01780689,0.01869392,0.02015334,0.02096671,0.02257714,0.02432555,
                             0.02654229,0.02890041,0.03275772,0.04085687,0.04920061,0.0568219,0.06522525,0.09740171,0.1110982,0.1383659,
                           0.06097223,0.06811177,0.07404343,0.07932669,0.08442349,0.0883944,0.09298591,0.0964361,0.09975088,0.1070759,
                             0.1121707,0.1199021,0.1349449,0.1580033,0.1787628,0.1968977,0.2150478,0.2727937,0.2969215,0.342652,
                           0.1227619,0.1360761,0.1473208,0.1565375,0.1659376,0.1737538,0.1820478,0.1874003,0.1911773,0.1991744,0.2086036,
                             0.2226894,0.245687,0.2794811,0.3077283,0.3285149,0.3587663,0.4260081,0.4548423,0.5025304,
                           0.2569734,0.2807185,0.3004544,0.3166578,0.3310639,0.3423149,0.3539376,0.3638423,0.3707911,0.3849804,
                             0.4027341,0.4267229,0.4548959,0.4978543,0.5393592,0.5693082,0.6173452,0.6846767,0.7185779,0.7757726,
                           0.4063705,0.4418796,0.4701133,0.4912607,0.510468,0.528245,0.5394942,0.5477942,0.5610035,0.5793259,0.6046962,
                             0.6299256,0.6682264,0.7182808,0.7674903,0.8059049,0.8593518,0.9323056,0.9668378,1.033236,
                           0.7807272,0.843123,0.8808574,0.9126672,0.9396796,0.960439,0.978473,0.9922291,1.008917,1.038971,1.071566,
                             1.107042,1.151322,1.220174,1.283471,1.330622,1.381731,1.465399,1.510273,1.570922,
                           1.342772,1.422975,1.47796,1.514661,1.556413,1.586022,1.608228,1.630456,1.652725,1.693232,1.725,
                             1.76851,1.824637,1.904191,1.980362,2.024488,2.089696,2.179498,2.226941,2.267269,
                           2.529015,2.652573,2.730094,2.785313,2.832947,2.862264,2.886662,2.909889,2.945636,2.987453,3.027188,
                             3.086178,3.145357,3.244686,3.322164,3.375402,3.441574,3.537447,3.594048,3.606768,
                           4.940525,5.100422,5.211391,5.288603,5.345469,5.379895,5.410683,5.44296,5.506372,5.538015,5.614381,
                             5.612352,5.736223,5.811657,5.860137,5.914662,6.229093,6.139462,6.136908,6.083464,
                           7.308872,7.500962,7.616301,7.699705,7.776036,7.817443,7.831599,7.91947,7.970223,8.03098,8.046011,
                             8.037557,8.229807,8.207456,8.27931,8.353512,8.695223,8.395307,8.444783,8.332947,
                           11.99451,12.21962,12.37417,12.44409,12.53803,12.5802,12.60703,12.69658,12.79327,12.80341,12.72086,
                             12.68056,12.90105,12.95976,12.90918,13.10627,12.74418,12.65242,12.97392,12.54703,
                           16.64509,16.87218,17.06891,17.19426,17.28926,17.30334,17.25235,17.70002,17.41296,17.31092,17.42144,
                             17.36365,17.58575,17.48862,18.22938,17.31067,17.21504,17.31552,16.94165,17.2159,
                           23.59017,23.81475,24.08688,24.22897,24.25407,24.15445,24.6323,24.20535,24.21415,24.25998,24.17379,
                             24.42562,24.20963,24.07561,24.17858,24.3309,24.21982,23.96657,23.18492,23.19917,
                           35.12209,35.38765,35.54676,35.93301,35.72465,37.15706,35.62476,35.65893,35.79171,35.74076,35.81368,
                             35.84007,35.68559,35.55025,35.65699,35.59003,36.08693,35.01463,34.55093,34.54058,
                           58.16783,58.45219,58.56312,59.10713,59.24837,58.97294,58.55741,58.68723,58.78797,58.88322,58.60795,
                             59.02749,58.99561,58.86774,58.73229,58.48511,57.83821,57.56131,57.55227,57.56808,
                           184.7882,185.2093,185.3756,185.2619,185.2936,185.2547,185.5489,185.53,185.5679,186.3204,185.5461,185.3334,
                             186.3847,184.6643,184.2206,184.1881,184.2057,184.22,184.2071,184.2049,
                           460.981,461.4302,461.0973,460.9452,461.0463,460.7887,460.7704,460.6247,460.6181,460.5617,460.5465,460.5333,
                             460.538,460.5324,460.5328,460.5398,460.51,460.52,460.531,460.5183),ncol=17),
          nList = c(2:10,12,15,20,30,50,80,120,180,300,500,1000),
          totalN = 100000)