2022.11.8初版

误入BioInfor的大黄鸭 --一个喜欢把教程写着写着写成科普的本科临床医学生

组间相关性主要用于对比两组指标之间的相关性，可探究一个基因集与另一个基因集各基因之间表达的相关性，也可以是基因集与其他的参数。此教程的代码可在图上显示出相关系数及p值的等级。

目录

• 1.输入文件的准备
• 2.包的准备
• 3.第一种方法
• 4.第二种方法

输入文件的准备

输入文件需要两个，因为是探究两组指标的相关性，所以需要在相同的样本里对比。

我们需要准备这样的两个数据

接下来我以m6A的Writers与免疫检查点相关基因进行组间相关性检验

包的准备

install.packages("psych")
install.packages("pheatmap")
install.packages("reshape")
install.packages("tidyverse")
library("psych")
library("pheatmap")
library("reshape")
library("tidyverse")

加载数据

table1 <- read.table("immcheck.txt",header = T,check.names = F,sep = "\t",row.names = 1)
table2 <- read.table("Writers.txt",header = T,check.names = F,sep = "\t",row.names = 1)
put="Writers-immcheck" #文件名

第一种方法

data.corr <- corr.test(table1, table2, method="pearson", adjust="fdr")
data.r <- data.corr$r
data.p <- data.corr$p
data.r.fmt <- matrix(sprintf("%.2f", data.r), nrow=nrow(data.p))
getSig <- function(dc) {
  sc <- ''
  if (dc < 0.001) sc <- '***'
  else if (dc < 0.01) sc <- '**'
  else if (dc < 0.05) sc <- '*'
  sc
  }
sig.mat <- matrix(sapply(data.p, getSig), nrow=nrow(data.p))

a=paste0(data.r.fmt,sig.mat)
sig.mat<-matrix(a,nrow = ncol(table1),ncol = ncol(table2),byrow = TRUE)

pdf(paste0(put,"cor.pdf"),height = 8,width = 8)
pheatmap(data.r, clustering_method="average",main="Writers-immcheck", cluster_tabel1=F,treeheight_col = 10, cluster_tabel2=TRUE,display_numbers=sig.mat) #这块地方自己命名
dev.off()

效果图

第二种方法

pp <- corr.test(table1,table2,method="pearson",adjust = "fdr")

cor <- pp$r
pvalue <- pp$p

data.r.fmt <- matrix(sprintf("%.2f", cor), nrow=nrow(pvalue))

getSig <- function(pval) {
  stars = ""
  if(pval <= 0.001)
    stars = "***"
  if(pval > 0.001 & pval <= 0.01)
    stars = "**"
  if(pval > 0.01 & pval <= 0.05)
    stars = "*"
  if(pval > 0.05 & pval <= 0.1)
    stars = ""
  stars
}

sig.mat <- matrix(sapply(pvalue, getSig), nrow=nrow(pvalue))
a=paste0(data.r.fmt,sig.mat)

pdf(paste0(put,"cor.pdf"),height = 8,width = 15)

heatmap <- melt(cor) %>% rename(replace=c("X1"="sample","X2"="gene",
                                          "value"="cor")) %>%
  mutate(pvalue=melt(pvalue)[,3]) %>%
  mutate(signif = a)

  ggplot(heatmap,aes(sample,gene,col=cor)) +
  geom_point(aes(size = abs(cor)),shape=16) +
  geom_text(aes(label=signif),size=4,color="black",
              hjust=0.5,vjust=0.7)+
  scale_color_gradientn(colors=LaCroixColoR::lacroix_palette("PeachPear", type="continuous"),
                        breaks=c(3.9,6,8,9.8))+ 
  scale_size_area(breaks=c(2808, 10000, 19688),label=scales::comma) +
  coord_fixed() +
  theme_minimal() +
  theme(plot.margin=margin(0.1,0.1,0.3,0.3,unit="cm"),
        panel.grid=element_line(linetype = "dotted", color="grey60",size=.2),
        legend.text=element_text(size=15,color="black"),
        legend.direction = "vertical",
        legend.box = "horizontal",    # 将图例并行放置
        axis.text=element_text(color="black",size=3)) +
  scale_size(range=c(1,12),guide=NULL)+
  guides(color=guide_colorbar(barheight = unit(14, "lines"),
                              barwidth = unit(.8, "lines"),order=1),
         size=guide_legend(override.aes = list(shape=21))) 
  #+labs(x=NULL,y=NULL,color=NULL,size=NULL)

dev.off()

效果图

本教程就先讲到这啦，后续随后更新，欢迎大家关注支持～大家关注一下我v：误入BioInfor的大黄鸭，回复“TCGACox”获取完整版的代码