## BTI Bioinformatic class 17th May 2016 # ## BELOW ALL THE DOCUMENTATION MLMM ON INPUTS, FUNCTIONS. ## TO GET STARTED GO LINE 90!!! ############################################################################################################################################## ###MLMM - Multi-Locus Mixed Model ###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH ####### # ##note: require EMMA #library(emma) #source('emma.r') # ##INPUT FILE DEFINITION # #PHENOTYPE - Y: a vector of length m, with names(Y)=individual names #GENOTYPE - X: a n by m matrix, where n=number of individuals, m=number of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names #KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names #each of these data being sorted in the same way, according to the individual name # ##FOR PLOTING THE GWAS RESULTS #SNP INFORMATION - snp_info: a data frame having at least 3 columns: # - 1 named 'SNP', with SNP names (same as colnames(X)), # - 1 named 'Chr', with the chromosome number to which belong each SNP # - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to. ####### # ##FUNCTIONS USE #save this file somewhere on your computer and source it! #source('path/mlmm.r') # ###FORWARD + BACKWARD ANALYSES #mygwas<-mlmm(Y,X,K,nbchunks,maxsteps) #X,Y,K as described above #nbchunks: an integer defining the number of chunks of X to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory #maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0, # 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. # It's value must be specified as an integer >= 3 # ###RESULTS # ##STEPWISE TABLE #mygwas$step_table # ##PLOTS # ##PLOTS FORM THE FORWARD TABLE #plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC')) # ##RSS PLOT #plot_step_RSS(mygwas) # ##GWAS MANHATTAN PLOTS # #FORWARD STEPS #plot_fwd_GWAS(mygwas,step,snp_info,pval_filt) #step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor) #snp_info as described above #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot # #OPTIMAL MODELS #Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria # #plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt) #snp_info as described above #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot # ##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST #plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2) #step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor) #snp_info as described above #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot #chrom is an integer specifying the chromosome on which the region of interest is #pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info # #plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2) #snp_info as described above #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot #chrom is an integer specifying the chromosome on which the region of interest is #pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info # ##QQPLOTS of pvalues #qqplot_fwd_GWAS(mygwas,nsteps) #nsteps=maximum number of forward steps to be displayed # #qqplot_opt_GWAS(mygwas,opt=c('extBIC','mbonf')) # ############################################################################################################################################## ## define path and load libraries## ################################### setwd('/bioinfo/Data/gwas') source("mlmm_cof.r") source("mlmm.r") source("emma.r") library(matrixcalc) library(corpcor) ## clear cache (if needed) rm(list=ls()) ###################### ## load input files ## ###################### # 1 # load list of SNP, chromosome and position# ################################################ snp_info = read.table("GWAS_snp_info.txt", sep="\t", header=T) # 2 # load genotypic data# ########################## X = read.table("GWAS_genotype.txt", sep=",", header = T, stringsAsFactors = F) # check dimensions and set as matrix dim(X) X <- as.matrix(X) class(X) # 3a # Relationship matrix, option 1: load pre-computed kinship matrix# ####################################################################### K = read.table("GWAS_kinship.txt", sep="\t", header = T, stringsAsFactors = F) dim(K) K.mod = as.matrix(K) # check if it is positive definite class(K.mod) is.positive.definite(K.mod) # if "FALSE" , make it: K.mod = make.positive.definite(K) is.positive.definite(K.mod) # 3b # Relationship matrix, option 2: compute it ################################################ X_ok <- 2*X Allele_p_freq <- colSums(X_ok, na.rm = T) / apply(X_ok, 2, function(x) {2 * length(which(!is.na(x)))}) Allele_q_freq <- 1 - Allele_p_freq X_cent <- sweep(X_ok, 2, 2*Allele_p_freq, "-") X_cent = as.matrix(X_cent) G_mat <- tcrossprod(X_cent) / (2 * sum(Allele_p_freq * Allele_q_freq)) # clean intermediate files rm(Allele_p_freq, Allele_q_freq, X2, X_ok, X_cent) # save the Relationship matrix write.table(G_mat, "Gmat_GWAS.txt", sep= "\t") # 4 # load phenotypic data ########################## Y_ok = read.table("GWAS_phenotypes.txt", sep="\t", header = T, stringsAsFactors = F) dim(Y_ok) # Display phenotype list names(Y_ok) # select a phenotypic trait for GWAS (here Malate) Y = as.numeric(Y_ok[,22]) head(Y) # 5 # load structure covariate (only if fixed effect (Q, K or SNP) included in the model)# ############################### Q = read.table("GWAS_structureQ1.txt", sep="\t", header = T, stringsAsFactors = F) Q = as.matrix(Q[,1]) ########################## ### RUN THE MIXED MODEL### ########################## # RUNNING MLMM WITHOUT COFACTOR# ################################ # with precompute kinship (option 3a) mygwas <-mlmm(Y, X, K.mod,5,10) # with computed kinship (option 3b) mygwas_G_mat <-mlmm(Y, X, G_mat,5,10) # view content of my gwas object names(mygwas) # RUNNING MLMM WITH COFACTOR# ############################# mygwas.FW.Q = mlmm_cof(Y, X, Q, G_mat, 2, 10) mygwas.FW.Q <- mlmm_cof(Y, X, Q, K.mod,5,10) ################# ##PLOTTING DATA## ################# # MANHATTAN PLOT FOR THE BEST MODEL (MBONF) plot_opt_GWAS(mygwas, opt = c('mbonf'), snp_info, 1) # COMPILE TRAIT DISTRIBUTION AND MANHATTAN FOR THE BEST MODELS (EBIC and MBONF) par(mfrow = c(3,1)) hist(Y, main = "Your trait", cex.axis=2) plot_opt_GWAS(mygwas, opt = c('mbonf'), snp_info, 1) plot_opt_GWAS(mygwas, opt = c('extBIC'), snp_info, 1) # COMPILE TRAIT DISTRIBUTION AND MANHATTAN FOR THE BEST MODELS WITH AND WITHOUT COFACTOR (MBONF) par(mfrow = c(3,1)) hist(Y, main = "Your trait", cex.axis=2) plot_opt_GWAS(mygwas, opt = c('mbonf'), snp_info, 1) plot_opt_GWAS(mygwas.FW.Q, opt = c('mbonf'), snp_info, 1) # PLOTTING THE LIKELIHOOD ACCORDING TO FOUR MODELS AND FOR EVERY STEP par(mfrow = c(2,2)) plot_step_table(mygwas,'h2') plot_step_table(mygwas,'maxpval') plot_step_table(mygwas,'BIC') plot_step_table(mygwas,'extBIC') # PLOTTING THE GENETIC VARIANCE EVOLUTION AT EVERY STEP plot_step_RSS(mygwas) plot_step_RSS(mygwas.FW.Q) # PLOT THE QQPLOT OF ALL STEPS IN A GRAPH qqplot_fwd_GWAS(mygwas,7) # PLOT THE MANAHATTAN PLOT FOR EACH STEP par(mfrow = c(2,3)) plot_fwd_GWAS(mygwas,1,snp_info,1) plot_fwd_GWAS(mygwas,2,snp_info,1) plot_fwd_GWAS(mygwas,3,snp_info,1) plot_fwd_GWAS(mygwas,4,snp_info,1) plot_fwd_GWAS(mygwas,5,snp_info,1) plot_fwd_GWAS(mygwas,6,snp_info,1) plot_fwd_GWAS(mygwas,7,snp_info,1) plot_fwd_GWAS(mygwas,8,snp_info,1) plot_fwd_GWAS(mygwas,9,snp_info,1) plot_fwd_GWAS(mygwas,10,snp_info,1) # PLOT THE QQPLOT DEPENDING OF THE CORRECTION (BEST MODEL) par(mfrow = c(2,1)) qqplot_opt_GWAS(mygwas,'extBIC') qqplot_opt_GWAS(mygwas,'mbonf') # PLOT PHYSICAL REGION plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2) #snp_info as described above #pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot #chrom is an integer specifying the chromosome on which the region of interest is #pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info plot_opt_region(mygwas,opt=c('mbonf'),snp_info, 1 ,6 ,40000000,43000000) ### RUN MLMM ON ALL TRAITS AT ONCE ### ###################################### Y_all <- Y_ok[, 2:ncol(Y_ok)] rownames(Y_ok) <- Y_ok[, 1] mygwas_all <- NULL mygwas_all <- list() for (i in 1:ncol(Y_all)) { mygwas_all[[i]] <- mlmm(Y_all[,i], X, G_mat, 2, 10) } # WRITE OUTPUT FILES FROM MLMM LOOP save(mygwas_all, file="mygwas_agro_2011only.r") # EXTRACT THE RESULTS (DEPENDS ON THE BEST METHOD OF CORRECTION MBONF OR EXTENDED BIC) ###################################################################################### ## EXTENDED BIC mygwas$opt_extBIC$cof res <- mygwas$opt_extBIC$out res.sort <- res[with(res, order(pval,SNP)),] head(res.sort) write.table(res.sort, "MLMM_YOUR_TRAIT_eBIC.txt", sep= "\t") ## MBONF mygwas$opt_mbonf$cof res <- mygwas$opt_mbonf$out res.sort <- res[with(res, order(pval,SNP)),] head(res.sort) write.table(res.sort, "MLMM_YOUR_TRAIT__mBonf.txt", sep= "\t") ########################################## ### EXTRACT DATA ON mygwas_all R OBJECT### ########################################## setwd('\\\\pac-sm-nas01/profils$/pc/gbauchet/Bureau/R_code_MLMM_cofactor/') source("mygwas_allmetabo_2012.r") mygwas_all = read.table("mygwas_agro_2011only.r") # DESCRIPTION OF "mygwas_all" R OBJECT CONTENT: # cofactors steps 1-20 #(mygwas_all[[i]][1]): [1] "step_table" $step $step_ $cof $ncof $h2 $maxpval $BIC $extBIC ##pvalues and estimates #(mygwas_all[[i]][2]): [1] "pval_step" $out $cof $coef -> Estimate Std. Error t value ##Heritabilities #(mygwas_all[[i]][3]): [1] "RSSout" expl_RSS h2_RSS unexpl_RSS ##Bonferronni threshold #(mygwas_all[[i]][4]): [1] "bonf_thresh" a single value #Extended BIC values, threshold and estimate #(mygwas_all[[i]][5]): [1] "opt_extBIC" $out $cof $coef #mBonferonni values, threshold and estimate #(mygwas_all[[i]][6]): [1] "opt_mbonf" (names(mygwas_all[[60]][6]$opt_mbonf)) $out $cof $coef ##COFACTOR STEPS DATA ###################### i = NULL step_out.tot = NULL step__out.tot = NULL cof_out.tot = NULL ncof_out.tot = NULL h2_out.tot = NULL maxpval_out.tot = NULL BIC_out.tot = NULL extBIC_out.tot = NULL for (i in 1:length(mygwas_all)) { z <- mygwas_all[[i]] #extract steps # step_out <- z$step_table[,1] length(step_out) <- 20 setp_out.tot <- cbind(step_out.tot, step_out) #extract steps fwd or bwd position step__out <- z$step_table[,2] length(step__out) <- 20 step__out.tot <- cbind(step__out.tot, step__out) #extract cof snp id cof_out <- as.vector(z$step_table[,3]) length(cof_out) <- 20 cof_out.tot <- cbind(cof_out.tot, cof_out) #extract cof setp # ncof_out <- z$step_table[,4] length(ncof_out) <- 20 ncof_out.tot <- cbind(ncof_out.tot, ncof_out) #extract h? for each step h2_out <- z$step_table[,5] length(h2_out) <- 20 h2_out.tot <- cbind(h2_out.tot, h2_out) #extract maxpval at each step maxpval_out <- z$step_table[,6] length(maxpval_out) <- 20 maxpval_out.tot <- cbind(maxpval_out.tot, maxpval_out) #extract BIC at each step BIC_out <- z$step_table[,7] length(BIC_out) <- 20 BIC_out.tot <- cbind(BIC_out.tot, BIC_out) #extract eBIC at each step extBIC_out <- z$step_table[,8] length(extBIC_out) <- 20 extBIC_out.tot <- cbind(extBIC_out.tot, extBIC_out) } write.table(step_out.tot, "mygwas_all_step_out.txt", sep= "\t") write.table(step__out.tot, "mygwas_step_out.txt", sep= "\t") write.table(cof_out.tot, "mygwas_cof_out.txt", sep= "\t") write.table(ncof_out.tot, "mygwas_ncof_out.txt", sep= "\t") write.table(h2_out.tot, "mygwas_h2_out.txt", sep= "\t") write.table(maxpval_out.tot, "mygwas_maxpval.txt", sep= "\t") write.table(BIC_out.tot, "mygwas_BIC.txt", sep= "\t") write.table(extBIC_out.tot, "mygwas_extBIC.txt", sep= "\t") ##HERITABILITIES ################# i = NULL expl_Rss_out.tot = NULL h2_Rss_out.tot = NULL unexpl_Rss_out.tot = NULL for (i in 1:length(mygwas_all)) { z <- mygwas_all[[i]] #extract explained variance expl_Rss_out <- z$RSSout[,1] length(expl_Rss_out) <- 20 expl_Rss_out.tot <- cbind(expl_Rss_out.tot, expl_Rss_out) #extract h? values h2_Rss_out <- z$RSSout[,2] length(h2_Rss_out) <- 20 h2_Rss_out.tot <- cbind(h2_Rss_out.tot, h2_Rss_out) #extract unexplained variance unexpl_Rss_out <- z$RSSout[,3] length(unexpl_Rss_out) <- 20 unexpl_Rss_out.tot <- cbind(unexpl_Rss_out.tot, unexpl_Rss_out) } write.table(expl_Rss_out.tot, "mygwas_all_expl_RSS.txt", sep= "\t") write.table(h2_Rss_out.tot, "mygwas_all_RSS.txt", sep= "\t") write.table(unexpl_Rss_out.tot, "mygwas_all_unexpl_RSS.txt", sep= "\t") ##BONFERONNI THRESHOLD ###################### write.table(mygwas_all$bonf_thresh, "mygwas_all_bonf_tresh.txt", sep= "\t") ##EXTENDED eBIC VALUE THRESHOLD AND ESTIMATE ########################################### i = NULL eBIC_out_out.tot = NULL eBIC_cof_out.tot = NULL eBIC_coef_out.tot = NULL eBIC_out_out.tot <- 9436 for (i in 1:length(mygwas_all)) { z <- mygwas_all[[i]] #extract pvalues with BIC correction eBIC_out_out <- z$opt_extBIC$out #sorting (need to add snp id for each column) #BIC_out_out_sort <- eBIC_out_out[with(BIC_out_out, order(pval,SNP)),] #length(eBIC_out_out.tot) <- 9381 eBIC_out_out.tot <- cbind(eBIC_out_out.tot, eBIC_out_out[-2], eBIC_out_out[-1]) #extract cofactor SNP ID with eBIC correction eBIC_cof_out <- z$opt_extBIC$cof length(eBIC_cof_out) <- 10 eBIC_cof_out.tot <- cbind(eBIC_cof_out.tot, eBIC_cof_out) #extract estimate and std error for each snp used as cofactors with eBIC correction eBIC_coef_out <- z$opt_extBIC$coef eBIC_coef_out.tot <- rbind(eBIC_coef_out.tot, eBIC_coef_out[1,]) } write.table(eBIC_out_out.tot, "mygwas_all_eBIC_out.txt", sep= "\t") write.table(eBIC_cof_out.tot, "mygwas_all_eBIC_cof.txt", sep= "\t") write.table(eBIC_coef_out.tot, "mygwas_all_eBIC_coef.txt", sep= "\t") ##MBONFERONNI VALUE THRESHOLD AND ESTIMATE ########################################## i = NULL mBonf_out_out.tot = NULL mBonf_cof_out.tot = NULL mBonf_coef_out.tot = NULL mBonf_out_out.tot <- 9436 for (i in 1:length(mygwas_all)) { z <- mygwas_all[[i]] #extract pvalues with mBonf correction mBonf_out_out <- z$opt_mbonf$out #sorting (need to add snp id for each column) #mBonf_out_out_sort <- mBonf_out_out[with(mBonf_out_out, order(pval,SNP)),] #length(mBonf_out_out.tot) <- 9381 #mBonf_out_out.tot <- cbind(mBonf_out_out.tot,rownames(mBonf_out_out)) #mBonf_out_out.tot <- cbind(rn = rownames(mBonf_out_out), mBonf_out_out[1], stack(mBonf_out_out[-1]), row.names = NULL) #mBonf_out_out.tot <- cbind(mBonf_out_out.tot, mBonf_out_out$SNP, mBonf_out_out$pval) mBonf_out_out.tot <- cbind(mBonf_out_out.tot, mBonf_out_out[-2], mBonf_out_out[-1]) #extract cofactor SNP ID with mBonf correction mBonf_cof_out <- z$opt_mbonf$cof length(mBonf_cof_out) <- 10 mBonf_cof_out.tot <- rbind(mBonf_cof_out.tot, mBonf_cof_out) #extract estimate and std error for each snp used as cofactors with mBonf correction mBonf_coef_out <- z$opt_mbonf$coef mBonf_coef_out.tot <- rbind(mBonf_coef_out.tot, mBonf_coef_out[1,]) } write.table(mBonf_out_out.tot, "mygwas_all_mBonf_out.txt", sep= "\t") write.table(mBonf_cof_out.tot, "mygwas_all_mBonf_cof.txt", sep= "\t") write.table(mBonf_coef_out.tot, "mygwas_all_mBonf_coef.txt", sep= "\t")