细胞免疫浸润网站使用

1.cibersort

数据输入格式：输入数据的要求： 1.不可以有负值和缺失值 2.不要取log 3.如果是芯片数据，昂飞芯片使用RMA标准化，Illumina 的Beadchip 和Agilent的单色芯片，用limma处理。 4.如果是RNA-Seq表达量，使用FPKM和TPM都很合适。芯片的要求可能把你唬住了，GEO常规的表达矩阵都是这样得到的，直接下载使用即可。注意有的表达矩阵下载下来就已经取过log，需要逆转回去。有的经过了标准化或者有负值，需要处理原始数据，前面写过介绍文：

https://www.jianshu.com/p/d7035ba8347b

https://www.jianshu.com/p/e3d734b2c404

做成cibersort要求的输入文件

这个算法并没有被写成R包，而是只有一个放着函数的脚本–CIBERSORT.R，把它下载下来放在工作目录即可。

需要两个输入文件:

一个是表达矩阵文件

一个是官网提供的LM22.txt，记录了22种免疫细胞的基因表达特征数据。三个输入文件的要求及形式，参考以下网址：

https://www.jianshu.com/p/0baac4c52ac8

由于CIBERSORT.R读取文件的代码比较粗暴，为了适应它，导出文件之前需要把行名变成一列。不然后面就会有报错。

exp2=as.data.frame(tpms)

exp2=rownames_to_column(exp2)

write.table(exp2,file="exp.txt",row.names=F,quote=F,sep="\t")

source("CIBERSORT.R")

if(F){

  TME.results = CIBERSORT("LM22.txt",

                          "exp.txt" ,

                          perm = 1000,

                          QN = T)

  save(TME.results,file = "ciber_CHOL.Rdata")

}

load("ciber_CHOL.Rdata")

TME.results[1:4,1:4]

## B cells naive B cells memory Plasma cells

## TCGA-W5-AA36-01A-11R-A41I-07 0.00000000 0.002351185 0.02550133

## TCGA-W5-AA2H-01A-31R-A41I-07 0.04512086 0.354414124 0.01961627

## TCGA-ZU-A8S4-11A-11R-A41I-07 0.00203370 0.000000000 0.04582565

## TCGA-WD-A7RX-01A-12R-A41I-07 0.15785229 0.000000000 0.01847074

## T cells CD8

## TCGA-W5-AA36-01A-11R-A41I-07 0.07766099

## TCGA-W5-AA2H-01A-31R-A41I-07 0.14262301

## TCGA-ZU-A8S4-11A-11R-A41I-07 0.09962641

## TCGA-WD-A7RX-01A-12R-A41I-07 0.13769951

re <- TME.results[,-(23:25)]

运行有些慢。计算出来的结果包含了22种免疫细胞的丰度，还有三列其他统计量，不管它们。

3.5. 经典的免疫细胞丰度热图

那些在一半以上样本里丰度为0的免疫细胞，就不展示在热图里了。我看了一下这个热图，从聚类的情况来看，normal和tumor没有很好的分开。

library(pheatmap)k<-apply(re,2,function(x){sum(x==0)<nrow(TME.results)/2})table(k)

## k## FALSE  TRUE ##    8    14

re2<-as.data.frame(t(re[,k]))an=data.frame(group=Group,row.names=colnames(exp))pheatmap(re2,scale="row",show_colnames=F,annotation_col=an,color=colorRampPalette(c("navy","white","firebrick3"))(50))

3.6. 直方图

可以展示出每个样本的免疫细胞比例

library(RColorBrewer)

mypalette<-colorRampPalette(brewer.pal(8,"Set1"))

dat<-re%>%as.data.frame()%>%

rownames_to_column("Sample")%>%

gather(key=Cell_type,value=Proportion,-Sample)

ggplot(dat,aes(Sample,Proportion,fill=Cell_type))+

geom_bar(stat="identity")+labs(fill="Cell Type",x="",y="Estiamted Proportion")+

theme_bw()+

theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),legend.position="bottom")+

scale_y_continuous(expand=c(0.01,0))+

scale_fill_manual(values=mypalette(22))

3.7 箱线图

展示免疫细胞之间的比较。

ggplot(dat,aes(Cell_type,Proportion,fill=Cell_type))+

geom_boxplot(outlier.shape=21,color="black")+

theme_bw()+labs(x="Cell Type",y="Estimated Proportion")+

theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),legend.position="bottom")+

scale_fill_manual(values=mypalette(22))

乱了点？那就让箱线图拥有顺序吧。

a = dat %>%  group_by(Cell_type) %>%  summarise(m = median(Proportion)) %>%  arrange(desc(m)) %>%  pull(Cell_type)dat$Cell_type = factor(dat$Cell_type,levels = a)

ggplot(dat,aes(Cell_type,Proportion,fill = Cell_type)) +  geom_boxplot(outlier.shape = 21,color = "black") +  

 theme_bw() +  

 labs(x = "Cell Type", y = "Estimated Proportion") +   

 theme(axis.text.x = element_blank(),        

axis.ticks.x = element_blank(),        

legend.position = "bottom") + 

 scale_fill_manual(values = mypalette(22))

既然我们已经把正常样本也算了，那就做个比较：

dat$Group = ifelse(as.numeric(str_sub(dat$Sample,14,15))<10,"tumor","normal")

library(ggpubr)

ggplot(dat,aes(Cell_type,Proportion,fill = Group)) +  geom_boxplot(outlier.shape = 21,color = "black") +  

 theme_bw() +  labs(x = "Cell Type", y = "Estimated Proportion") +  theme(legend.position = "top") +  

 theme(axis.text.x = element_text(angle=80,vjust = 0.5))+  scale_fill_manual(values = mypalette(22)[c(6,1)])+ stat_compare_means(aes(group = Group,label = ..p.signif..),method = "kruskal.test")

2.TImer2.0

肿瘤微环境中免疫细胞的组成和丰度是影响肿瘤进展和免疫治疗效果的重要因素。由于直接测量方法的局限性，计算算法通常用于从肿瘤转录组图谱推断免疫细胞组成。这些估计的肿瘤免疫浸润人群与肿瘤的基因组和转录组学变化相关，提供了对肿瘤免疫相互作用的见解。今天我们就一起来看这个肿瘤免疫分析的又一神器——TIMER 2.0数据库（开始的页面就已经有数据格式要求，注意看）。具体操作如下：

https://zhuanlan.zhihu.com/p/367363990

http://www.sci666.net/59270.html

https://www.jingege.wang/2021/09/27/%E7%A0%94%E7%A9%B6%E5%85%8D%E7%96%AB%E6%B5%B8%E6%B6%A6%E7%9A%84%E7%A5%9E%E5%99%A8timer%E6%9B%B4%E6%96%B0-timer2-0/

3.xcell

https://www.jianshu.com/p/a20a534233aa

4.ESTIMATE

肿瘤微环境中，免疫细胞和基质细胞是两种主要肺肿瘤细胞类型；而ESTIMATE，利用表达谱数据来预测基质细胞和免疫细胞评分，进而预测这两种细胞的含量；最后可以计算每个肿瘤样本里面的肿瘤纯度了，也就是说，基质细胞和免疫细胞含量多，那肿瘤纯度低，反之肿瘤纯度就高了。

可能对RNA-Seq数据的预测不是很准，芯片数据还可以。

代码参考：https://www.jianshu.com/p/b64e9f338f8b

另一种分析方法网页+代码参考：https://zhuanlan.zhihu.com/p/404707114

https://zhuanlan.zhihu.com/p/136747705

整理的代码及注意事项

https://www.jianshu.com/p/ec5307256ca5

library(utils)rforge<-"http://r-forge.r-project.org"

install.packages("estimate",repos=rforge,dependencies=TRUE)

library(estimate)

初步运行自带的测试数据

library(estimate)

OvarianCancerExpr<system.file("extdata","sample_input.txt",package="estimate")

filterCommonGenes(input.f=OvarianCancerExpr,output.f="OV_10412genes.gct",id="GeneSymbol")

estimateScore(input.ds="OV_10412genes.gct",output.ds="OV_estimate_score.gct",platform="affymetrix")

plotPurity(scores="OV_estimate_score.gct",samples="s519",platform="affymetrix")

从测试数据可以看出，设置好工作路径后，将.txt文件直接放在当前工作目录下，就可以不跑OvarianCancerExpr，直接跑filterCommonGenes，因为可以设置input.f=xxx.txt。

数据格式如图1：

图. 1

完整代码

setwd("E:/1.RStudio/geo_ovarian_GPL96_570/geo_array_ovarian/OVCP3")

library(utils)

rforge <- "http://r-forge.r-project.org"

install.packages("estimate", repos=rforge, dependencies=TRUE)

library(estimate)

###提前将Exp0放在当前工作目录下

filterCommonGenes(input.f="Exp0.txt",

                  output.f="OV_6934genes.gct",

                  id="GeneSymbol")

estimateScore(input.ds = "OV_6934genes.gct",

              output.ds="OV_estimate_score.gct",

              platform="affymetrix")

plotPurity(scores="OV_estimate_score.gct",

          platform="affymetrix")

###########结果读取和导出方式

scores=read.table("OV_estimate_score.gct",skip = 2,header = T)

rownames(scores)=scores[,1]

scores=t(scores[,3:ncol(scores)])

View(scores)

save(scores,file = 'BRCA_estimate_score.Rdata')