|
0
|
1 ##############################################################################################################################################
|
|
|
2 ###MLMM - Multi-Locus Mixed Model
|
|
|
3 ###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH
|
|
|
4 #######
|
|
|
5 #
|
|
|
6 ##note: require EMMA
|
|
|
7 #library(emma)
|
|
|
8 #source('emma.r')
|
|
|
9 #
|
|
|
10 ##REQUIRED DATA & FORMAT
|
|
|
11 #
|
|
|
12 #PHENOTYPE - Y: a vector of length m, with names(Y)=individual names
|
|
|
13 #GENOTYPE - X: a n by m matrix, where n=number of individuals, m=n umber of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names
|
|
|
14 #KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names
|
|
|
15 #each of these data being sorted in the same way, according to the individual name
|
|
|
16 #
|
|
|
17 ##FOR PLOTING THE GWAS RESULTS
|
|
|
18 #SNP INFORMATION - snp_info: a data frame having at least 3 columns:
|
|
|
19 # - 1 named 'SNP', with SNP names (same as colnames(X)),
|
|
|
20 # - 1 named 'Chr', with the chromosome number to which belong each SNP
|
|
|
21 # - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.
|
|
|
22 #######
|
|
|
23 #
|
|
|
24 ##FUNCTIONS USE
|
|
|
25 #save this file somewhere on your computer and source it!
|
|
|
26 #source('path/mlmm.r')
|
|
|
27 #
|
|
|
28 ###FORWARD + BACKWARD ANALYSES
|
|
|
29 #mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)
|
|
|
30 #X,Y,K as described above
|
|
|
31 #nbchunks: an integer defining the number of chunks of to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory
|
|
|
32 #maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,
|
|
|
33 # however to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.
|
|
|
34 # It's value must be specified as an integer >= 3
|
|
|
35 #
|
|
|
36 ###RESULTS
|
|
|
37 #
|
|
|
38 ##STEPWISE TABLE
|
|
|
39 #mygwas$step_table
|
|
|
40 #
|
|
|
41 ##PLOTS
|
|
|
42 #
|
|
|
43 ##PLOTS FORM THE FORWARD TABLE
|
|
|
44 #plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))
|
|
|
45 #
|
|
|
46 ##RSS PLOT
|
|
|
47 #plot_step_RSS(mygwas)
|
|
|
48 #
|
|
|
49 ##GWAS MANHATTAN PLOTS
|
|
|
50 #
|
|
|
51 #FORWARD STEPS
|
|
|
52 #plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)
|
|
|
53 #step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)
|
|
|
54 #snp_info as described above
|
|
|
55 #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
|
|
|
56 #
|
|
|
57 #OPTIMAL MODELS
|
|
|
58 #Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria
|
|
|
59 #
|
|
|
60 #plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)
|
|
|
61 #snp_info as described above
|
|
|
62 #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
|
|
|
63 #
|
|
|
64 ##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST
|
|
|
65 #plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)
|
|
|
66 #step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)
|
|
|
67 #snp_info as described above
|
|
|
68 #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
|
|
|
69 #chrom is an integer specifying the chromosome on which the region of interest is
|
|
|
70 #pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info
|
|
|
71 #
|
|
|
72 #plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)
|
|
|
73 #snp_info as described above
|
|
|
74 #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
|
|
|
75 #chrom is an integer specifying the chromosome on which the region of interest is
|
|
|
76 #pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info
|
|
|
77 #
|
|
|
78 ##QQPLOTS of pvalues
|
|
|
79 #qqplot_fwd_GWAS(mygwas,nsteps)
|
|
|
80 #nsteps=maximum number of forward steps to be displayed
|
|
|
81 #
|
|
|
82 #qqplot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'))
|
|
|
83 #
|
|
|
84 ##############################################################################################################################################
|
|
|
85
|
|
|
86 mlmm<-function(Y,X,K,nbchunks,maxsteps) {
|
|
|
87
|
|
|
88 n<-length(Y)
|
|
|
89 m<-ncol(X)
|
|
|
90
|
|
|
91 stopifnot(ncol(K) == n)
|
|
|
92 stopifnot(nrow(K) == n)
|
|
|
93 stopifnot(nrow(X) == n)
|
|
|
94 stopifnot(nbchunks >= 2)
|
|
|
95 stopifnot(maxsteps >= 3)
|
|
|
96
|
|
|
97 #INTERCEPT
|
|
|
98
|
|
|
99 Xo<-rep(1,n)
|
|
|
100
|
|
|
101 #K MATRIX NORMALISATION
|
|
|
102
|
|
|
103 K_norm<-(n-1)/sum((diag(n)-matrix(1,n,n)/n)*K)*K
|
|
|
104 rm(K)
|
|
|
105
|
|
|
106 #step 0 : NULL MODEL
|
|
|
107 cof_fwd<-list()
|
|
|
108 cof_fwd[[1]]<-as.matrix(Xo)
|
|
|
109 colnames(cof_fwd[[1]])<-'Xo'
|
|
|
110
|
|
|
111 mod_fwd<-list()
|
|
|
112 mod_fwd[[1]]<-emma.REMLE(Y,cof_fwd[[1]],K_norm)
|
|
|
113
|
|
|
114 herit_fwd<-list()
|
|
|
115 herit_fwd[[1]]<-mod_fwd[[1]]$vg/(mod_fwd[[1]]$vg+mod_fwd[[1]]$ve)
|
|
|
116
|
|
|
117 RSSf<-list()
|
|
|
118 RSSf[[1]]<-'NA'
|
|
|
119
|
|
|
120 RSS_H0<-list()
|
|
|
121 RSS_H0[[1]]<-'NA'
|
|
|
122
|
|
|
123 df1<-1
|
|
|
124 df2<-list()
|
|
|
125 df2[[1]]<-'NA'
|
|
|
126
|
|
|
127 Ftest<-list()
|
|
|
128 Ftest[[1]]<-'NA'
|
|
|
129
|
|
|
130 pval<-list()
|
|
|
131 pval[[1]]<-'NA'
|
|
|
132
|
|
|
133 fwd_lm<-list()
|
|
|
134
|
|
|
135 cat('null model done! pseudo-h=',round(herit_fwd[[1]],3),'\n')
|
|
|
136
|
|
|
137 #step 1 : EMMAX
|
|
|
138
|
|
|
139 M<-solve(chol(mod_fwd[[1]]$vg*K_norm+mod_fwd[[1]]$ve*diag(n)))
|
|
|
140 Y_t<-crossprod(M,Y)
|
|
|
141 cof_fwd_t<-crossprod(M,cof_fwd[[1]])
|
|
|
142 fwd_lm[[1]]<-summary(lm(Y_t~0+cof_fwd_t))
|
|
|
143 Res_H0<-fwd_lm[[1]]$residuals
|
|
|
144 Q_<-qr.Q(qr(cof_fwd_t))
|
|
|
145
|
|
|
146 RSS<-list()
|
|
|
147 for (j in 1:(nbchunks-1)) {
|
|
|
148 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[1]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
|
|
|
149 RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
150 rm(X_t)}
|
|
|
151 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[1]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[1]])-1))])
|
|
|
152 RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
153 rm(X_t,j)
|
|
|
154
|
|
|
155 RSSf[[2]]<-unlist(RSS)
|
|
|
156 RSS_H0[[2]]<-sum(Res_H0^2)
|
|
|
157 df2[[2]]<-n-df1-ncol(cof_fwd[[1]])
|
|
|
158 Ftest[[2]]<-(rep(RSS_H0[[2]],length(RSSf[[2]]))/RSSf[[2]]-1)*df2[[2]]/df1
|
|
|
159 pval[[2]]<-pf(Ftest[[2]],df1,df2[[2]],lower.tail=FALSE)
|
|
|
160
|
|
|
161 cof_fwd[[2]]<-cbind(cof_fwd[[1]],X[,colnames(X) %in% names(which(RSSf[[2]]==min(RSSf[[2]]))[1])])
|
|
|
162 colnames(cof_fwd[[2]])<-c(colnames(cof_fwd[[1]]),names(which(RSSf[[2]]==min(RSSf[[2]]))[1]))
|
|
|
163 mod_fwd[[2]]<-emma.REMLE(Y,cof_fwd[[2]],K_norm)
|
|
|
164 herit_fwd[[2]]<-mod_fwd[[2]]$vg/(mod_fwd[[2]]$vg+mod_fwd[[2]]$ve)
|
|
|
165 rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)
|
|
|
166
|
|
|
167 cat('step 1 done! pseudo-h=',round(herit_fwd[[2]],3),'\n')
|
|
|
168
|
|
|
169 #FORWARD
|
|
|
170
|
|
|
171 for (i in 3:(maxsteps)) {
|
|
|
172 if (herit_fwd[[i-2]] < 0.01) break else {
|
|
|
173
|
|
|
174 M<-solve(chol(mod_fwd[[i-1]]$vg*K_norm+mod_fwd[[i-1]]$ve*diag(n)))
|
|
|
175 Y_t<-crossprod(M,Y)
|
|
|
176 cof_fwd_t<-crossprod(M,cof_fwd[[i-1]])
|
|
|
177 fwd_lm[[i-1]]<-summary(lm(Y_t~0+cof_fwd_t))
|
|
|
178 Res_H0<-fwd_lm[[i-1]]$residuals
|
|
|
179 Q_ <- qr.Q(qr(cof_fwd_t))
|
|
|
180
|
|
|
181 RSS<-list()
|
|
|
182 for (j in 1:(nbchunks-1)) {
|
|
|
183 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[i-1]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
|
|
|
184 RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
185 rm(X_t)}
|
|
|
186 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[i-1]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[i-1]])-1))])
|
|
|
187 RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
188 rm(X_t,j)
|
|
|
189
|
|
|
190 RSSf[[i]]<-unlist(RSS)
|
|
|
191 RSS_H0[[i]]<-sum(Res_H0^2)
|
|
|
192 df2[[i]]<-n-df1-ncol(cof_fwd[[i-1]])
|
|
|
193 Ftest[[i]]<-(rep(RSS_H0[[i]],length(RSSf[[i]]))/RSSf[[i]]-1)*df2[[i]]/df1
|
|
|
194 pval[[i]]<-pf(Ftest[[i]],df1,df2[[i]],lower.tail=FALSE)
|
|
|
195
|
|
|
196 cof_fwd[[i]]<-cbind(cof_fwd[[i-1]],X[,colnames(X) %in% names(which(RSSf[[i]]==min(RSSf[[i]]))[1])])
|
|
|
197 colnames(cof_fwd[[i]])<-c(colnames(cof_fwd[[i-1]]),names(which(RSSf[[i]]==min(RSSf[[i]]))[1]))
|
|
|
198 mod_fwd[[i]]<-emma.REMLE(Y,cof_fwd[[i]],K_norm)
|
|
|
199 herit_fwd[[i]]<-mod_fwd[[i]]$vg/(mod_fwd[[i]]$vg+mod_fwd[[i]]$ve)
|
|
|
200 rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)}
|
|
|
201 cat('step ',i-1,' done! pseudo-h=',round(herit_fwd[[i]],3),'\n')}
|
|
|
202 rm(i)
|
|
|
203
|
|
|
204 ##gls at last forward step
|
|
|
205 M<-solve(chol(mod_fwd[[length(mod_fwd)]]$vg*K_norm+mod_fwd[[length(mod_fwd)]]$ve*diag(n)))
|
|
|
206 Y_t<-crossprod(M,Y)
|
|
|
207 cof_fwd_t<-crossprod(M,cof_fwd[[length(mod_fwd)]])
|
|
|
208 fwd_lm[[length(mod_fwd)]]<-summary(lm(Y_t~0+cof_fwd_t))
|
|
|
209
|
|
|
210 Res_H0<-fwd_lm[[length(mod_fwd)]]$residuals
|
|
|
211 Q_ <- qr.Q(qr(cof_fwd_t))
|
|
|
212
|
|
|
213 RSS<-list()
|
|
|
214 for (j in 1:(nbchunks-1)) {
|
|
|
215 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[length(mod_fwd)]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
|
|
|
216 RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
217 rm(X_t)}
|
|
|
218 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[length(mod_fwd)]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[length(mod_fwd)]])-1))])
|
|
|
219 RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
220 rm(X_t,j)
|
|
|
221
|
|
|
222 RSSf[[length(mod_fwd)+1]]<-unlist(RSS)
|
|
|
223 RSS_H0[[length(mod_fwd)+1]]<-sum(Res_H0^2)
|
|
|
224 df2[[length(mod_fwd)+1]]<-n-df1-ncol(cof_fwd[[length(mod_fwd)]])
|
|
|
225 Ftest[[length(mod_fwd)+1]]<-(rep(RSS_H0[[length(mod_fwd)+1]],length(RSSf[[length(mod_fwd)+1]]))/RSSf[[length(mod_fwd)+1]]-1)*df2[[length(mod_fwd)+1]]/df1
|
|
|
226 pval[[length(mod_fwd)+1]]<-pf(Ftest[[length(mod_fwd)+1]],df1,df2[[length(mod_fwd)+1]],lower.tail=FALSE)
|
|
|
227 rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)
|
|
|
228
|
|
|
229 ##get max pval at each forward step
|
|
|
230 max_pval_fwd<-vector(mode="numeric",length=length(fwd_lm))
|
|
|
231 max_pval_fwd[1]<-0
|
|
|
232 for (i in 2:length(fwd_lm)) {max_pval_fwd[i]<-max(fwd_lm[[i]]$coef[2:i,4])}
|
|
|
233 rm(i)
|
|
|
234
|
|
|
235 ##get the number of parameters & Loglikelihood from ML at each step
|
|
|
236 mod_fwd_LL<-list()
|
|
|
237 mod_fwd_LL[[1]]<-list(nfixed=ncol(cof_fwd[[1]]),LL=emma.MLE(Y,cof_fwd[[1]],K_norm)$ML)
|
|
|
238 for (i in 2:length(cof_fwd)) {mod_fwd_LL[[i]]<-list(nfixed=ncol(cof_fwd[[i]]),LL=emma.MLE(Y,cof_fwd[[i]],K_norm)$ML)}
|
|
|
239 rm(i)
|
|
|
240
|
|
|
241 cat('backward analysis','\n')
|
|
|
242
|
|
|
243 ##BACKWARD (1st step == last fwd step)
|
|
|
244
|
|
|
245 dropcof_bwd<-list()
|
|
|
246 cof_bwd<-list()
|
|
|
247 mod_bwd <- list()
|
|
|
248 bwd_lm<-list()
|
|
|
249 herit_bwd<-list()
|
|
|
250
|
|
|
251 dropcof_bwd[[1]]<-'NA'
|
|
|
252 cof_bwd[[1]]<-as.matrix(cof_fwd[[length(mod_fwd)]][,!colnames(cof_fwd[[length(mod_fwd)]]) %in% dropcof_bwd[[1]]])
|
|
|
253 colnames(cof_bwd[[1]])<-colnames(cof_fwd[[length(mod_fwd)]])[!colnames(cof_fwd[[length(mod_fwd)]]) %in% dropcof_bwd[[1]]]
|
|
|
254 mod_bwd[[1]]<-emma.REMLE(Y,cof_bwd[[1]],K_norm)
|
|
|
255 herit_bwd[[1]]<-mod_bwd[[1]]$vg/(mod_bwd[[1]]$vg+mod_bwd[[1]]$ve)
|
|
|
256 M<-solve(chol(mod_bwd[[1]]$vg*K_norm+mod_bwd[[1]]$ve*diag(n)))
|
|
|
257 Y_t<-crossprod(M,Y)
|
|
|
258 cof_bwd_t<-crossprod(M,cof_bwd[[1]])
|
|
|
259 bwd_lm[[1]]<-summary(lm(Y_t~0+cof_bwd_t))
|
|
|
260
|
|
|
261 rm(M,Y_t,cof_bwd_t)
|
|
|
262
|
|
|
263 for (i in 2:length(mod_fwd)) {
|
|
|
264 dropcof_bwd[[i]]<-(colnames(cof_bwd[[i-1]])[2:ncol(cof_bwd[[i-1]])])[which(abs(bwd_lm[[i-1]]$coef[2:nrow(bwd_lm[[i-1]]$coef),3])==min(abs(bwd_lm[[i-1]]$coef[2:nrow(bwd_lm[[i-1]]$coef),3])))]
|
|
|
265 cof_bwd[[i]]<-as.matrix(cof_bwd[[i-1]][,!colnames(cof_bwd[[i-1]]) %in% dropcof_bwd[[i]]])
|
|
|
266 colnames(cof_bwd[[i]])<-colnames(cof_bwd[[i-1]])[!colnames(cof_bwd[[i-1]]) %in% dropcof_bwd[[i]]]
|
|
|
267 mod_bwd[[i]]<-emma.REMLE(Y,cof_bwd[[i]],K_norm)
|
|
|
268 herit_bwd[[i]]<-mod_bwd[[i]]$vg/(mod_bwd[[i]]$vg+mod_bwd[[i]]$ve)
|
|
|
269 M<-solve(chol(mod_bwd[[i]]$vg*K_norm+mod_bwd[[i]]$ve*diag(n)))
|
|
|
270 Y_t<-crossprod(M,Y)
|
|
|
271 cof_bwd_t<-crossprod(M,cof_bwd[[i]])
|
|
|
272 bwd_lm[[i]]<-summary(lm(Y_t~0+cof_bwd_t))
|
|
|
273 rm(M,Y_t,cof_bwd_t)}
|
|
|
274
|
|
|
275 rm(i)
|
|
|
276
|
|
|
277 ##get max pval at each backward step
|
|
|
278 max_pval_bwd<-vector(mode="numeric",length=length(bwd_lm))
|
|
|
279 for (i in 1:(length(bwd_lm)-1)) {max_pval_bwd[i]<-max(bwd_lm[[i]]$coef[2:(length(bwd_lm)+1-i),4])}
|
|
|
280 max_pval_bwd[length(bwd_lm)]<-0
|
|
|
281
|
|
|
282 ##get the number of parameters & Loglikelihood from ML at each step
|
|
|
283 mod_bwd_LL<-list()
|
|
|
284 mod_bwd_LL[[1]]<-list(nfixed=ncol(cof_bwd[[1]]),LL=emma.MLE(Y,cof_bwd[[1]],K_norm)$ML)
|
|
|
285 for (i in 2:length(cof_bwd)) {mod_bwd_LL[[i]]<-list(nfixed=ncol(cof_bwd[[i]]),LL=emma.MLE(Y,cof_bwd[[i]],K_norm)$ML)}
|
|
|
286 rm(i)
|
|
|
287
|
|
|
288 cat('creating output','\n')
|
|
|
289
|
|
|
290 ##Forward Table: Fwd + Bwd Tables
|
|
|
291 #Compute parameters for model criteria
|
|
|
292 BIC<-function(x){-2*x$LL+(x$nfixed+1)*log(n)}
|
|
|
293 extBIC<-function(x){BIC(x)+2*lchoose(m,x$nfixed-1)}
|
|
|
294
|
|
|
295 fwd_table<-data.frame(step=ncol(cof_fwd[[1]])-1,step_=paste('fwd',ncol(cof_fwd[[1]])-1,sep=''),cof='NA',ncof=ncol(cof_fwd[[1]])-1,h2=herit_fwd[[1]]
|
|
|
296 ,maxpval=max_pval_fwd[1],BIC=BIC(mod_fwd_LL[[1]]),extBIC=extBIC(mod_fwd_LL[[1]]))
|
|
|
297 for (i in 2:(length(mod_fwd))) {fwd_table<-rbind(fwd_table,
|
|
|
298 data.frame(step=ncol(cof_fwd[[i]])-1,step_=paste('fwd',ncol(cof_fwd[[i]])-1,sep=''),cof=paste('+',colnames(cof_fwd[[i]])[i],sep=''),ncof=ncol(cof_fwd[[i]])-1,h2=herit_fwd[[i]]
|
|
|
299 ,maxpval=max_pval_fwd[i],BIC=BIC(mod_fwd_LL[[i]]),extBIC=extBIC(mod_fwd_LL[[i]])))}
|
|
|
300
|
|
|
301 rm(i)
|
|
|
302
|
|
|
303 bwd_table<-data.frame(step=length(mod_fwd),step_=paste('bwd',0,sep=''),cof=paste('-',dropcof_bwd[[1]],sep=''),ncof=ncol(cof_bwd[[1]])-1,h2=herit_bwd[[1]]
|
|
|
304 ,maxpval=max_pval_bwd[1],BIC=BIC(mod_bwd_LL[[1]]),extBIC=extBIC(mod_bwd_LL[[1]]))
|
|
|
305 for (i in 2:(length(mod_bwd))) {bwd_table<-rbind(bwd_table,
|
|
|
306 data.frame(step=length(mod_fwd)+i-1,step_=paste('bwd',i-1,sep=''),cof=paste('-',dropcof_bwd[[i]],sep=''),ncof=ncol(cof_bwd[[i]])-1,h2=herit_bwd[[i]]
|
|
|
307 ,maxpval=max_pval_bwd[i],BIC=BIC(mod_bwd_LL[[i]]),extBIC=extBIC(mod_bwd_LL[[i]])))}
|
|
|
308
|
|
|
309 rm(i,BIC,extBIC,max_pval_fwd,max_pval_bwd,dropcof_bwd)
|
|
|
310
|
|
|
311 fwdbwd_table<-rbind(fwd_table,bwd_table)
|
|
|
312
|
|
|
313 #RSS for plot
|
|
|
314 mod_fwd_RSS<-vector()
|
|
|
315 mod_fwd_RSS[1]<-sum((Y-cof_fwd[[1]]%*%fwd_lm[[1]]$coef[,1])^2)
|
|
|
316 for (i in 2:length(mod_fwd)) {mod_fwd_RSS[i]<-sum((Y-cof_fwd[[i]]%*%fwd_lm[[i]]$coef[,1])^2)}
|
|
|
317 mod_bwd_RSS<-vector()
|
|
|
318 mod_bwd_RSS[1]<-sum((Y-cof_bwd[[1]]%*%bwd_lm[[1]]$coef[,1])^2)
|
|
|
319 for (i in 2:length(mod_bwd)) {mod_bwd_RSS[i]<-sum((Y-cof_bwd[[i]]%*%bwd_lm[[i]]$coef[,1])^2)}
|
|
|
320 expl_RSS<-c(1-sapply(mod_fwd_RSS,function(x){x/mod_fwd_RSS[1]}),1-sapply(mod_bwd_RSS,function(x){x/mod_bwd_RSS[length(mod_bwd_RSS)]}))
|
|
|
321 h2_RSS<-c(unlist(herit_fwd),unlist(herit_bwd))*(1-expl_RSS)
|
|
|
322 unexpl_RSS<-1-expl_RSS-h2_RSS
|
|
|
323 plot_RSS<-t(apply(cbind(expl_RSS,h2_RSS,unexpl_RSS),1,cumsum))
|
|
|
324
|
|
|
325 #GLS pvals at each step
|
|
|
326 pval_step<-list()
|
|
|
327 pval_step[[1]]<-list(out=data.frame('SNP'=colnames(X),'pval'=pval[[2]]),cof='')
|
|
|
328 for (i in 2:(length(mod_fwd))) {pval_step[[i]]<-list(out=rbind(data.frame(SNP=colnames(cof_fwd[[i]])[-1],'pval'=fwd_lm[[i]]$coef[2:i,4]),
|
|
|
329 data.frame(SNP=colnames(X)[-which(colnames(X) %in% colnames(cof_fwd[[i]]))],'pval'=pval[[i+1]])),cof=colnames(cof_fwd[[i]])[-1])}
|
|
|
330
|
|
|
331 #GLS pvals for best models according to extBIC and mbonf
|
|
|
332
|
|
|
333 opt_extBIC<-fwdbwd_table[which(fwdbwd_table$extBIC==min(fwdbwd_table$extBIC))[1],]
|
|
|
334 opt_mbonf<-(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),])[which(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),]$ncof==max(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),]$ncof))[1],]
|
|
|
335 bestmodel_pvals<-function(model) {if(substr(model$step_,start=0,stop=3)=='fwd') {
|
|
|
336 pval_step[[as.integer(substring(model$step_,first=4))+1]]} else if (substr(model$step_,start=0,stop=3)=='bwd') {
|
|
|
337 cof<-cof_bwd[[as.integer(substring(model$step_,first=4))+1]]
|
|
|
338 mixedmod<-emma.REMLE(Y,cof,K_norm)
|
|
|
339 M<-solve(chol(mixedmod$vg*K_norm+mixedmod$ve*diag(n)))
|
|
|
340 Y_t<-crossprod(M,Y)
|
|
|
341 cof_t<-crossprod(M,cof)
|
|
|
342 GLS_lm<-summary(lm(Y_t~0+cof_t))
|
|
|
343 Res_H0<-GLS_lm$residuals
|
|
|
344 Q_ <- qr.Q(qr(cof_t))
|
|
|
345 RSS<-list()
|
|
|
346 for (j in 1:(nbchunks-1)) {
|
|
|
347 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof)])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
|
|
|
348 RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
349 rm(X_t)}
|
|
|
350 X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof)])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof)-1))])
|
|
|
351 RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
|
|
|
352 rm(X_t,j)
|
|
|
353 RSSf<-unlist(RSS)
|
|
|
354 RSS_H0<-sum(Res_H0^2)
|
|
|
355 df2<-n-df1-ncol(cof)
|
|
|
356 Ftest<-(rep(RSS_H0,length(RSSf))/RSSf-1)*df2/df1
|
|
|
357 pval<-pf(Ftest,df1,df2,lower.tail=FALSE)
|
|
|
358 list(out=rbind(data.frame(SNP=colnames(cof)[-1],'pval'=GLS_lm$coef[2:(ncol(cof)),4]),
|
|
|
359 data.frame('SNP'=colnames(X)[-which(colnames(X) %in% colnames(cof))],'pval'=pval)),cof=colnames(cof)[-1])} else {cat('error \n')}}
|
|
|
360 opt_extBIC_out<-bestmodel_pvals(opt_extBIC)
|
|
|
361 opt_mbonf_out<-bestmodel_pvals(opt_mbonf)
|
|
|
362
|
|
|
363 list(step_table=fwdbwd_table,pval_step=pval_step,RSSout=plot_RSS,bonf_thresh=-log10(0.05/m),opt_extBIC=opt_extBIC_out,opt_mbonf=opt_mbonf_out)}
|
|
|
364
|
|
|
365 plot_step_table<-function(x,type){
|
|
|
366 if (type=='h2') {plot(x$step_table$step,x$step_table$h2,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='h2')
|
|
|
367 abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
|
|
|
368 else if (type=='maxpval'){plot(x$step_table$step,-log10(x$step_table$maxpval),type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='-log10(max_Pval)')
|
|
|
369 abline(h=x$bonf_thresh,lty=2)
|
|
|
370 abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
|
|
|
371 else if (type=='BIC'){plot(x$step_table$step,x$step_table$BIC,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='BIC')
|
|
|
372 abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
|
|
|
373 else if (type=='extBIC'){plot(x$step_table$step,x$step_table$extBIC,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='EBIC')
|
|
|
374 abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
|
|
|
375 else {cat('error! \n argument type must be one of h2, maxpval, BIC, extBIC')}}
|
|
|
376
|
|
|
377 plot_step_RSS<-function(x){
|
|
|
378 op<-par(mar=c(5, 5, 2, 2))
|
|
|
379 plot(0,0,xlim=c(0,nrow(x$RSSout)-1),ylim=c(0,1),xlab='step',ylab='%var',col=0)
|
|
|
380 polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,3],0,0), col='brown1', border=0)
|
|
|
381 polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,2],0,0), col='forestgreen', border=0)
|
|
|
382 polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,1],0,0), col='dodgerblue4', border=0)
|
|
|
383 abline(v=(nrow(x$step_table)/2-0.5),lty=2)
|
|
|
384 par(op)}
|
|
|
385
|
|
|
386 plot_GWAS<-function(x) {
|
|
|
387 output_<-x$out[order(x$out$Pos),]
|
|
|
388 output_ok<-output_[order(output_$Chr),]
|
|
|
389
|
|
|
390 maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200)))
|
|
|
391 plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2))
|
|
|
392 # plot_col<-c('blue','darkgreen','red','cyan','purple')
|
|
|
393 size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x
|
|
|
394 a<-rep(maxpos[1],size[1])
|
|
|
395 b<-rep(plot_col[1],size[1])
|
|
|
396 for (i in 2:max(unique(output_ok$Chr))){
|
|
|
397 a<-c(a,rep(maxpos[i],size[i]))
|
|
|
398 b<-c(b,rep(plot_col[i],size[i]))}
|
|
|
399
|
|
|
400 output_ok$xpos<-output_ok$Pos+a
|
|
|
401 output_ok$col<-b
|
|
|
402 output_ok$col[output_ok$SNP %in% x$cof]<-'red'
|
|
|
403
|
|
|
404 d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2
|
|
|
405
|
|
|
406 plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome')
|
|
|
407 axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))
|
|
|
408 abline(h=x$bonf_thresh,lty=3,col='black')
|
|
|
409
|
|
|
410
|
|
|
411 if (length(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]) > 0) {
|
|
|
412 text(output_ok$xpos[-log10(output_ok$pval) > x$bonf_thresh], -log10(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7)
|
|
|
413 legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
|
|
|
414 } else {
|
|
|
415 legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
|
|
|
416 }
|
|
|
417 }
|
|
|
418
|
|
|
419 plot_region<-function(x,chrom,pos1,pos2){
|
|
|
420 region<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)
|
|
|
421 region$col<- if (chrom %% 2 == 0) {'gray60'} else {'gray10'}
|
|
|
422 region$col[which(region$SNP %in% x$cof)]<-'red'
|
|
|
423 plot(region$Pos,-log10(region$pval),type='p',pch=20,main=paste('chromosome',chrom,sep=''),xlab='position (bp)',ylab='-log10(pval)',col=region$col,xlim=c(pos1,pos2))
|
|
|
424 abline(h=x$bonf_thresh,lty=3,col='black')}
|
|
|
425
|
|
|
426
|
|
|
427 plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {
|
|
|
428 stopifnot(step<=length(x$pval_step))
|
|
|
429 output<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by='SNP'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)
|
|
|
430 plot_GWAS(output)}
|
|
|
431
|
|
|
432 plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {
|
|
|
433 stopifnot(step<=length(x$pval_step))
|
|
|
434 output<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by='SNP'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)
|
|
|
435 plot_region(output,chrom,pos1,pos2)}
|
|
|
436
|
|
|
437
|
|
|
438 plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {
|
|
|
439 if (opt=='extBIC') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by='SNP'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)
|
|
|
440 plot_GWAS(output)}
|
|
|
441 else if (opt=='mbonf') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by='SNP'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)
|
|
|
442 plot_GWAS(output)}
|
|
|
443 else {cat('error! \n opt must be extBIC or mbonf')}}
|
|
|
444
|
|
|
445 plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {
|
|
|
446 if (opt=='extBIC') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by='SNP'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)
|
|
|
447 plot_region(output,chrom,pos1,pos2)}
|
|
|
448 else if (opt=='mbonf') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by='SNP'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)
|
|
|
449 plot_region(output,chrom,pos1,pos2)}
|
|
|
450 else {cat('error! \n opt must be extBIC or mbonf')}}
|
|
|
451
|
|
|
452
|
|
|
453 qqplot_fwd_GWAS<-function(x,nsteps){
|
|
|
454 stopifnot(nsteps<=length(x$pval_step))
|
|
|
455 e<--log10(ppoints(nrow(x$pval_step[[1]]$out)))
|
|
|
456 ostep<-list()
|
|
|
457 ostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))
|
|
|
458 for (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}
|
|
|
459
|
|
|
460 maxp<-ceiling(max(unlist(ostep)))
|
|
|
461
|
|
|
462 plot(e,ostep[[1]],type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))
|
|
|
463 abline(0,1,col="dark grey")
|
|
|
464
|
|
|
465 for (i in 2:nsteps) {
|
|
|
466 par(new=T)
|
|
|
467 plot(e,ostep[[i]],type='b',pch=20,cex=0.8,col=i,axes='F',xlab='',ylab='',xlim=c(0,max(e)+1),ylim=c(0,maxp))}
|
|
|
468 legend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),'cof',sep=' '))
|
|
|
469 }
|
|
|
470
|
|
|
471 qqplot_opt_GWAS<-function(x,opt){
|
|
|
472 if (opt=='extBIC') {
|
|
|
473 e<--log10(ppoints(nrow(x$opt_extBIC$out)))
|
|
|
474 o<--log10(sort(x$opt_extBIC$out$pval))
|
|
|
475 maxp<-ceiling(max(o))
|
|
|
476 plot(e,o,type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste('optimal model according to extBIC'))
|
|
|
477 abline(0,1,col="dark grey")}
|
|
|
478 else if (opt=='mbonf') {
|
|
|
479 e<--log10(ppoints(nrow(x$opt_mbonf$out)))
|
|
|
480 o<--log10(sort(x$opt_mbonf$out$pval))
|
|
|
481 maxp<-ceiling(max(o))
|
|
|
482 plot(e,o,type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste('optimal model according to mbonf'))
|
|
|
483 abline(0,1,col="dark grey")}
|
|
|
484 else {cat('error! \n opt must be extBIC or mbonf')}}
|
|
|
485
|
|
|
486
|
