一、生信分析的目的

前面010-1.目标基因在高低表达组免疫差异分析，讲了免疫细胞在我们要研究的目标基因高低表达组之间的关系。但是对于单基因研究来说这个分析显然是不足的，因为它仅仅展示了存在差异的免疫细胞。对基因与免疫细胞的关系仍需进一步描述，就是在010.1中我们只是说明了这个免疫细胞和目标基因是在高低组中是存在差异，但是相关性无法用具体的值来描述，这一节的代码就是用具体的数值来描述二者之间的相关性。

二、具体的代码

remove(list = ls()) ##清空当前环境 
setwd("D:/生信代码复现/12.基因表达与免疫浸润关系") ##💚💚💚💚💚设置路径


#引用包
library(limma)
library(reshape2)
library(ggpubr)
library(vioplot)
library(ggExtra)

expFile="geneExp.txt"              #💚💚💚表达数据文件，这个是我们之前筛选的目的基因的表达量文件
immFile="CIBERSORT-Results.txt"    #💚💚💚免疫细胞浸润的结果文件，这个相当于是免疫细胞浸润的表达量文件
pFilter=0.05            #💚💚💚💚免疫细胞浸润结果的过滤条件，一般默认0.05


#读取表达数据文件
rt=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)
gene=colnames(rt)[1]

#删掉正常样品
tumorData=rt[rt$Type=="Tumor",1,drop=F]
tumorData=as.matrix(tumorData)
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData))
data=avereps(tumorData)

#根据目标基因表达量对样品进行分组
data=as.data.frame(data)
data$gene=ifelse(data[,gene]>median(data[,gene]), "High", "Low")

#读取免疫细胞结果文件，并对数据进行整理
immune=read.table(immFile, header=T, sep="\t", check.names=F, row.names=1)
immune=immune[immune[,"P-value"]<pFilter,]
immune=as.matrix(immune[,1:(ncol(immune)-3)])

#删除正常样品
group=sapply(strsplit(row.names(immune),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
immune=immune[group==0,]
row.names(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", row.names(immune))
immune=avereps(immune)

#数据合并
sameSample=intersect(row.names(immune), row.names(data))
rt=cbind(immune[sameSample,,drop=F], data[sameSample,,drop=F])


#把数据转换成ggplot2输入文件
data=rt[,-(ncol(rt)-1)]
data=melt(data,id.vars=c("gene"))
colnames(data)=c("gene", "Immune", "Expression")
##########绘制相关性散点图##########
outTab=data.frame()
for(i in colnames(rt)[1:(ncol(rt)-2)]){
  x=as.numeric(rt[,gene])
  y=as.numeric(rt[,i])
  if(sd(y)==0){y[1]=0.00001}
  cor=cor.test(x, y, method="spearman")
  outVector=cbind(Cell=i, cor=cor$estimate, pvalue=cor$p.value)
  outTab=rbind(outTab,outVector)
  if(cor$p.value<0.05){
    outFile=paste0("cor.", i, ".pdf")
    df1=as.data.frame(cbind(x,y))
    p1=ggplot(df1, aes(x, y)) + 
      xlab(paste0(gene, " expression")) + ylab(i)+
      geom_point() + geom_smooth(method="lm",formula = y ~ x) + theme_bw()+
      stat_cor(method = 'spearman', aes(x =x, y =y))
    p2=ggMarginal(p1, type="density", xparams=list(fill = "orange"), yparams=list(fill = "blue"))
    #相关性图形
    pdf(file=outFile, width=5.2, height=5)
    print(p2)
    dev.off()
  }
}
#输出相关性的结果文件
write.table(outTab,file="cor.result.txt",sep="\t",row.names=F,quote=F)

三、最终的效果

①R语言中的效果：

②文件效果：

四、遇到的问题

出现错误警告：

> remove(list = ls()) ##清空当前环境 
> setwd("D:/生信代码复现/12.基因表达与免疫浸润关系") ##💚💚💚💚💚设置路径
> #引用包
> library(limma)
> library(reshape2)
Warning message:
程辑包‘reshape2’是用R版本4.3.3 来建造的 
> library(ggpubr)
载入需要的程辑包：ggplot2
Warning message:
程辑包‘ggplot2’是用R版本4.3.3 来建造的 
> library(vioplot)
载入需要的程辑包：sm
Package 'sm', version 2.2-6.0: type help(sm) for summary information
载入需要的程辑包：zoo

载入程辑包：‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Warning messages:
1: 程辑包‘vioplot’是用R版本4.3.3 来建造的 
2: 程辑包‘sm’是用R版本4.3.3 来建造的 
3: 程辑包‘zoo’是用R版本4.3.3 来建造的 
> library(ggExtra)
Warning message:
程辑包‘ggExtra’是用R版本4.3.3 来建造的 
> expFile="geneExp.txt"              #💚💚💚表达数据文件，这个是我们之前筛选的目的基因的表达量文件
> immFile="CIBERSORT-Results.txt"    #💚💚💚免疫细胞浸润的结果文件，这个相当于是免疫细胞浸润的表达量文件
> pFilter=0.05            #💚💚💚💚免疫细胞浸润结果的过滤条件，一般默认0.05
> #读取表达数据文件
> rt=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)
> gene=colnames(rt)[1]
> #删掉正常样品
> tumorData=rt[rt$Type=="Tumor",1,drop=F]
> tumorData=as.matrix(tumorData)
> rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData))
> data=avereps(tumorData)
> #根据目标基因表达量对样品进行分组
> data=as.data.frame(data)
> data$gene=ifelse(data[,gene]>median(data[,gene]), "High", "Low")
> #读取免疫细胞结果文件，并对数据进行整理
> immune=read.table(immFile, header=T, sep="\t", check.names=F, row.names=1)
> immune=immune[immune[,"P-value"]<pFilter,]
> immune=as.matrix(immune[,1:(ncol(immune)-3)])
> #删除正常样品
> group=sapply(strsplit(row.names(immune),"\\-"), "[", 4)
> group=sapply(strsplit(group,""), "[", 1)
> group=gsub("2", "1", group)
> immune=immune[group==0,]
> row.names(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", row.names(immune))
> immune=avereps(immune)
> #数据合并
> sameSample=intersect(row.names(immune), row.names(data))
> rt=cbind(immune[sameSample,,drop=F], data[sameSample,,drop=F])
> #把数据转换成ggplot2输入文件
> data=rt[,-(ncol(rt)-1)]
> data=melt(data,id.vars=c("gene"))
> colnames(data)=c("gene", "Immune", "Expression")
> ##########绘制相关性散点图##########
> outTab=data.frame()
> for(i in colnames(rt)[1:(ncol(rt)-2)]){
+ x=as.numeric(rt[,gene])
+ y=as.numeric(rt[,i])
+ if(sd(y)==0){y[1]=0.00001}
+ cor=cor.test(x, y, method="spearman")
+ outVector=cbind(Cell=i, cor=cor$estimate, pvalue=cor$p.value)
+ outTab=rbind(outTab,outVector)
+ if(cor$p.value<0.05){
+ outFile=paste0("cor.", i, ".pdf")
+ df1=as.data.frame(cbind(x,y))
+ p1=ggplot(df1, aes(x, y)) + 
+ xlab(paste0(gene, " expression")) + ylab(i)+
+ geom_point() + geom_smooth(method="lm",formula = y ~ x) + theme_bw()+
+ stat_cor(method = 'spearman', aes(x =x, y =y))
+ p2=ggMarginal(p1, type="density", xparams=list(fill = "orange"), yparams=list(fill = "blue"))
+ #相关性图形
+ pdf(file=outFile, width=5.2, height=5)
+ print(p2)
+ dev.off()
+ }
+ }
There were 22 warnings (use warnings() to see them)
> #输出相关性的结果文件
> write.table(outTab,file="cor.result.txt",sep="\t",row.names=F,quote=F)
> > 
> warnings()
警告信息：
1: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
2: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
3: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
4: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
5: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
6: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
7: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
8: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
9: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
10: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
11: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
12: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
13: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
14: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
15: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
16: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
17: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
18: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
19: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
20: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
21: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值
22: In cor.test.default(x, y, method = "spearman") : 无法给连结计算精確p值

五、参考

①代码具体来源——公众号：《叉叉滴同学的生信笔记》