生信可视化(part4)--相关性图

学习笔记，仅供参考，有错必纠

文章目录

• ​​生信可视化​​

• ​​相关性散点图​​
• ​​相关性热图​​
• ​​相关性网络图​​

生信可视化

相关性散点图

相关性散点图可以展示两个gene之间的相关性.

输入数据格式：

其中，行表示gene，列表示样品.

代码：

library(ggplot2)
library(ggpubr)
library(ggExtra)

inputFile="./22.cor/input.txt"      
gene1="MTMR14"             #第一个基因名字
gene2="PRKCD"              #第二个基因名字

#读取输入文件，提取基因表达量
rt=read.table(inputFile,sep="\t",header=T,check.names=F,row.names=1)
x=as.numeric(rt[gene1,])
y=as.numeric(rt[gene2,])

#相关性分析
df1=as.data.frame(cbind(x,y))
corT=cor.test(x,y,method="spearman")
cor=corT$estimate
pValue=corT$p.value
ggplot(df1, aes(x, y)) + 
      xlab(gene1)+ylab(gene2)+
      geom_point()+ geom_smooth(method="lm",formula = y ~ x) + theme_bw()+
      stat_cor(method = 'spearman', aes(x =x, y =y))

下图中，x轴表示MTMR14基因的表达，y轴表示PRKCD基因的表达，图中蓝色的线表示回归线. 图左上方的表示相关系数，则表示相关系数的P值.

相关性热图

相关性热图可以在一张图中展示多个gene之间的相关性.

输入数据格式：

其中，行表示gene，列表示样品.

代码：

library(corrplot)
inputFile="./23.corrplot/input.txt"       

rt=read.table(inputFile,sep="\t",header=T,row.names=1)      #读取文件
rt=t(rt)      #数据转置
M=cor(rt)     #相关型矩阵

#绘制相关性图形
corrplot(M,
         method = "circle",
         order = "hclust", #聚类
         type = "upper",
         col=colorRampPalette(c("green", "white", "red"))(50)
         )

下图中，圆圈越红则两个gene之间的相关系数越偏向于1，圆圈越绿则两个gene之间的相关系数越偏向于-1.

相关性网络图

输入数据格式：

其中，行表示gene，列表示样品.

代码：

library(igraph)
library(reshape2)

inputFile="./25.corNetwork/input.txt"
cutoff=0.4  #相关性阈值

#读取输入文件
data=read.table(inputFile,header=T,sep="\t",row.names=1,check.names=F)
cordata=cor(t(data))

#保留相关性矩阵的一半
mydata = cordata
upper = upper.tri(mydata)
mydata[upper] = NA

#把相关性矩阵转换为数据框
df = data.frame(gene=rownames(mydata),mydata)
dfmeltdata = melt(df,id="gene")
dfmeltdata = dfmeltdata[!is.na(dfmeltdata$value),]
dfmeltdata = dfmeltdata[dfmeltdata$gene!=dfmeltdata$variable,]
dfmeltdata = dfmeltdata[abs(dfmeltdata$value)>cutoff,]

#定义网络图的节点和边
corweight = dfmeltdata$value
weight = corweight+abs(min(corweight))+5
d = data.frame(p1=dfmeltdata$gene,p2=dfmeltdata$variable,weight=dfmeltdata$value)
g = graph.data.frame(dfmeltdata,directed = FALSE)

#设置颜色，节点大小，字体大小
E(g)$color = ifelse(corweight>0,rgb(254/255,67/255,101/255,abs(corweight)),rgb(0/255,0/255,255/255,abs(corweight)))
V(g)$size = 8
V(g)$shape = "circle"
V(g)$lable.cex = 1.2
V(g)$color = "white"
E(g)$weight = weight

#可视化
layout(matrix(c(1,1,1,0,2,0),byrow=T,nc=3),height=c(6,1),width=c(3,4,3))
par(mar=c(1.5,2,2,2))
vertex.frame.color = NA
plot(g,layout=layout_nicely,vertex.label.cex=V(g)$lable.cex,edge.width = E(g)$weight,edge.arrow.size=0,vertex.label.color="black",vertex.frame.color=vertex.frame.color,edge.color=E(g)$color,vertex.label.cex=V(g)$lable.cex,vertex.label.font=2,vertex.size=V(g)$size,edge.curved=0.4)

#绘制图例
color_legend = c(rgb(254/255,67/255,101/255,seq(1,0,by=-0.01)),rgb(0/255,0/255,255/255,seq(0,1,by=0.01)))
par(mar=c(2,2,1,2),xpd = T,cex.axis=1.6,las=1)
barplot(rep(1,length(color_legend)),border = NA, space = 0,ylab="",xlab="",xlim=c(1,length(color_legend)),horiz=FALSE,
        axes = F, col=color_legend,main="")
axis(3,at=seq(1,length(color_legend),length=5),c(1,0.5,0,-0.5,-1),tick=FALSE)

下图中的每个节点表示一个gene，若两个gene之间有连线，则他们具有共表达的关系. 正相关用红色表示，负相关用蓝色表示. 边的颜色越深，则表示两个gene之间的相关系数越大.