2
$\begingroup$

I am doing a gene-gene pairwise correlation, I get a plot and I do see positive and negative correlation correlation plot

My code:

library(corrplot)
library(Hmisc)
gen = read.csv('NEW_RBP/CORPLOT/TF_NORMAL_NEW_40.txt',row.names = 1)
head(gen)
dim(gen)
gendat <- t(gen)
#gendat <- gen
head(gendat)
dim(gendat)


macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
p_mat <- cor_5$P
row_names <- rownames(M)

pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")


corrplot(M, method = "color",is.corr = FALSE,
         tl.srt = 45,diag = TRUE,tl.pos='n'
         ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dev.off()

How do I extract those gene pairs? One way what I have read and done is:

library(dplyr)
flat_cor_mat <- function(cor_r, cor_p){
    #This function provides a simple formatting of a correlation matrix
    #into a table with 4 columns containing :
    # Column 1 : row names (variable 1 for the correlation test)
    # Column 2 : column names (variable 2 for the correlation test)
    # Column 3 : the correlation coefficients
    # Column 4 : the p-values of the correlations
    library(tidyr)
    library(tibble)
    cor_r <- rownames_to_column(as.data.frame(cor_r), var = "row")
    cor_r <- gather(cor_r, column, cor, -1)
    cor_p <- rownames_to_column(as.data.frame(cor_p), var = "row")
    cor_p <- gather(cor_p, column, p, -1)
    cor_p_matrix <- left_join(cor_r, cor_p, by = c("row", "column"))
    cor_p_matrix
}


#cor_3 <- rcorr(as.matrix(gendat))

my_cor_matrix <- flat_cor_mat(cor_5$r, cor_5$P)
head(my_cor_matrix)

dim(my_cor_matrix)

Delete.na <- function(DF, n=0) {
      DF[rowSums(is.na(DF)) <= n,]
    }

TF_NORMAL_CORR = Delete.na(my_cor_matrix)
dim(RBP_DISEASE_40_CORR)

TF_NORMAL_HIGH = TF_NORMAL_CORR[ TF_NORMAL_CORR$cor > 0.9 & RBP_DISEASE_40_CORR$p < 0.00000001, ]

TF_NORMAL_LOW = TF_NORMAL_CORR[TF_NORMAL_CORR$cor < -0.9 & TF_NORMAL_CORR$p < 0.00000001, ]

Now even if I do that sort of high cut-off filtering I still get considerably high number of correlation pairs like few thousands if positive correlation is considered.

So few of my issues which I would like to get input or help

  • How do I extract the cluster of regions with positive correlation or negative correlation ?
  • If I found certain gene pair/pairs how do I label them ? should I make a vector or map or any other way like if say more than 0.8 correlation I would like to put a label on those regions only
  • Considering label the gene order is important to be considered, what I have read and found is this how to find the order of the label
tempCor <- cor(gendat)
tempCor <- data.frame(tempCor)
dim(tempCor)
names(tempCor) <- c(1:741)

out <- corrplot(t(tempCor),method = "color",is.corr = FALSE,
                 tl.srt = 45,diag = TRUE,tl.pos='n'
                 ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dimnames(out)

Now with all these how can I extract those patches of high and low correlation pairs and add selected label or highlight them and add label

I m not sure if all these can be done with corrplot! like adding selected label if yes i would be really glad to use that..

My data file data used

"In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples."

Its symmetric so would it matter if i use upper left or bottom right? I guess anyone would do good.

Format - a list for each group if i understand if there is one gene with many positive correlation then that would be considered as group?is that you are saying ?

I have 16 samples in this case but here I using gene-gene. I'm bit curious with a naive question: Is it possible to do gene to sample correlation or if i know what genes are positively or negatively correlated can I trace back without have to plot another heatmap of found genes which are positively or negatively correlated to see their expression in sample ?

I have given the data file it might help you better with my context..

$\endgroup$
2
  • 1
    $\begingroup$ In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples. $\endgroup$
    – llrs
    Dec 18, 2018 at 21:52
  • $\begingroup$ I have updated my question with your queries ..hope it helps $\endgroup$
    – kcm
    Dec 19, 2018 at 4:49

1 Answer 1

4
$\begingroup$

Here is a solution using the pheatmap library to cluster and visualise the correlation matrix, then extract the groups from the cluster dendrograms:

gen <- read.csv(file = "TF_NORMAL_NEW_40.txt", header = TRUE, row.names = 1, sep = ",")
gendat <- t(gen)
library(corrplot)
library(Hmisc)
macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
str(M)
p_mat <- cor_5$P
row_names <- rownames(M)

library(pheatmap)
pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")
plot_data_pheatmap <- pheatmap(M, color = rev(macolor), clustering_method = "complete", fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

jpeg version of pdf

tree_cut <- cutree(plot_data_pheatmap$tree_row, h = 20)
tc <- data.frame(tip = names(tree_cut), clust_membership = as.character(unname(tree_cut)))
row.names(tc) <- tc$tip
tc <- tc['clust_membership']
head(tc)
str(tc)
class(tc$clust_membership)
pdf("TF_NORMAL_NEW_40_gene_clustcut.pdf",20,15,paper = "a4r")
pheatmap(M, color = rev(macolor), clustering_method = "complete", annotation_row = tc, annotation_col = tc, fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

clusters at height 20 jpeg version

$\endgroup$
4
  • $\begingroup$ okay! but how about the question i asked regarding high and low correlation? you saying i should be giving the corelation matrix as input to the pheatmap ? $\endgroup$
    – kcm
    Dec 19, 2018 at 4:40
  • $\begingroup$ @krushnachChandra Yes. Input the correlation matrix to pheatmap. I have edited my answer. $\endgroup$ Dec 19, 2018 at 5:08
  • $\begingroup$ okay i will give it a try so its the same thing as plotting in pheatmap there won't be any change i suppose.. $\endgroup$
    – kcm
    Dec 19, 2018 at 6:19
  • 1
    $\begingroup$ @krushnachChandra I have updated my answer to use your data and replicate what you have done to explain my method $\endgroup$ Dec 19, 2018 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.