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”获取完整版的代码
