提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  其他

当前栏目

【画图】如何批量展现基因表达相关性？

批量 如何 展现 基因 画图 表达 相关性
2023-06-13 09:15:24 时间

为什么要画这个图？

现在已经有明确的实验证明，跟SARS病毒一样，新冠状病毒2019-nCoV与宿主细胞的ACE2受体结合[1]。上次教程已经给大家演示了，GTEx数据库有人各组织中基因表达谱数据，下载整理这个数据可以绘制出ACE2受体在人体组织中的表达量情况以及可能的功能有哪些。

【画图】冠状病毒结合的宿主细胞受体ACE2在人组织中的表达情况

【画图】ACE2在TCGA肺癌数据的表达情况（请不要过度解读这个图的结果！）

我们平时在做分析的时候，找到目标基因预测功能以后要看看这个基因与其他基因表达相关情况，这就需要批量分析ACE2与其他基因表达量的相关性。

啰嗦几句+致谢

今天要画的图中结果可以看到一些lncRNA的表达与ACE2在肺组织中的正相关的那么课题就来了，至于没写国自然的小伙伴，这个干货可要抓住了！

涉及到批量最好用并行，下面代码涉及到并行的部分请诸位一定要改代码，电脑过热发烧什么的，可别怪站长没提醒哦！

郑重感谢小丫画图！大家可以在微店里搜索各种生信美图的代码，下面代码绝大部分是修改自Ya37，清洗数据部分做了优化以完成一对多的画图。

画图

1. 获得人肺组织表达谱数据批量并行计算

lungTMP<-read.csv("Lung.csv",header = T,check.names = F,row.names = 1)
pathway.score <- function(data,gene){
  y <- as.numeric(data[gene,])
  rownames <- rownames(data)
  dataframes<-do.call(rbind,lapply(cl=cl,rownames, function(x){
    dd  <- cor.test(as.numeric(data[x,]),y,type="spearman")
    data.frame(gene=gene,Genesets=x,cor=dd$estimate,p.value=dd$p.value )
  }))
  dataframes
}
ACE2cor<-pathway.score(data=lungTMP,gene="ENSG00000130234.10")

2. 整理画图需要的数据 

ACE2corG<-na.omit(AnnoENSG(data=ACE2cor,col="Genesets"))
ACE2corG_sig<-subset(ACE2corG,ACE2corG$p.value<0.05 & abs(ACE2corG$cor)>0.3)
ACE2corG_sig$gene<-"ACE2"
ACE2corG_sig_lnc<-Filterlnc(ACE2corG_sig)
ACE2corG_sig_lncSort<-ACE2corG_sig_lnc[order(ACE2corG_sig_lnc$cor,decreasing = T),][1:5,]
ACE2corG_sig_Sort<-ACE2corG_sig[order(ACE2corG_sig$cor),][1:5,]
ACE2corG_circle<-rbind(ACE2corG_sig_lncSort,ACE2corG_sig_Sort)

3. 整理画图数据

Links<-data.frame(Gene_1=ACE2corG_circle$gene,
                  Gene_2=ACE2corG_circle$SYMBOL,
                  Correlation=ACE2corG_circle$cor)
Corr <- data.frame(rbind(data.frame(Gene=Links[,1], Correlation=Links[,3]), 
                         data.frame(Gene=Links[,2], Correlation=Links[,3])),stringsAsFactors = F)
Corr$Index <- seq(1,nrow(Corr),1)
Corr <- Corr[order(Corr[,1]),]
corrsp<-split(Corr,Corr$Gene)
corrspe<-lapply(corrsp, function(x){x$Gene_Start<-0
                                    if (nrow(x)==1){x$Gene_End<-1}else{
                                      x$Gene_End<-sum(abs(x$Correlation))}
                                    x})
GeneID<-do.call(rbind,corrspe)
GeneID<-GeneID[!duplicated(GeneID$Gene),]
library(ggsci)
mycol <- pal_d3("category20c")(20)
n<-nrow(GeneID)
GeneID$Color<-mycol[1:n]
Corr[,2] <- abs(Corr[,2]) 
corrsl<-split(Corr,Corr$Gene)
aaaaa<-c()
corrspl<-lapply(corrsl,function(x){nn<-nrow(x)
                                    for (i in 1:nn){
                                      aaaaa[1]<-0
                                      aaaaa[i+1]<-x$Correlation[i]+aaaaa[i]}
                                      bbbbb<-data.frame(V4=aaaaa[1:nn],V5=aaaaa[2:(nn+1)])
                                      bbbbbb<-cbind(x,bbbbb)
                                      bbbbbb
                                      })
Corr<-do.call(rbind,corrspl)
Corr <- Corr[order(Corr$Index),]
Links$start_1 <- Corr$V4[1:(nrow(Corr)/2)]
Links$end_1 <- Corr$V5[1:(nrow(Corr)/2)]
Links$start_2 <- Corr$V4[(nrow(Corr)/2 + 1):nrow(Corr)]
Links$end_2 <- Corr$V5[(nrow(Corr)/2 + 1):nrow(Corr)]
color <- data.frame(colorRampPalette(c("#67BE54", "#FFFFFF", "#F82C2B"))(201))
for (i in 1:nrow(Links)){
  Links[i,8] <- substring(color[Links[i,3] * 100 + 101, 1], 1, 7)
}
names(Links)[8] <- "color"

4.画图

library(circlize)
pdf("correlation.pdf",width = 6,height = 6)
par(mar=rep(0,4))
circos.clear()
circos.par(start.degree = 90, 
           gap.degree = 5, 
           track.margin = c(0,0.23),
           cell.padding = c(0,0,0,0)
)
circos.initialize(factors = GeneID$Gene,
                  xlim = cbind(GeneID$Gene_Start, GeneID$Gene_End))
circos.trackPlotRegion(ylim = c(0, 1), factors = GeneID$Gene, 
                       track.height = 0.05,
                       panel.fun = function(x, y) {
                         name = get.cell.meta.data("sector.index") 
                         i = get.cell.meta.data("sector.numeric.index") 
                         xlim = get.cell.meta.data("xlim")
                         ylim = get.cell.meta.data("ylim")
                         circos.text(x = mean(xlim), y = 1,
                                     labels = name,
                                     cex = 0.7, 
                                     niceFacing = TRUE, 
                                     facing = "bending",
                                     adj = c(0.5, -2.8), 
                                     font = 2  
                         )
                         circos.rect(xleft = xlim[1], 
                                     ybottom = ylim[1],
                                     xright = xlim[2], 
                                     ytop = ylim[2],
                                     col = GeneID$Color[i],
                                     border = GeneID$Color[i])

                         circos.axis(labels.cex = 0.7, 
                                     direction = "outside"
                         )})

for(i in 1:nrow(Links)){
  circos.link(sector.index1 = Links$Gene_1[i], 
              point1 = c(Links[i, 4], Links[i, 5]),
              sector.index2 = Links$Gene_2[i], 
              point2 = c(Links[i, 6], Links[i, 7]),
              col = paste(Links$color[i], "C9", sep = ""), 
              border = FALSE, 
              rou = 0.7
              )}
i <- seq(0,0.995,0.005)
rect(-1+i/2, 
     -1, 
     -0.9975+i/2, 
     -0.96, 
     col = paste(as.character(color[,1]), "FF", sep = ""),
     border = paste(as.character(color[,1]), "FF", sep = ""))
text(-0.97, -1.03, "-1")
text(-0.51, -1.03, "1")

备注：上面annoE是站长自己写的注释基因的包。

画图素材：

1、在GTEx上下载其中人肺组织表达谱数据

2、需要annoE包和提取lncRNA的FilterLnc代码

3、想快速的运行的需要并行处理的代码

上面的素材，大家可以自行准备，当然如果嫌麻烦

相关文章