转录组生信参考ontology

RNA 12. SCI 文章中肿瘤免疫浸润计算方法之 CIBER

2022-03-06  本文已影响0人  桓峰基因

免疫浸润也是近几年肿瘤研究的一个重要方向。通过表达数据即可推算出这个整体样本中究竟有哪些免疫细胞。下面我们就基于数据库数据来看下整个流程分析！

前言

我们介绍了CIBERSORT，一种从其基因表达谱表征复杂组织的细胞组成的方法。当应用于计数来自新鲜、冷冻和固定组织(包括实体肿瘤)的RNA混合物中的造血亚群时，CIBERSORT在噪声、未知混合物含量和密切相关的细胞类型方面优于其他方法。CIBERSORT将使细胞生物标志物和治疗靶标的RNA混合物大规模分析成为可能。

从事免疫相关工作的研究人员，目前只需常规普通的转录组测序数据，就能拿到该样本中各类免疫细胞如DC细胞、NK细胞、CD4+ T细胞等所占比例。

例如肿瘤微环境主要由肿瘤细胞、成纤维细胞、免疫细胞、各种信号分子和细胞外基质及特殊的理化特征等共同组成，肿瘤微环境显著影响着肿瘤的诊断、生存结局和临床治疗敏感性。其中免疫浸润也是近几年肿瘤研究的一个重要方向。

所以我们要清楚一个概念那就是肿瘤组织中并不是100%的细胞是肿瘤细胞，不同肿瘤组织的微环境都有着各自的特点。

那么简单肿瘤组织中存在着那么多不同类型的细胞，但是传统的转录组混池测序方法（也叫Bulk RNA-seq）是将组织整体的RNA表达水平进行检测，我们并不能有效区分究竟哪些细胞表达了哪些基因。

实例解析

1. 数据读取

在R语言中运行Cibersort共需要三个文件，分别是：

1）官方提供的22种细胞基因集"LM22.txt"；

2）基因表达矩阵；

3）Cibersort代码。

下面我们就仔细述说一下这三个文件如何来获得，如下：

1）LM22.txt获取方法：

在Cibersort论文中（https://www.nature.com/articles/nmeth.3337#MOESM207）下载Supplementry table 1，得到如下矩阵，另存为制表符分割的txt（"LM22.txt"）

2）基因表达矩阵

表达矩阵必须如下：

  1. 第一列是基因名，并且基因名必须为 SYMBOL 号模式；

  2. 第一行是样品名，不能有重复基因名；

  3. 第一列列名不能空白。

对于矩阵中表达量输入数据的要求，如下：

1.不可以有负值和缺失值

2.不要取log

3.如果是芯片数据，昂飞芯片使用RMA标准化，Illumina 的Beadchip 和Agilent的单色芯片，用limma处理。

4.如果是RNA-Seq表达量，使用FPKM和TPM都合适。

3）Cibersort代码

在R中新建R Script，复制以下网址中代码，保存为"Cibersort.R"

https://rdrr.io/github/singha53/amritr/src/R/supportFunc_cibersort.R

数据的读取我们仍然使用的是 TCGA-COAD 的数据集，表达数据的读取以及临床信息分组的获得我们上期已经提过，我们使用的是edgeR 软件包计算出来的差异表达结果，提取上调基因 2832 的 ENSEMBL 号，

###########基因列表
DEG=read.table("DEG-resdata.xls",sep="\t",check.names=F,header = T)
geneList<-DEG[DEG$sig=="Up",]$Row.names
table(DEG$sig)
## 
## Down   Up 
## 1296 2832

直接用之前的转录组差异分析后的数据来演示，数据格式如下

mergedataTP<-DEG[,8:ncol(DEG)]
rownames(mergedataTP)=DEG$Row.names
mergedataTP[1:5,1:3]
##                 TCGA-3L-AA1B-01A-11R-A37K-07 TCGA-4N-A93T-01A-11R-A37K-07
## ENSG00000142959                           20                           15
## ENSG00000163815                          175                          108
## ENSG00000107611                           49                           13
## ENSG00000162461                           55                           89
## ENSG00000163959                          153                          259
##                 TCGA-4T-AA8H-01A-11R-A41B-07
## ENSG00000142959                           49
## ENSG00000163815                           59
## ENSG00000107611                            6
## ENSG00000162461                           48
## ENSG00000163959                          102

Count 转 TPM

首先，提取基因长度，经过一波操作之后，得到每个基因的长度，如下：

##########提取基因长度
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("F:/demo script/gencode.v22.annotation.gtf.gz",format = "gtf")
exons <- exonsBy(txdb,by="gene")
exons_length<-lapply(exons,function(x){sum(width(reduce(x)))})
exons_length<-as.data.frame(exons_length)
exons_length1<-t(exons_length)
exons_length1<-as.data.frame(exons_length1)
colnames(exons_length1)="Length"
Gene<-gsub("\\.(\\.?\\d*)","",rownames(exons_length1))
exons_length1$Gene=as.data.frame(Gene)[,1]
exons_length1[1:5,1:2]
##                    Length            Gene
## ENSG00000000003.13   4535 ENSG00000000003
## ENSG00000000005.5    1610 ENSG00000000005
## ENSG00000000419.11   1207 ENSG00000000419
## ENSG00000000457.12   6883 ENSG00000000457
## ENSG00000000460.15   5967 ENSG00000000460

计算 TPM，同样经过一波操作，我们计算得到TPM的矩阵，如下：

##############计算FPKM/RPKM/TPM
le = exons_length1[match(rownames(mergedataTP),exons_length1$Gene),"Length"]
head(le)
## [1]  2047  1105 13168  3195  5293  1209
#这个函数是现成的。
countToTpm <- function(counts, effLen)
{
  rate <- log(counts) - log(effLen)
  denom <- log(sum(exp(rate)))
  exp(rate - denom + log(1e6))
}

#install.packages("tibble")
library(tibble)
tpms <- apply(mergedataTP,2,countToTpm,le)
tpms[1:3,1:3]
##                 TCGA-3L-AA1B-01A-11R-A37K-07 TCGA-4N-A93T-01A-11R-A37K-07
## ENSG00000142959                    13.540046                     11.98080
## ENSG00000163815                   219.474345                    159.79893
## ENSG00000107611                     5.156847                      1.61412
##                 TCGA-4T-AA8H-01A-11R-A41B-07
## ENSG00000142959                    61.137762
## ENSG00000163815                   136.370688
## ENSG00000107611                     1.163758
exp2 = as.data.frame(tpms)
exp2 = rownames_to_column(exp2)
exp2[1:5,1:3]
##           rowname TCGA-3L-AA1B-01A-11R-A37K-07
## 1 ENSG00000142959                    13.540046
## 2 ENSG00000163815                   219.474345
## 3 ENSG00000107611                     5.156847
## 4 ENSG00000162461                    23.856120
## 5 ENSG00000163959                    40.058762
##   TCGA-4N-A93T-01A-11R-A37K-07
## 1                     11.98080
## 2                    159.79893
## 3                      1.61412
## 4                     45.54404
## 5                     80.00373

但是由于我们下载是Count 矩阵基因ID 是 ENSEMBL 号，所以我们需要将其转换为SYMBOL的标识，如下：

############ENSEMBL to SYMBOL
library(org.Hs.eg.db)
library(clusterProfiler)
genes<-exp2$rowname
eg <- bitr(genes, 
           fromType="ENSEMBL", 
           toType=c("ENTREZID","ENSEMBL",'SYMBOL'),
           OrgDb="org.Hs.eg.db")
head(eg)
##           ENSEMBL ENTREZID   SYMBOL
## 1 ENSG00000142959   266675    BEST4
## 2 ENSG00000163815     7123   CLEC3B
## 3 ENSG00000107611     8029     CUBN
## 4 ENSG00000162461   284723 SLC25A34
## 5 ENSG00000163959   200931   SLC51A
## 6 ENSG00000144410   130749      CPO

合并数据，并且去掉重复的基因行，最后获得第一列为基因symbol， 第一行为样本编号，中间数值为TPM，如下：

exp3<-merge(eg,exp2,by.x="ENSEMBL",by.y="rowname",all.x=TRUE)
exp4<-exp3[,3:ncol(exp3)]
exp5<-exp4[!duplicated(exp4$SYMBOL),]
exp5[1:5,1:4]
##   SYMBOL TCGA-3L-AA1B-01A-11R-A37K-07 TCGA-4N-A93T-01A-11R-A37K-07
## 1 DBNDD1                  145.7848372                     932.6654
## 2  HSPB6                  463.5507914                     125.3295
## 3   PDK4                  199.3329095                     223.2314
## 4  ABCB5                    0.7485274                       0.0000
## 5 PRSS22                  441.2341966                    1428.1244
##   TCGA-4T-AA8H-01A-11R-A41B-07
## 1                    477.12059
## 2                     35.59667
## 3                    241.33307
## 4                      0.00000
## 5                    247.16721
write.table(exp5,file = "exp.txt",row.names = F,quote = F,sep = "\t")

2. CIBERSORT 分析结果

我们必要保证 'Cibersort.R'，'LM22.txt' 和 'exp.txt' 三个文件必须在同一个目录下，并且在同一文件夹下可以得到运算结果（"CIBERSORT-Results.txt"）。Cibersort结果的默认文件名为CIBERSORT-Results.txt，在同一文件夹下进行第二次运算会覆盖第一次得到的文件，建议在每一次运算之后对文件重命名。如果表达矩阵中基因不能完全覆盖LM22.txt中的基因，Cibersort 同样可以正常运行，但不能少于 "LM22.txt" 中所需基因的一半。

一条命令计算完毕，如下：

source('Cibersort.R')
CIBERSORT('LM22.txt','exp.txt', perm = 1000, QN = T) #perm置换次数=1000，QN分位数归一化=TRUE

读取结果文件 "CIBERSORT-Results.txt"，我们看下获得的结果，如下：

results<-read.table("CIBERSORT-Results.txt",header=TRUE,row.names = 1,check.names = FALSE,sep="\t")
results[1:5,]
##                              B cells naive B cells memory Plasma cells
## TCGA-3L-AA1B-01A-11R-A37K-07             0     0.01395401   0.03348333
## TCGA-4N-A93T-01A-11R-A37K-07             0     0.03180062   0.11197878
## TCGA-4T-AA8H-01A-11R-A41B-07             0     0.02393778   0.11489327
## TCGA-5M-AAT4-01A-11R-A41B-07             0     0.01961242   0.01678773
## TCGA-5M-AAT5-01A-21R-A41B-07             0     0.03367153   0.10927705
##                              T cells CD8 T cells CD4 naive
## TCGA-3L-AA1B-01A-11R-A37K-07   0.6226358                 0
## TCGA-4N-A93T-01A-11R-A37K-07   0.4945048                 0
## TCGA-4T-AA8H-01A-11R-A41B-07   0.5156207                 0

我们看到后面几列的解释，如下：

  1. Sample：样本

  2. T Cell..：细胞浸润比例

  3. P-value : 置换检验（蒙特卡罗方法），越小越可信

  4. Correlation ：原表达矩阵乘以细胞占比后的数据矩阵与原表达矩阵的相关性

  5. RMSE：均方根误差(Root Mean Squared Error)，越小效果越好。

3. 结果展示

对于肿瘤的免疫浸润分析结果，有多种表现形式，我们这里主要说几种 SCI 文件中常见的图形，如热图，直方图，箱线图等。

2.1 经典的免疫细胞丰度热图

library(pheatmap)
re <- results[,-(23:25)]
k <- apply(re,2,function(x) {sum(x == 0) < nrow(results)/2})
re2 <- as.data.frame(t(re[,k]))
##########构造注释文件
an = data.frame(Group = group$Group,
                row.names = group$Sample)
###
#####修改色带，breaks
bk <- c(seq(-15,-5,by=1),seq(-4.9,4.9,by=0.2),seq(5,15,by=1))
pheatmap(re2,scale = "row",
         show_colnames = F,
         cluster_cols = F,
         annotation_col = an,
         drop_levels = TRUE,
         color = c(rep("blue",11),colorRampPalette(colors = c("blue","white","red"))(50),
                   rep("red",11)),
         breaks = bk,
         #legend_breaks = c(-5,-2,0,2,5)
         )

2.2 直方图

##########直方图
library(RColorBrewer)
library(tidyr)
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))

2.3 箱式图

#####箱式图
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))

2.4 排序的箱式图

###########排序
library(dplyr)
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))

2.5 添加分组显著性

############tumor -vs- normal
library(stringr)
library(ggpubr)
dat$Group = ifelse(as.numeric(str_sub(dat$Sample,14,15))<10,"TP","NT")

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")

个人感觉总结的已经很全面了，不足的地方请大家指正。

关注公众号，每日更新，扫码进群交流不停歇，马上就出视频版，关注我，您最佳的选择！


桓峰基因 生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你 49篇原创内容 --> 公众号

References:

  1. Newman AM, Alizadeh AA. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019 Jul;37(7):773-782. doi: 10.1038/s41587-019-0114-2. Epub 2019 May 6. PMID: 31061481; PMCID: PMC6610714.

  2. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol. 2018;1711:243-259. doi: 10.1007/978-1-4939-7493-1_12. PMID: 29344893; PMCID: PMC5895181.

  3. Li B, Cui Y, Nambiar DK, Sunwoo JB, Li R. The Immune Subtypes and Landscape of Squamous Cell Carcinoma. Clin Cancer Res. 2019 Jun 15;25(12):3528-3537. doi: 10.1158/1078-0432.CCR-18-4085. Epub 2019 Mar 4. PMID: 30833271; PMCID:

  4. Newman, A., Liu, C., Green, M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 12, 453--457 (2015).

上一篇下一篇

猜你喜欢

热点阅读