#Summary: # --- R FUNCTIONS: ------- # (1) for simulation: # (i) gen.f # generate a two-locus penetrance matrix # (ii) gen.data # generate two-locus case-control data, data is in a list # the two components are for cases and controsl, respectively, # each is a 3x3 matrix of counts on the two loci. # (2) Main function: pc() # (i) IT #interaction test, it is used in function pc. # (ii) pc # partitioning chisq. # (iii) there are other two supporting functions: mychisq(), expected(), used in main function. # ---- DATA FORMAT: ------- # 2x9 table (two rows for case/control, 9 columns for 9 genotype combinations # example: #> data1 # aa bb Aa bb AA bb aa Bb Aa Bb AA Bb aa BB Aa BB AA BB #case 179 240 81 132 207 61 35 45 20 #control 189 238 71 161 185 67 38 40 11 # ----- RUN THE PROGRAM ------ # pc(data) ## run the pc program. #################################################################################### ## Generate a 3x3 penetrance matrix (two-locus penetrance model) using logistic model: gen.f=function(alpha=-2, beta1=0.2, beta2=0.1, gamma=0.3, x1=0, x2=0){ # penetrance model is # logit(f)=alpha+beta1*x1+beta2*x2+gamma*x1*x2. # x1,x2: trend score for Aa,Bb # For locus 1, it is recessive when x1=0, dominant if x1=1, additive if x1=0.5 # similarly for x2 (for locus 2). # beta1, beta2: main effect at locus 1 and 2; gamma: interaction # alpha: logarithm of population risk z1=c(0,x1, 1); z2=c(0,x2,1) z12=z1%*%t(z2) alpha+outer(beta1*z1,beta2*z2, "+")+gamma*z12 ->tmp f=exp(tmp)/(1+exp(tmp)) f } # end of generating penetrance f. #################################################################################### ## Generate two-locus genotype counts for case and controls: gen.data =function(n1=1000, n0=1000, f=matrix(0.1, 3,3), p1=0.3, p2=0.4){ ## n1: numb case; n0:control; ## F1, F0: inbreeding coef for case and control ##f: penetrance matrix 3x3 ## p1 ,p2 population allele freq for locus 1 and locus2 #################################################### rmn=function(n,size,prob){ U=matrix( runif(n*size),n, size) f.local=function(u, prob){ k=length(prob) numb=1:k cp=c(0, cumsum(prob)) for (i in 1:k){ numb[i]=sum(cp[i]tmp p.case =tmp/sum(tmp) as.vector(p.case)->p.case rmn(1, n1, p.case) -> geno.case matrix(geno.case,nrow=3) ->geno.case R*(1-f)->tmp p.case =tmp/sum(tmp) as.vector(p.case)->p.con rmn(1, n1, p.case) -> geno.con matrix(geno.con,nrow=3) ->geno.con list(case=geno.case, control=geno.con) } #end of gen.data function #################################################################### ## supporting functions.. ## chi-square test for any 2-way table: mychisq=function(x){ n=sum(x) apply(x,1,sum)->row.sum apply(x,2,sum)->col.sum e=row.sum%*%t(col.sum)/n e[e==0]=0.5 sum((x-e)^2/e) } #### ## expected counts for a contingency table: expected= function(x){ n=sum(x) apply(x,1, sum)/n ->p.r apply(x,2, sum)/n ->p.c p.r%*%t(p.c)*n->e delta=x-e e/n } # end of expected. ######################################################################## ################## Interaction test #################################### ######################################################################## ## Interaction test (IT), but also provide total chisq, Marginal tests at Locus 1 ## (L1) and Locus 2 (L2), also provide main effect test ## main effect test = total chisq - interaction test. IT=function(data , alpha1=0.1) { n1=sum(data[[1]]) n0=sum(data[[2]]) case=mychisq(data[[1]]); control=mychisq(data[[2]]); combined=mychisq(data[[1]]+data[[2]]) ## deviation from independence: d=delta = p(ij)-pi*pj: exp.case=expected(data[[1]]); exp.con =expected(data[[2]]) data[[1]]/n1-exp.case -> r.case; data[[2]]/n0-exp.con -> r.con d=r.case-r.con v=(exp.case/n1+exp.con/n0) d/sqrt(v)->residuals interact = d^2/v interaction = sum(interact) # fulldata=rbind(as.vector(data[[1]]), as.vector(data[[2]])) X2=mychisq(fulldata) # totoal chisq main=X2-interaction ## marginal data: mdata1=rbind(apply(data[[1]],1, sum), apply(data[[2]],1,sum)) mdata2=rbind(apply(data[[1]],2, sum), apply(data[[2]],2,sum)) x2.L1=mychisq(mdata1) x2.L2=mychisq(mdata2) ## ## alltest = c(total=X2,main=main, L1=x2.L1,L2= x2.L2, interaction=interaction) return(alltest) }# end IT ######################################################################## ################ MAIN FUNCTION ######################################### ######################################################################## ## MAIN FUNCTION: for partitioning chi-squares.. pc =function(data){ # data can be either a list of two components for cases and controls respectively # or a 2x9 table, first row for cases and second row for controls for the 9 genotype combinations. if (!is.list(data)) data=list(case=matrix(data[1,],3,3), control=matrix(data[2,],3,3)) case=data[[1]];cont=data[[2]] case1=case[1:2,1:2];cont1 =cont[1:2,1:2] par1=IT( list(case1, cont1)) i=1:2 case2=cbind(case[i,1]+case[i,2], case[i,3]) cont2=cbind(cont[i,1]+cont[i,2], cont[i,3]) par2=IT(list(case2, cont2)) i=1:2 case3=cbind(case[1,i]+case[2,i], case[3,i]) cont3=cbind(cont[1,i]+cont[2,i], cont[3,i]) par3=IT(list(case3, cont3)) i=1:2 case4 =matrix( c(sum(case[i,i]), case[3,1]+case[3,2], case[1,3]+case[2,3],case[3,3]),2,2) cont4 =matrix( c(sum(cont[i,i]), cont[3,1]+cont[3,2], cont[1,3]+cont[2,3],cont[3,3]),2,2) par4=IT(list(case4, cont4)) tmp = rbind(par1,par2,par3,par4) IT(data)->total tmp[,5]->pars c( total, pars )->all.chisq all.chisq=all.chisq[c(1,2,5, 6:9,3:4)] return(all.chisq) } #end of the main function pc ############# END #########################################################