We begin with:
multiple .annot files that have multiple columns of 1 and 0 information for 22 chromosomes
Our end goal:
23 correlation plots
produce one correlation plot for each chromosome (22 plots - there's 22 chromosomes total in this project)
Average out the 22 plots and just make one average correlation plot (1 plot)
we begin with data (.annot files) that looks similar to this:
end goal correlational plot that looks something like this:
Before opening R
I scp all the files and manually saved them on my local computer location (/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/"). See ssh page for how to scp
location command used for scp username@ip.adress
BDir (base)
:/space/chen-syn01/1/data/chen/bin/ldsc/1000G_Phase3_baselineLD_v2.2_ldscores/*.annot
DDir (d, n, p, u1, u2, us1, us2)
:/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/DMR/invar/annot/*.annot
:/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/DMR/invar/annot/union/*.annot
:/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/motor/annot/union_sum/union_sum_chr*.annot
DDir1 (cm1, cm2)
/space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/ldsc/partitioned_h2/Mar2020/motor_annot/chromatin_motor_chr*.annot
/space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/ldsc/partitioned_h2/Mar2020/motor_annot/annot/union_sum/union_sum_chr*.annot
Aftet obtaining all your files, open R
Step by step R code expatiation:
Setting up work directory:
#setwd("/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/DMR/invar/corMat/")
#DDir="/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/DMR/invar/annot/"
DDir="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/"
DDir1="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/cmakowski/"
BDir="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/"
ODir="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/"
Note: I used scp and fetched the files to my local computer and did everything locally.
Basically we have set up 4 different directory
DDir and DDir1 are data directories; BDir is where the base files are stored; ODir is where we want to keep the output files. You can add more directories as you pull data from different locations.
Loop begin- looping through all 22 chromosomes one at a time.
for (i in 1:22){
base=read.table(paste0(BDir,'baselineLD.',i,'.annot'), h=T)
n=read.table(paste0(DDir,'Neanderthal_chr',i,'.annot'),h=T)
d=read.table(paste0(DDir,'Denisovan_chr',i,'.annot'),h=T)
p=read.table(paste0(DDir,'Presentday_Human_chr',i,'.annot'),h=T)
u1=read.table(paste0(DDir,'Unclassified_chr',i,'.annot'),h=T)
u2=read.table(paste0(DDir,'union_chr',i,'.annot'),h=T)
loading in input files we want to use.
other then the base file all other files are just one collumn of 1 or 0 (see example below)
You can add more command lines as you enter data from different places.
you need to use paste0() function for files names to be successfully recognized because we are inside a loop
The base data looks something like this:
Example of what n data looks like (all the files other then 'base' all look similar to this)
Combine the data set into one data frame
data=cbind(base,d,n,p,u1,u2,cm1,us1,cm2,us2)
data1=data[,-c(1:5)]
cormat=cor(data1)
write.table(cormat,file=paste0(ODir,'corMat',i,'.txt'),quote=F)
we don't need the first 5 columns of the data hence : data[,-c(1:5)]
cor(data1) calculates the correlation of items in data 1; the data is also transformed from a dataframe to a matrix
the last line of code saves the correlation data into .txt file in your ODir. This file will be used later when we are want to generate the average plot. In the case of this project you will output 22 files because we have files from 22 chromosomes in this project.
install.packages("corrplot")
library(corrplot)
plot.new()
tiff(paste0(file="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/chr", i, "_testplot.tiff"),
width=2500,height=2000)
corrplot(cormat, method="circle", tl.col="black")
dev.off()
} #end of for loop
This will generate 22 correlational plots for one of each chromosome. The plots will saved tiff files
corrplot(data, method="", tl.col='black'
the data has to be matrix
methods you can change to many other things some commonly useful ones:
circle
upper
tl.color changes the text label color. It is changed to black here because the default color is red.
Loop ends.
Plots will look something like this
22 plots will be generated for this project. The number of variables will depend on how many variable are in your data. In this example we have 101 variables.
We now have 22 of those plots but we want a summary of the 22 plots. So below we will proceed to make one correlation plot that is the mean of those 22 plots.
Essentially what we want to do is to add up the value that is on the same spot in each plot and divide them by the number of plots you have to get an average value.
first <- read.table("/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/corMat1.txt", header=TRUE, quote="", comment.char="")
for (i in 2:22){ #22 chromosomes
data1 <- read.table(paste0("/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/corMat",i,".txt"), header=TRUE, quote="", comment.char="")
for (j in 1:145){ #145 vraiables, update to number of variables
first[j,]<- first[j,] + data1[j,]
}
}
# first[1,1] #for double checking (should be the number of chromosomes)
averagedf <- first/22 #22 for 22 chromosomes
class(averagedf)
# averagedf[1,1] # for double checking (should be 1)
avaeragematirx <- data.matrix(averagedf)
#write.table(avaeragematirx,file='corMatAverage.txt',quote=F)
first we input the very first correlation matrix as a base
then we build a for loop
we will be looping through matrix file of chromosome 2-22 (we want to add them to the matrix)
the inner for loop (j loop) is going through every column of a matrix and adding it to the same column of the previous data set.
ie making sure column 1 in dataB is being added to column 1 in dataA
each time it loops it adds on the the previous number (ie it accumulates)
After the j loop you can check if the right amount of time is added by calling for first[1,1] the first row of the first column in correlation data should be 1 for every data set. So first[1,1] should return how ever many layers you have added to the matrix, in this case first[1,1] will return 22.
then you divide the total sum data frame by 22 (or what ever number of chromosome data you have) to fine the average
you can double check if your calculations were sucessful y by calling for averagedf[1,1] #if you did things right it should return 1
Convert the avaerage data to a matrix (beccause we want to make a corrplot later and corrplot only accept things in matrix formate)
You can save the average corrrlation into a .txt if you want
Now that we have the matrix set up from the provous step we can just take that matrix and simply make a correlation plot out of that matrix. Don't over think this!
#making average plot
library(corrplot)
plot.new()
tiff(paste0(file="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/avaerage_plot(update3).tiff"),
width=2500,height=2000)
corrplot(avaeragematirx, method="circle", tl.col="black", type='upper')
dev.off()
we are making this plot using corrplot
make sure you call for a new plot if you want to save it properly
tiff() is saving it as a tiff file. You can replace this with png or other image format you prefer
corrplot(data, method=, type= ) is where the actual correlation plot is made. Make sure the data is a matrix other wise you will run into an error:
Error in matrix(if (is.null(value)) logical() else value, nrow = nr, dimnames = list(rn, :
length of 'dimnames' [2] not equal to array extent
Error in matrix(if (is.null(value)) logical() else value, nrow = nr, dimnames = list(rn, :
length of 'dimnames' [2] not equal to array extent
This error is happening because your data is not a matrix. corrplot(data, method =, type=), the data part only accepts things in matrix format
rm(list=ls())
#setwd("/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/DMR/invar/corMat/")
DDir="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/"
#original location of files: /space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/ldsc/partitioned_h2/Mar2020/motor_annot/chromatin_motor_chr*.annot#original location of files: /space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/ldsc/partitioned_h2/Mar2020/motor_annot/annot/union_sum/union_sum_chr*.annotDDir1="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/cmakowski/"
#original location of file: BDir="/space/chen-syn01/1/data/chen/bin/ldsc/1000G_Phase3_baselineLD_v2.2_ldscores/"BDir="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/data/"
#original location of file: ODir ="/space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/DMR/invar/corMat/"ODir="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/"
for (i in 1:22){
base=read.table(paste0(BDir,'baselineLD.',i,'.annot'), h=T)
n=read.table(paste0(DDir,'Neanderthal_chr',i,'.annot'),h=T)
d=read.table(paste0(DDir,'Denisovan_chr',i,'.annot'),h=T)
p=read.table(paste0(DDir,'Presentday_Human_chr',i,'.annot'),h=T)
u1=read.table(paste0(DDir,'Unclassified_chr',i,'.annot'),h=T)
u2=read.table(paste0(DDir,'union_chr',i,'.annot'),h=T)
# /space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/motor/annot/chromatin_motor_chr*.annotcm1=read.table(paste0(DDir,'chromatin_motor_chr',i,'.annot'),h=T)
# /space/chen-syn01/1/data/chen/gwas_gclust/UKB/LDSC/make_annot/motor/annot/union_sum/union_sum_chr*.annotus1=read.table(paste0(DDir,'union_sum_chr',i,'.annot'),h=T)
#/space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/ldsc/partitioned_h2/Mar2020/motor_annot/chromatin_motor_chr*.annotcm2=read.table(paste0(DDir1,'chromatin_motor_chr',i,'.annot'),h=T)
#space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/ldsc/partitioned_h2/Mar2020/motor_annot/annot/union_sum/union_sum_chr*.annotus2=read.table(paste0(DDir1,'union_sum_chr',i,'.annot'),h=T)
data=cbind(base,d,n,p,u1,u2,cm1,us1,cm2,us2)
data1=data[,-c(1:5)]
cormat=cor(data1)
write.table(cormat,file=paste0(ODir,'corMat',i,'.txt'),quote=F)
#install.packages("corrplot")
library(corrplot)
plot.new()
tiff(paste0(file="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/chr", i, "_testplot.tiff"),
width=2500,height=2000)
corrplot(cormat, method="circle", tl.col="black")
dev.off()
}
#####average data
rm(list=ls())
first <- read.table("/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/corMat1.txt", header=TRUE, quote="", comment.char="")
for (i in 2:22){ #22 chromosomes
data1 <- read.table(paste0("/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/corMat",i,".txt"), header=TRUE, quote="", comment.char="")
for (j in 1:145){ #145 vraiables, update to number of variables
first[j,]<- first[j,] + data1[j,]
}
}
averagedf <- first/22 #22 for 22 chromosomes
avaeragematirx <- data.matrix(averagedf)
write.table(avaeragematirx,file='corMatAverage.txt',quote=F)
#making average plot
library(corrplot)
plot.new()
tiff(paste0(file="/Users/niniliu/Desktop/2020lab/ch-email/correlation matrix plots/avaerage_plot(update3).tiff"),
width=2500,height=2000)
corrplot(averagedf, method="circle", tl.col="black", type='upper')
dev.off()