肿瘤免疫细胞浸润与临床相关性分析
比如文章：
Profiles of immune infiltration in colorectal cancer and theirclinical significant: A gene expression- based study
肿瘤免疫细胞浸润是指免疫细胞从血液中移向肿瘤组织，开始发挥它的作用，可以从肿瘤组织中分离出的浸润免疫细胞。
肿瘤中免疫细胞的浸润与临床结果密切相关，肿瘤中浸润的免疫细胞最有可能作为药物靶标来提高患者的生存率。
分析用到的软件是CIBERSORT，这个软件的原理就是把我们输入的基因矩阵转换为免疫细胞的矩阵，这个矩阵的列名是细胞名称，行名是样品。当然我们需要根据p值对样品进行筛选。
CIBERSORT，是由Binder G等开发的反卷积算法，可以基于标准化的基因表达数据来计算复杂组织的细胞组成，改方法能量化特定细胞类型的丰度。通过流式细胞分选术很好地验证并成功评估了乳腺癌及肝癌组织中免疫细胞的组成。具体可浏览官网：https://cibersortx.stanford.edu/
下面我们就介绍肿瘤免疫浸润与临床相关性分析
一.表达数据准备
首先，我们需要准备一个表达矩阵，这里我们的案例用TCGA-LUAD的表达数据，我们之前已经下载TCGA数据库33个Project的RNA-Seq转录组数据上传到百度网盘，你可以去下载，也可以自己准备。差异表达分析的文章中也有介绍了从TCGA数据库下载数据的视频教程：一文就会TCGA数据库基因表达差异分析。
options(stringsAsFactors = F)
library(stringr)
library(limma)
library(survival)
library(forestplot)
library(pheatmap)
library(ggplot2)
library(limma)
library(corrplot)
library(limma)
#处理表达数据，作为免疫浸润计算的输入
exp = read.table("TCGA-COAD-Exp.txt",header = T,sep="\t",check.names = F)
#########重复的基因取平均值
delDuplicateGene = function(exp){
geneNames= exp[,"gene_id"]
exp = exp[,-1]
exp = as.matrix(exp)
rownames(exp) = as.vector(geneNames)
exp = avereps(exp)# 需要limma包
exp = as.data.frame(exp)
return(exp)
}
####################删除0方差列
del0varcol = function(exp){
# 从数据集中删除零方差列，可以使用相同的apply 表达式，设置方差不等于零。
exp =t(exp)
exp = exp[,which(apply(exp, 2, var)!=0)]
exp = t(exp)
return(exp)
}
expdata = delDuplicateGene(exp)
expdata = del0varcol(expdata)
expdata = data.frame(gene = rownames(expdata),expdata,check.names = F)
#写出作为免疫浸润的输入文件
write.table(expdata,file = "as.ImmueInputExp.txt",row.names = F,sep = "\t")上面代码中TCGA-COAD-Exp.txt文件，没有去掉重复的基因，所以我们需要处理一下，也同时删除没有表达的基因。输出的as.ImmueInputExp.txt文件就作为计算免疫细胞浸润的比例文件，第一列是基因名称，其他列是样本名称。
二.免疫浸润计算
登录网址，没有账号得先注册
https://cibersortx.stanford.edu/
可以按照下面流程上传数据，具体不解释啦。
开始计算
几分钟后计算结束
点击进去看
我们可以下载txt或者csv等其他格式的文件，用于后续的分析。
其实，我们也可以下载CIBERSORT的R代码，自己在本地运算，如果没有自己的服务器，就用人家的在先工具吧。要不然计算时间无法想象。
三.免疫浸润结果可视化
TCGAProject = "TCGA-COAD"
immuneFile = ".\\TCGA-COAD-Mixture.txt"
###创建文件夹
dir.create(TCGAProject)
setwd(TCGAProject)
dir.create("ImmuneInfiltration",showWarnings = F)
setwd("..\\")
outputPath = paste0(TCGAProject,"\\ImmuneInfiltration\\")
############一个函数，一个Barcode向量，区分正常组织和肿瘤组织
getSampleInfo <- function(Barcode){
TumorBarcode <-c()
NormalBarcode <-c()
for (id in Barcode) {
ids <- str_split(id,'-',simplify = T)[4]
ids <- as.numeric(substr(ids,1,1))
if(ids == 0){
TumorBarcode <- c(TumorBarcode,id)
}else if (ids == 1) {
NormalBarcode <-c(NormalBarcode,id)
}
}
sampleNum <- data.frame(NormalBarcode = length(NormalBarcode),TumorBarcode = length(TumorBarcode))
sampleList <- data.frame(ID = c(NormalBarcode,TumorBarcode),Type = c(rep("Normal",length(NormalBarcode)),
rep("Tumor",length(TumorBarcode))))
return(list(BarcodeSort = c(NormalBarcode,TumorBarcode),
sampleNum = sampleNum ,sampleList = sampleList))
}
1.过滤数据
就是删除p值小于0.05的样本，还有就是免疫细胞比例为0的列。
immuedata = read.table(immuneFile,header = T,row.names = 1,
sep = "\t",check.names = F)
immue = immuedata[immuedata$`P-value` < 0.05,]
immue = immue[,-c((ncol(immue)-3):ncol(immue))]#删除后3列不是免疫细胞的数据
SamplesInfo = getSampleInfo(rownames(immue))
sampleList = SamplesInfo$sampleList
normalBarcode = sampleList[sampleList$Type=="Normal","ID"] %>% as.vector()
tumorBarcode = sampleList[sampleList$Type=="Tumor","ID"] %>% as.vector()
immue = immue[c(normalBarcode,tumorBarcode),]
##其实上面的5行代码有点多余，因为我们给的表达数据文件已经排序过
immue = immue[,which(apply(immue, 2, var)!=0)]#删除0方差列的细胞2.绘制barplot
####################barplot绘制
barplotData = immue
barplotData=t(barplotData)
col=rainbow(nrow(barplotData),s=0.7,v=0.7)
pdf(paste0(outputPath,"barplot.pdf"),height=20,width=100)
par(las=1,mar=c(8,4,4,15))
a1 = barplot(barplotData,col=col,yaxt="n",ylab="Relative Percent",xaxt="n")
a2=axis(2,tick=F,labels=F)
axis(2,a2,paste0(a2*100,"%"))
axis(1,a1,labels=F)
par(srt=60,xpd=T);text(a1,-0.02,colnames(barplotData),adj=1,cex=0.6);par(srt=0)
ytick2 = cumsum(barplotData[,ncol(barplotData)])
ytick1 = c(0,ytick2[-length(ytick2)])
#text(par('usr')[2],(ytick1+ytick2)/2,rownames(barplotData),cex=0.6,adj=0)
legend(par('usr')[2]*0.98,par('usr')[4],legend=rownames(barplotData),col=col,pch=15,bty="n",cex=1.3)
dev.off()
###
图片太长，只是作为探索性的看一下免疫细胞在不同样本中的比例，可分组取均值作柱状图。
3.热图
heatmapdata = immue
heatmapdata = t(heatmapdata)
heatmapdata = heatmapdata[,c(normalBarcode,tumorBarcode)]
library(pheatmap)
Type=c(rep("Normal",length(normalBarcode)),rep("Tumor",length(tumorBarcode))) #修改对照和处理组样品数目
names(Type)=colnames(heatmapdata)
Type=as.data.frame(Type)
pdf(paste0(outputPath,"heatmap.pdf"),height=8,width=50)
pheatmap(heatmapdata,
annotation=Type,
color = colorRampPalette(c("green", "black", "red"))(50),
cluster_cols =F,
fontsize = 8,
fontsize_row=7,
fontsize_col=6)
dev.off()
4.相关性热图
###############相关性热图
corheatplotdata = immue
library(corrplot)
pdf(paste0(outputPath,"corHeatmap.pdf"),height=13,width=13) #保存图片的文件名称
corrplot(corr=cor(corheatplotdata),
method = "color",
order = "hclust",
tl.col="black",
addCoef.col = "black",
number.cex = 0.8,
col=colorRampPalette(c("blue", "white", "red"))(50),
)
dev.off()
也就是看一下免疫细胞之间的相关性
5.小提琴图
#############小提琴图
library(vioplot) #引用包
normal= SamplesInfo[["sampleNum"]][["NormalBarcode"]]#正常样品数目
tumor= SamplesInfo[["sampleNum"]][["TumorBarcode"]]#肿瘤样品数目
vioplotdata = immue
pdf(paste0(outputPath,"vioplot.pdf"),height=8,width=15)#保存图片的文件名称
par(las=1,mar=c(10,6,3,3))
x=c(1:ncol(vioplotdata))
y=c(1:ncol(vioplotdata))
plot(x,y,
xlim=c(0,(ncol(vioplotdata)-1)*3),#xlim=c(0,63)这个与细胞个数(列数)有关
#xlim=c(0,57)可以设置成列数-1，再×3
ylim=c(min(vioplotdata),max(vioplotdata)+0.02),
main="",xlab="", ylab="Fraction",
pch=21,
col="white",
xaxt="n")
#对每个免疫细胞循环，绘制vioplot，正常用green表示，肿瘤用blue表示
for(i in 1:ncol(vioplotdata)){
normalData=vioplotdata[1:normal,i]
tumorData=vioplotdata[(normal+1):(normal+tumor),i]
vioplot(normalData,at=3*(i-1),lty=1,add = T,col = 'green')
vioplot(tumorData,at=3*(i-1)+1,lty=1,add = T,col = 'blue')
wilcoxTest=wilcox.test(normalData,tumorData)
p=round(wilcoxTest$p.value,3)
mx=max(c(normalData,tumorData))
lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))
text(x=3*(i-1)+0.5,y=mx+0.02,labels=ifelse(p<0.001,paste0("p<0.001"),paste0("p=",p)),cex = 0.8)
text(seq(1,(ncol(vioplotdata)-1)*3+1,3),-0.05,xpd = NA,labels=colnames(vioplotdata),cex = 1,srt = 45,pos=2)
}
dev.off()
四.免疫浸润与临床相关性分析
得到免疫浸润结果后，就需要探讨免疫浸润与临床相关性，比如与肿瘤分期什么关系，我们演示的临床数据之前也上传到了网盘，具体看文章：TCGA数据库33个Project的RNA-Seq转录组数据为你整理打包好了。你也可以自己下载相应的TCGA-COAD临床数据，可参考文章：一文解决TCGA数据库临床数据下载与整理。这样演示的是TCGA-COAD的临床数据文件：TCGA-COAD-clinical.tsv
1.临床数据的处理
本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！