Correlation between 2 files

Slack CH April 8th 2021 - April 17th 2021

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

Original instructions

Hi Nini, Can you help calculating correlation to evaluate consistency between two files?

"/space/chen-syn01/1/data/chen/brain_age/00README"

[cinliu@ip99 data]$ cd /space/chen-syn01/1/data/chen/brain_age/

[cinliu@ip99 brain_age]$ cat 00README

# copy Hao's xgboost results using a symbolic link

ln -s /space/chen-syn01/1/data/haowang/ukb_age/phenotype/UKB /space/chen-syn01/1/data/chen/brain_age/xgboost

# I copied Tuo's age trajectory results from dropbox and saved in the "residuals" subdirectory.

# evaluate consistency of the above two results by calculating correlation for each region

# For example, for this region, occipital_thickness, you can find the values here

#/space/chen-syn01/1/data/haowang/ukb_age/phenotype/UKB/UKB_gclus_xgb_result_female.txt # the variable is ChenT_occipital_thickness_brainAGER_Female # look for brainAGER

# /space/chen-syn01/1/data/chen/brain_age/residuals/Thickness_nonadjust_female.csv # the variable is ChenT_bh_occipital_thickness # look for nonadjust and sex stratified spreadsheets.

# Hao's results only have the testing sample (N=1931) (training/testing: 90/10 split), a subset of Tuo's sample

Results in corplot form:

consistency_eval_female.R

rm(list=ls())#consistency_evallibrary(corrplot)library(RColorBrewer)
setwd("/Users/nini/Desktop/2021lab/CH/consistency_eval")
####FEMALE AREA#######dat1_ukb <- read.delim("~/Desktop/2021lab/CH/consistency_eval/UKB/UKB_gclus_xgb_result_female.txt")dat1_area <- dat1_ukb[,163:174]#163-174 area "ChenA_region_area_brainAGER_Female"
#colnames(dat1_area) <- gsub("_brainAGER_Female","",colnames(dat1_area))#colnames(dat1_area) <- gsub("ChenA_","",colnames(dat1_area))dat1_area_cor <- cor(dat1_area)

###dat2 dat2 <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Area_nonadjust_female.csv")dat2_area <- dat2[,3:14]colnames(dat2_area) <- gsub("ChenA_bh_", "", colnames(dat2_area))dat2_area_cor <- cor(dat2_area)
#################two plot in same page jpeg(height=1800, width=3200, pointsize = 15,"image/Female_area.jpg")par(mfrow=c(1,2))corrplot(dat1_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Hao's data\nFile: UKB_gclus_xgb_result_female.txt\nRegion name simplified from: ChenA_region_area_brainAGER_Female", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))corrplot(dat2_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Tuo's data\nFile: Area_nonadjust_female.csv\nRegion name simplified from: ChenA_bh_region_area", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))dev.off()
#####FEMALE THICKNESS######
dat1_thickness <- dat1_ukb[,175:186]#175:186 thickness "ChenT_region_thickness_brainAGER_Female"colnames(dat1_thickness) <- gsub("_brainAGER_Female","",colnames(dat1_thickness))colnames(dat1_thickness) <- gsub("ChenT_","",colnames(dat1_thickness))dat1_area_cor <- cor(dat1_thickness)
###dat2 dat2 <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Thickness_nonadjust_female.csv")dat2_thickness <- dat2[,3:14]colnames(dat2_thickness) <- gsub("ChenT_bh_", "", colnames(dat2_thickness))dat2_area_cor <- cor(dat2_thickness)
#################two plot in same page jpeg(height=1800, width=3200, pointsize = 15,"image/Female_thickness.jpg")par(mfrow=c(1,2))corrplot(dat1_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Hao's data\nFile: UKB_gclus_xgb_result_female.txt\nRegion name simplified from: ChenT_region_thickness_brainAGER_Female", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))corrplot(dat2_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Tuo's data\nFile: Thickness_nonadjust_female.csv\nRegion name simplified from: ChenT_bh_region_area", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))dev.off()

consistency_eval_male.R

rm(list=ls())#consistency_evallibrary(corrplot)library(RColorBrewer)
setwd("/Users/nini/Desktop/2021lab/CH/consistency_eval")
####Male AREA#######dat1_ukb <- read.delim("~/Desktop/2021lab/CH/consistency_eval/UKB/UKB_gclus_xgb_result_male.txt")dat1_area <- dat1_ukb[,163:174]#163-174 area "ChenA_region_area_brainAGER_Male"colnames(dat1_area) <- gsub("_brainAGER_Male","",colnames(dat1_area))colnames(dat1_area) <- gsub("ChenA_","",colnames(dat1_area))dat1_area_cor <- cor(dat1_area)
###dat2 dat2 <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Area_nonadjust_male.csv")dat2_area <- dat2[,3:14]colnames(dat2_area) <- gsub("ChenA_bh_", "", colnames(dat2_area))dat2_area_cor <- cor(dat2_area)
#################two plot in same page jpeg(height=1800, width=3200, pointsize = 15,"image/Male_area.jpg")par(mfrow=c(1,2))corrplot(dat1_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Hao's data\nFile: UKB_gclus_xgb_result_male.txt\nRegion name simplified from: ChenA_region_area_brainAGER_Male", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))corrplot(dat2_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Tuo's data\nFile: Area_nonadjust_male.csv\nRegion name simplified from: ChenA_bh_region_area", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))dev.off()
#####Male THICKNESS######
dat1_thickness <- dat1_ukb[,175:186]#175:186 thickness "ChenT_region_thickness_brainAGER_Male"colnames(dat1_thickness) <- gsub("_brainAGER_Male","",colnames(dat1_thickness))colnames(dat1_thickness) <- gsub("ChenT_","",colnames(dat1_thickness))dat1_area_cor <- cor(dat1_thickness)
###dat2 dat2 <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Thickness_nonadjust_male.csv")dat2_thickness <- dat2[,3:14]colnames(dat2_thickness) <- gsub("ChenT_bh_", "", colnames(dat2_thickness))dat2_area_cor <- cor(dat2_thickness)
#################two plot in same page jpeg(height=1800, width=3200, pointsize = 15,"image/Male_thickness.jpg")par(mfrow=c(1,2))corrplot(dat1_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Hao's data\nFile: UKB_gclus_xgb_result_male.txt\nRegion name simplified from: ChenT_region_thickness_brainAGER_Male", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))corrplot(dat2_area_cor, order="alphabet", type = "upper",tl.cex = 2,method = "color", addCoef.col = "black",number.cex=2, title = "Tuo's data\nFile: Thickness_nonadjust_male.csv\nRegion name simplified from: ChenT_bh_region_area", mar=c(0,0,4,0), col=brewer.pal(n=8, name="RdYlBu"))dev.off()

Other things tried

looking at correlation with Rmarkdown

correlation.Rmd

---title: "Consistency Check"output: html_document---Pearson's corrlation (cor.test(method = pearson))```{r setup, include=FALSE}knitr::opts_chunk$set(echo = TRUE)```## AREA## Female Area TableTuo's data = residuals/Area_nonadjust_female.csv ; Hao's data = UKB_gclus_xgb_result_female.txt
ChenA_bh_region_area v.s ChenA_region_area_brainAGER_Female ```{r female area, include=FALSE}#female arearemove(list=ls())library(data.table)library(dplyr)library(formattable)library(tidyr)#data from Tuo Area_nonadjust_female <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Area_nonadjust_female.csv")FID_female <- read.csv("~/Desktop/2021lab/CH/consistency_eval/Dropbox_GWASstudy/residuals/FID_female.csv")Tuo_dat <- cbind(FID_female[,2],Area_nonadjust_female[,3:14]) #Tuo female thickness colnames(Tuo_dat)[1] <- "FID"###hao dat UKB_gclus_xgb_result_female <- read.delim("~/Desktop/2021lab/CH/consistency_eval/UKB/UKB_gclus_xgb_result_female.txt")Hao_dat <- UKB_gclus_xgb_result_female[,c(1,163:174)] #163:174 area AGER, 175:186 thickness AGER#colnames(Hao_dat)###reorder rows/make FID match Tuo_dat$order <- match(Tuo_dat$FID,Hao_dat$FID)Tuo_dat<- Tuo_dat[!(is.na(Tuo_dat$order)),] #since hao's data is a subset of tuo's we need to remove the data that Hao doesn't have##double check ## identical(Tuo_dat$FID,Hao_dat$FID) == TRUE ###reorder column orders list <- gsub("ChenA_bh_","",colnames(Tuo_dat))num <- c(1:13)for (i in 1:13) { num[i] <- grep(list[i],colnames(Hao_dat)) #warning will appear about length,it's ok #print(num)} Hao_dat <- Hao_dat[,num]
#correlation calculation results = list()counter = 0df <- as.data.frame(list)
k = 1for (i in 1:13) { a <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") # print(paste(colnames(Hao_dat)[i], " est:", a$estimate, " p=value:", a$p.value)) df$cor[i] <- round(a$estimate,5) df$pval[i] <- a$p.value df$t[i] <- a$statistic df$df[i] <- a$parameter df$CI[i] <- paste(round(a$conf.int[1],5),round(a$conf.int[2],5)) results[[i]] <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") #print(results[[i]]) }
df$p_details <- "p-value < 2.2e-16" #manually adding this note with quick visual inspection from print(results[[i]]) df <- df[c(2:13),]
``````{r pressure, echo=FALSE}formattable(df, list( cor = color_tile("lightpink", "transparent") ))``````{r correlation summary1}summary(df$cor) #female area correlation summary```## Male Area TableTuo data = residuals/Area_nonadjust_male.csv;Hao data = UKB_gclus_xgb_result_male.txt
ChenA_bh_region_area v.s ChenA_region_area_brainAGER_Male ```{r male_area, include=FALSE}#male arearemove(list=ls())library(data.table)library(dplyr)library(formattable)library(tidyr)#data from Tuo Area_nonadjust_male <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Area_nonadjust_male.csv")FID_male <- read.csv("~/Desktop/2021lab/CH/consistency_eval/Dropbox_GWASstudy/residuals/FID_male.csv")Tuo_dat <- cbind(FID_male[,2],Area_nonadjust_male[,3:14]) #Tuo male thickness colnames(Tuo_dat)[1] <- "FID"###hao dat UKB_gclus_xgb_result_male <- read.delim("~/Desktop/2021lab/CH/consistency_eval/UKB/UKB_gclus_xgb_result_male.txt")Hao_dat <- UKB_gclus_xgb_result_male[,c(1,163:174)] #163:174 area AGER, 175:186 thickness AGERcolnames(Hao_dat)###reorder rows/make FID match Tuo_dat$order <- match(Tuo_dat$FID,Hao_dat$FID)Tuo_dat<- Tuo_dat[!(is.na(Tuo_dat$order)),] #since hao's data is a subset of tuo's we need to remove the data that Hao doesn't have##double check ## identical(Tuo_dat$FID,Hao_dat$FID) == TRUE ###reorder column orders list <- gsub("ChenA_bh_","",colnames(Tuo_dat))num <- c(1:13)for (i in 1:13) { num[i] <- grep(list[i],colnames(Hao_dat)) #warning will appear about length,it's ok #print(num) } Hao_dat <- Hao_dat[,num]
#correlation calculation results = list()counter = 0df <- as.data.frame(list)
k = 1for (i in 1:13) { a <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") # print(paste(colnames(Hao_dat)[i], " est:", a$estimate, " p=value:", a$p.value)) df$cor[i] <- round(a$estimate,5) df$pval[i] <- a$p.value df$t[i] <- a$statistic df$df[i] <- a$parameter df$CI[i] <- paste(round(a$conf.int[1],5),round(a$conf.int[2],5)) results[[i]] <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") print(results[[i]])}
df$p <- "p-value < 2.2e-16" #manually adding this note with quick visual inspection from print(results[[i]]) df <- df[c(2:13),]
``````{r male_area_plot, echo=FALSE}formattable(df, list( cor = color_tile("lightpink", "transparent") ))``````{r correlation summary2}summary(df$cor) #male area correlation summary```## THICKNESS ## Female Thickness Table
Tuo data = residuals/Thickness_nonadjust_female.csvHao data = UKB_gclus_xgb_result_female.txt;
ChenT_bh_region_area v.s ChenT_region_area_brainAGER_Female ```{r female thickness, include=FALSE}#female thicknessremove(list=ls())library(data.table)library(dplyr)library(formattable)library(tidyr)#data from Tuo Thickness_nonadjust_female <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Thickness_nonadjust_female.csv")FID_female <- read.csv("~/Desktop/2021lab/CH/consistency_eval/Dropbox_GWASstudy/residuals/FID_female.csv")Tuo_dat <- cbind(FID_female[,2],Thickness_nonadjust_female[,3:14]) #Tuo female thickness colnames(Tuo_dat)[1] <- "FID"###hao dat UKB_gclus_xgb_result_female <- read.delim("~/Desktop/2021lab/CH/consistency_eval/UKB/UKB_gclus_xgb_result_female.txt")Hao_dat <- UKB_gclus_xgb_result_female[,c(1,175:186)] #163:174 area AGER, 175:186 thickness AGERcolnames(Hao_dat)###reorder rows/make FID match Tuo_dat$order <- match(Tuo_dat$FID,Hao_dat$FID)Tuo_dat<- Tuo_dat[!(is.na(Tuo_dat$order)),] #since hao's data is a subset of tuo's we need to remove the data that Hao doesn't have##double check ## identical(Tuo_dat$FID,Hao_dat$FID) == TRUE ###reorder column orders list <- gsub("ChenT_bh_","",colnames(Tuo_dat))num <- c()for (i in 1:13) { num[i] <- grep(list[i],colnames(Hao_dat)) #warning will appear about length,it's ok #print(num)} Hao_dat <- Hao_dat[,num]
#correlation calculation cor.test(Tuo_dat[,2],Hao_dat[,2],method = "pearson")
results = list()counter = 0df <- as.data.frame(list)
k = 1for (i in 1:13) { a <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") # print(paste(colnames(Hao_dat)[i], " est:", a$estimate, " p=value:", a$p.value)) df$cor[i] <- round(a$estimate,5) df$pval[i] <- a$p.value df$t[i] <- a$statistic df$df[i] <- a$parameter df$CI[i] <- paste(round(a$conf.int[1],5),round(a$conf.int[2],5)) results[[i]] <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") print(results[[i]]) }
df$p_details <- "p-value < 2.2e-16" #manually adding this note with quick visual inspection from print(results[[i]]) df <- df[c(2:13),]

``````{r female_thickness, echo=FALSE}formattable(df, list( cor = color_tile("lightpink", "transparent")))``````{r correlation summary3}summary(df$cor) #female thickness correlation summary```## Male Thickness Table
Tuo data = residuals/Thickness_nonadjust_male.csvHao data = UKB_gclus_xgb_result_male.txt;
ChenT_bh_region_area v.s ChenT_region_area_brainAGER_Female ```{r male_thickness, include=FALSE}#male thickness#male thicknessremove(list=ls())library(data.table)library(dplyr)library(formattable)library(tidyr)#data from Tuo Thickness_nonadjust_male <- read.csv("~/Desktop/2021lab/CH/consistency_eval/residuals/Thickness_nonadjust_male.csv")FID_male <- read.csv("~/Desktop/2021lab/CH/consistency_eval/Dropbox_GWASstudy/residuals/FID_male.csv")Tuo_dat <- cbind(FID_male[,2],Thickness_nonadjust_male[,3:14]) #Tuo male thickness colnames(Tuo_dat)[1] <- "FID"###hao dat UKB_gclus_xgb_result_male <- read.delim("~/Desktop/2021lab/CH/consistency_eval/UKB/UKB_gclus_xgb_result_male.txt")Hao_dat <- UKB_gclus_xgb_result_male[,c(1,175:186)] #163:174 area AGER, 175:186 thickness AGERcolnames(Hao_dat)###reorder rows/make FID match Tuo_dat$order <- match(Tuo_dat$FID,Hao_dat$FID)Tuo_dat<- Tuo_dat[!(is.na(Tuo_dat$order)),] #since hao's data is a subset of tuo's we need to remove the data that Hao doesn't have##double check ## identical(Tuo_dat$FID,Hao_dat$FID) == TRUE ###reorder column orders list <- gsub("ChenT_bh_","",colnames(Tuo_dat))num <- c()for (i in 1:13) { num[i] <- grep(list[i],colnames(Hao_dat)) #warning will appear about length,it's ok} Hao_dat <- Hao_dat[,num]
#correlation calculation results = list()counter = 0df <- as.data.frame(list)
k = 1for (i in 1:13) { a <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") df$cor[i] <- round(a$estimate,5) df$pval[i] <- a$p.value df$t[i] <- a$statistic df$df[i] <- a$parameter df$CI[i] <- paste(round(a$conf.int[1],5),round(a$conf.int[2],5)) results[[i]] <- cor.test(Tuo_dat[,i],Hao_dat[,i],method = "pearson") print(results[[i]])}
df$p_details <- "p-value < 2.2e-16" #manually adding this note with quick visual inspection from print(results[[i]]) df <- df[c(2:13),]``````{r male_thickness_table, echo=FALSE}formattable(df, list( cor = color_tile("lightpink", "transparent")))``````{r correlation summary4}summary(df$cor) #male thickness correlation summary```

Page updated

Report abuse