pheno_geno_corr

email title: 'help with Miami plots' Carolina (Aug 4th, 2021)

Aug 13th 2021

Location on server: /space/chen-syn01/1/data/cinliu/plots/pheno_geno_corr

Goal plot:

Original instructions:

The second figure is the phenotypic and genetic correlation matrix. You can find the script, input files and output in this dropbox folder - I just tested it so you should be able to run the R script easily once you just update your path on lines 24 and 116. The main issue with the figure is the diagonals - it holds statistics for the genetic correlations per region between GWASr and GWASr+g analyses, and that information is getting lost in the current figure. We discussed putting these values of the diagonals in a separate column altogether to the side of the matrix (this would require modification of line 108 of the code, where we could just have empty values along the diagonal instead). If you want to play around with anything manually though, there is also a ppt file in that folder where I have added the actual correlation values along the diagonal. All of these values have a highly significant p-value- so, if we were to put a box or something like that around sig values, they would all be "boxed". Let's prioritize putting these values as a separate column, but if you have any other ideas for how to best represent this information, I trust your judgment!

globalanalysis_pheno_geno_corrmat.R

####creating pheno/geno correlation matrix between gclust phenotypes##pheno_cor: lower triangle is pheno corr with globals regressed out##pheno2_cor: upper right triangle is pheno corr NOT regressing out globals##diagonals hold geno correlations (rg), although we might alter so this is in a separate column to the side#load packagesrm(list=ls())
library(psych)library(tableone)library(tidyverse)library(Hmisc)library(pastecs)library(data.table)library(ggplot2)library(gplots)library(RColorBrewer)library(Matrix)library(plyr)library(dplyr)
#part 1: pheno_corr, load already generated phenotypic correlation matrix for phenotypes pre-residualized for globals
###-----update path------dir="/Users/nini/Desktop/2021lab/carolina/pheno_geno_corr/input_files/"###-------------------
pheno_cor=read.table(paste0(dir,'UKB_pheno_matrix.txt'),h=T, sep='\t', row.names = 1)
#reorder based on Peng et al 2016 paper and Makowski et al Science 2021
reorder<-c("motorpremotor_area","dorsolateralprefrontal_area","dorsomedialfrontal_area","orbitofrontal_area", "parsopercularis_area","superiortemporal_area","posterolateraltemporal_area","anteromedialtemporal_area", "inferiorparietal_area","superiorparietal_area","precuneus_area","occipital_area","motor_premotor_SMA_thickness", "superiorparietal_thickness","inferiorparietal_thickness","perisylvian_thickness","occipital_thickness", "ventromedialoccipital_thickness","ventralfrontal_thickness","temporalpole_thickness","medialtemporal_thickness" , "middletemporal_thickness" , "dorsolateralprefrontal_thickness","medialprefrontal_thickness")
pheno_cor_reorder<-pheno_cor[reorder,reorder]
pheno_cor_reorder_number = pheno_cor_reorderrow.names(pheno_cor_reorder_number)<-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area", "5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area", "9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness", "2.superiorparietal_thickness","3.inferiorparietal_thickness","4.perisylvian_thickness","5.occipital_thickness", "6.ventromedialoccipital_thickness","7.ventralfrontal_thickness","8.temporalpole_thickness","9.medialtemporal_thickness" , "10.middletemporal_thickness" , "11.dorsolateralprefrontal_thickness","12.medialprefrontal_thickness")
names(pheno_cor_reorder_number)<-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area", "5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area", "9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness", "2.superiorparietal_thickness","3.inferiorparietal_thickness","4.perisylvian_thickness","5.occipital_thickness", "6.ventromedialoccipital_thickness","7.ventralfrontal_thickness","8.temporalpole_thickness","9.medialtemporal_thickness" , "10.middletemporal_thickness" , "11.dorsolateralprefrontal_thickness","12.medialprefrontal_thickness")

#part 2: pheno2_corr, load already generated phenotypic correlation matrix for phenotypes NOT pre-residualized for globals
pheno2_cor=read.table(paste0(dir,'UKB39k_24pheno_presidwithoutglobals.txt'),h=T, sep='\t', row.names = 1)pheno2_cor_reorder<-pheno2_cor[reorder,reorder]
pheno2_cor_reorder_number = pheno2_cor_reorderrow.names(pheno2_cor_reorder_number)<-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area", "5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area", "9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness", "2.superiorparietal_thickness","3.inferiorparietal_thickness","4.perisylvian_thickness","5.occipital_thickness", "6.ventromedialoccipital_thickness","7.ventralfrontal_thickness","8.temporalpole_thickness","9.medialtemporal_thickness" , "10.middletemporal_thickness" , "11.dorsolateralprefrontal_thickness","12.medialprefrontal_thickness")
names(pheno2_cor_reorder_number)<-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area", "5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area", "9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness", "2.superiorparietal_thickness","3.inferiorparietal_thickness","4.perisylvian_thickness","5.occipital_thickness", "6.ventromedialoccipital_thickness","7.ventralfrontal_thickness","8.temporalpole_thickness","9.medialtemporal_thickness" , "10.middletemporal_thickness" , "11.dorsolateralprefrontal_thickness","12.medialprefrontal_thickness")
##########################combining the two pheno matrices into one, then will add in diagonal info separately##########################lowertri will be phenotypic correlation with globals, upper tri will be without globals#diagonal will be genetic correlation for each phenotype with and without globalsa<-tril(as.matrix(pheno_cor_reorder_number),-1) #lower triangular matrix, omit diagonalsb<-triu(as.matrix(pheno2_cor_reorder_number),1) #upper triangular matrixcc<-a+bcc<-as.data.frame(as.matrix(cc))
## re-order again based on hclust results of reordered dendrogram based on Makowski et al 2020
rereorder <-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area", "5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area", "9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness", "2.superiorparietal_thickness","3.inferiorparietal_thickness","5.occipital_thickness","4.perisylvian_thickness","8.temporalpole_thickness", "9.medialtemporal_thickness" ,"10.middletemporal_thickness" , "6.ventromedialoccipital_thickness", "11.dorsolateralprefrontal_thickness","7.ventralfrontal_thickness","12.medialprefrontal_thickness")
cc_reorder <- cc[rereorder,rereorder]

#load in rg_globals diagonal in order aboverg_diag=read.table(paste0(dir,'diagonals_rg.txt'),h=T, sep='\t', row.names = 1)rereorder2 <-c("motorpremotor_area","dorsolateralprefrontal_area","dorsomedialfrontal_area","orbitofrontal_area", "parsopercularis_area","superiortemporal_area","posterolateraltemporal_area","anteromedialtemporal_area", "inferiorparietal_area","superiorparietal_area","precuneus_area","occipital_area","motor_premotor_SMA_thickness", "superiorparietal_thickness","inferiorparietal_thickness","occipital_thickness","perisylvian_thickness","temporalpole_thickness", "medialtemporal_thickness" ,"middletemporal_thickness" , "ventromedialoccipital_thickness", "dorsolateralprefrontal_thickness","ventralfrontal_thickness","medialprefrontal_thickness")rg_diag_reorder=(rg_diag[rereorder2,])
#add genetic correlations for diagonaldiag(cc_reorder) <- rg_diag_reorder$rg
write.table(cc_reorder,file=paste0(dir,'pheno_matrix_reorder.txt'), quote=F)

#new correlation plot
####---update path---setwd("/Users/nini/Desktop/2021lab/carolina/pheno_geno_corr/")###--------------
library(RColorBrewer)library("gplots")
df = read.table(paste0(dir,'pheno_matrix_reorder.txt'), h=T)dfMatrix <- data.matrix(df)
col1<-colorRampPalette(c("blue2","white","magenta2"))
#color side bar 1 png(paste0("1Peng_reordered_globalslowtri_noglobalsuptri_rgdiag_matrix.png"),width = 15*300, height = 12*300, res = 300, pointsize = 10)par(oma=c(2,0,0,20),xpd=TRUE)# xpd allows Draw outside plot areaheatmap.2(dfMatrix,dendrogram ='none', trace = "none",density.info = "none", Colv = NA, Rowv = NA, labCol = FALSE,cexRow =3, offsetRow=11,breaks=seq(-1, 1, length.out=200.1), key.title="",key.xlab="", key.par=list(mar=c(4,0,1,2), cex=1.0, cex.lab=2.0, cex.axis=2.0), col= col1(200), RowSideColors = c( # grouping row-variables into different colors rep("paleturquoise4", 12), # area darker turquoise rep("paleturquoise2", 12)), #thickness lmat=rbind(c(3,0,0,4,5), c(0,0,2,1,0)), #specifying the bar to be on the right lhei=c(0.5,4), lwid=c(0.1,0.1,4,0.25,1))dev.off()
#####color side bar 2png(paste0("2Peng_reordered_globalslowtri_noglobalsuptri_rgdiag_matrix.png"),width = 15*300, height = 12*300, res = 300, pointsize = 10)par(oma=c(2,0,0,20),xpd=TRUE)# xpd allows Draw outside plot areaheatmap.2(dfMatrix,dendrogram ='none', trace = "none",density.info = "none", Colv = NA, Rowv = NA, labCol = FALSE,cexRow =3, offsetRow=11,breaks=seq(-1, 1, length.out=200.1), key.title="",key.xlab="", key.par=list(mar=c(4,0,1,2), cex=1.0, cex.lab=2.0, cex.axis=2.0), col= col1(200), RowSideColors = c( # grouping row-variables into different colors rep("blue3", 6), #dark blue rep("#ef062d", 6), #red rep("blue", 5), #blue rep("#e131b5", 3), #pink rep("blue", 1), #blue rep("#e131b5", 3)), #pink lmat=rbind(c(3,0,0,4,5), c(0,0,2,1,0)), lhei=c(0.5,4), lwid=c(0.1,0.1,4,0.25,1))dev.off()

#####make individual plot for diaginal data - genetic correlations per region between GWASr and GWASr+g analyses
m <- as.data.frame(matrix(0, ncol = 1, nrow = 24))####make plot for middle value only for (i in 1:24) { m[i,] <- dfMatrix[i,i]}m$order <- letters[1:24]mMatrix <- as.matrix(m)m$name <- colnames(df)
q <- ggplot(m, aes(x = order,y = 1, fill = V1)) + geom_tile() + scale_fill_gradientn(name="",limits=c(0.25,0.75), breaks =c(0.3,0.5,0.7),colours=c("deepskyblue1","dodgerblue4")) + theme(legend.position = "top", legend.key.height = unit(0.9, 'cm'), #change legend key height legend.key.width = unit(0.7, 'cm'), legend.text=element_text(size=14))
ggsave(paste0(dir,"bluesingle_column.png"),plot = last_plot(), width = 30,height = 7,units ="cm",dpi = 300)
#####Make diaginal blank plot for (i in 1:24) { dfMatrix[i,i] <- 0 }png(paste0("3Peng_reordered_globalslowtri_noglobalsuptri_rgdiag_matrix.png"),width = 15*300, height = 12*300, res = 300, pointsize = 10)par(oma=c(2,0,0,20),xpd=TRUE)# xpd allows Draw outside plot areaheatmap.2(dfMatrix,dendrogram ='none', trace = "none",density.info = "none", Colv = NA, Rowv = NA, labCol = FALSE,cexRow =3, offsetRow=11,breaks=seq(-1, 1, length.out=200.1), key.title="",key.xlab="", key.par=list(mar=c(4,0,1,2), cex=1.0, cex.lab=2.0, cex.axis=2.0), col= col1(200), RowSideColors = c( # grouping row-variables into different colors rep("paleturquoise4", 12), # area darker turquoise rep("paleturquoise2", 12)), #thickness lmat=rbind(c(3,0,0,4,5), c(0,0,2,1,0)), #specifying the bar to be on the right lhei=c(0.5,4), lwid=c(0.1,0.1,4,0.25,1))dev.off()

Plot put together in powerpoint: /space/chen-syn01/1/data/cinliu/plots/pheno_geno_corr/phencorr_rg_matrix.pptx

Page updated

Report abuse