免疫亚型的鉴定与可视化

测试开头

今天是生信星球陪你的第830天 免疫亚型的鉴定与可视化

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我!

1.背景知识

这是一篇2018年的22分文章，《The Immune Landscape of Cancer》将肿瘤样本划分为六种免疫亚型:wound healing, IFN-γ dominant, inflammatory, lymphocyte depleted, immunologically quiet, and TGF-β dominant,R包ImmuneSubtypeClassifier可以根据基因表达矩阵来对样本划分亚型。

多说一句，TCGA 样本的免疫亚型不需要自己计算，在文章附表2中已经有了。那就用于其他来源的表达矩阵咯。

2.示例数据和代码

出自R包的readme文件

#devtools::install_github("Gibbsdavidl/ImmuneSubtypeClassifier")

library(readr)
library(ImmuneSubtypeClassifier)

download.file(url = 'https://raw.githubusercontent.com/CRI-iAtlas/shiny-iatlas/develop/data/ebpp_test1_1to20.tsv', destfile = 'ebpp_test.tsv')
dat <- read_tsv('ebpp_test.tsv')

dat2 <- as.data.frame(dat[!duplicated(dat$GeneID),])
Xmat <- dat2[,-1]
rownames(Xmat) <- dat2[,1]

Xmat[101:104,1:4]

##            XY1     XY2      XY3      XY4
## ABCE1 1375.600 1202.25  953.933  909.373
## ABCF1 1239.670 1664.88 1192.710  635.941
## ABCF2  978.866 1590.79 2497.410 1848.180
## ABCF3 8834.770 3333.38 2428.680 5329.980

Xmat就是作者最终整理成的输入数据了。是一个未经log的表达矩阵

res0 <- callEnsemble(X = Xmat, geneids = 'symbol')
res0

##    SampleIDs BestCall            1            2            3            4
## 1        XY1        4 1.121173e-07 1.873900e-06 6.311234e-02 0.9027497470
## 2        XY2        3 4.243225e-02 2.309184e-06 5.495564e-01 0.0084770960
## 3        XY3        4 3.095071e-04 1.029907e-06 1.635270e-01 0.9063920975
## 4        XY4        4 1.154656e-04 3.823049e-07 2.888787e-02 0.8831390142
## 5        XY5        4 1.193054e-07 8.741189e-06 5.060996e-01 0.9260828793
## 6        XY6        4 2.479082e-04 5.528710e-05 1.183165e-04 0.9923758209
## 7        XY7        6 9.421768e-03 1.079501e-03 1.000559e-01 0.0084094275
## 8        XY8        3 6.143126e-04 2.594000e-06 8.872822e-01 0.1306741871
## 9        XY9        4 1.001003e-04 4.479314e-06 3.989228e-03 0.9793768525
## 10      XY10        2 4.624279e-06 9.888145e-01 7.359737e-06 0.0053910189
## 11      XY11        4 5.837736e-05 7.877929e-06 1.970483e-03 0.9895575941
## 12      XY12        3 1.944198e-06 2.005060e-06 3.616153e-01 0.5258071870
## 13      XY13        4 9.002434e-07 9.563100e-04 2.673559e-06 0.9931332767
## 14      XY14        4 1.235332e-05 3.907399e-05 2.785671e-03 0.9972651005
## 15      XY15        4 4.831095e-05 1.127145e-04 1.685999e-04 0.9935480356
## 16      XY16        3 2.731656e-03 1.722274e-05 9.742068e-01 0.0007705171
## 17      XY17        4 3.431090e-06 1.327718e-03 1.570120e-05 0.9971910715
## 18      XY18        3 3.816116e-07 5.000733e-03 8.685111e-01 0.0030849737
## 19      XY19        4 4.787456e-05 2.856209e-06 5.303099e-01 0.9858489633
## 20      XY20        4 9.208488e-04 1.361713e-04 3.989896e-04 0.9504880905

3.对输入数据的探索

这个包的资料特别少，也没有对应的文章。翻了帮助文档和github，各种教程，没有找到对输入数据的要求说明，暂不知道用的是TPM，FPKM，cpm还是别的什么类型。只能看出没有取log。如果后续有进展，我会在语雀中更新，点击文末原文链接就能看到。

我想试试看log后的表达矩阵是否也可以作为输入

res1 <- callEnsemble(X = log2(Xmat+1), geneids = 'symbol')
res1

##    SampleIDs BestCall            1            2            3            4
## 1        XY1        4 1.121173e-07 1.873900e-06 6.311234e-02 0.9027497470
## 2        XY2        3 4.243225e-02 2.309184e-06 5.495564e-01 0.0084770960
## 3        XY3        4 3.095071e-04 1.029907e-06 1.635270e-01 0.9063920975
## 4        XY4        4 1.154656e-04 3.823049e-07 2.888787e-02 0.8831390142
## 5        XY5        4 1.193054e-07 8.741189e-06 5.060996e-01 0.9260828793
## 6        XY6        4 2.479082e-04 5.528710e-05 1.183165e-04 0.9923758209
## 7        XY7        6 9.421768e-03 1.079501e-03 1.000559e-01 0.0084094275
## 8        XY8        3 6.143126e-04 2.594000e-06 8.872822e-01 0.1306741871
## 9        XY9        4 1.001003e-04 4.479314e-06 3.989228e-03 0.9793768525
## 10      XY10        2 4.624279e-06 9.888145e-01 7.359737e-06 0.0053910189
## 11      XY11        4 5.837736e-05 7.877929e-06 1.970483e-03 0.9895575941
## 12      XY12        3 1.944198e-06 2.005060e-06 3.616153e-01 0.5258071870
## 13      XY13        4 9.002434e-07 9.563100e-04 2.673559e-06 0.9931332767
## 14      XY14        4 1.235332e-05 3.907399e-05 2.785671e-03 0.9972651005
## 15      XY15        4 4.831095e-05 1.127145e-04 1.685999e-04 0.9935480356
## 16      XY16        3 2.731656e-03 1.722274e-05 9.742068e-01 0.0007705171
## 17      XY17        4 3.431090e-06 1.327718e-03 1.570120e-05 0.9971910715
## 18      XY18        3 3.816116e-07 5.000733e-03 8.685111e-01 0.0030849737
## 19      XY19        4 4.787456e-05 2.856209e-06 5.303099e-01 0.9858489633
## 20      XY20        4 9.208488e-04 1.361713e-04 3.989896e-04 0.9504880905

identical(res0$BestCall,res1$BestCall)

## [1] TRUE

两次计算出来的结果一样，那就是可以咯

再试试标准化后的数据行不行

res2 <- callEnsemble(X = t(scale(t(log2(Xmat+1)))), geneids = 'symbol')
res2

##    SampleIDs BestCall            1            2            3            4
## 1        XY1        2 1.055616e-04 9.624696e-01 1.735410e-01 4.088800e-01
## 2        XY2        3 5.200884e-02 2.582083e-01 2.077989e-01 6.790864e-04
## 3        XY3        3 1.050717e-01 8.588433e-05 1.228831e-01 1.287199e-02
## 4        XY4        2 2.210707e-03 1.608360e-01 3.885327e-04 3.327382e-02
## 5        XY5        3 3.972800e-06 9.615958e-03 6.854322e-01 2.701341e-01
## 6        XY6        4 1.715342e-04 3.799206e-01 4.371099e-06 8.692024e-01
## 7        XY7        6 1.250667e-01 1.097162e-01 4.140158e-03 6.115460e-04
## 8        XY8        3 2.792812e-04 2.527673e-04 8.357488e-01 1.747175e-02
## 9        XY9        4 9.675489e-01 8.905539e-02 1.780803e-04 9.530877e-01
## 10      XY10        2 1.756920e-07 9.999766e-01 5.824814e-05 1.829826e-03
## 11      XY11        4 2.118116e-07 1.764702e-02 1.047677e-04 6.112984e-01
## 12      XY12        3 2.059111e-07 1.993342e-01 6.549236e-01 1.061971e-02
## 13      XY13        2 4.458124e-07 9.985121e-01 7.005165e-08 1.838380e-02
## 14      XY14        4 9.115437e-06 2.641302e-03 6.308785e-06 9.922896e-01
## 15      XY15        2 7.724871e-07 9.421750e-01 9.746566e-07 9.982958e-01
## 16      XY16        3 3.725465e-06 2.441265e-03 5.236077e-01 3.568041e-06
## 17      XY17        2 8.232909e-06 9.958895e-01 4.845522e-07 4.058149e-02
## 18      XY18        3 5.211974e-07 3.527417e-01 9.039713e-01 3.747310e-03
## 19      XY19        2 9.518531e-06 1.477052e-01 2.665034e-03 1.926993e-01
## 20      XY20        4 1.538433e-01 1.374537e-01 2.018919e-06 3.616710e-01

identical(res0$BestCall,res2$BestCall)

## [1] FALSE

head(res0$BestCall)

## [1] 4 3 4 4 4 4

head(res2$BestCall)

## [1] 2 3 3 2 3 4

结果不一样了，那就是不行呗。

4.errorreport

区分免疫亚型用到了485个关键的基因。geneMatchErrorReport告诉你有多少基因是在输入数据中不存在的，如果比例过高，那最终计算出来的亚型就不那么可靠。

data(ebpp_gene)
head(ebpp_genes_sig)

##      Symbol Entrez         Ensembl
## 235  ACTL6A     86 ENSG00000136518
## 294   ADAM9   8754 ENSG00000168615
## 305 ADAMTS1   9510 ENSG00000154734
## 322    ADAR    103 ENSG00000160710
## 340   ADCY7    113 ENSG00000121281
## 479   AIMP2   7965 ENSG00000106305

geneMatchErrorReport(X=Xmat, geneid='symbol')

## $matchError
## [1] 0
## 
## $missingGenes
## [1] Symbol  Entrez  Ensembl
## <0 rows> (or 0-length row.names)

5.可视化

原数据是没有分组的，为了体现不同组的亚型比较，这里加上一列分组，仅做例子。

dat = res0
dat$BestCall = factor(dat$BestCall)
dat$group = rep(c("a","b"),each = 10)
library(ggplot2)
library(paletteer)
ggplot(data = dat,aes(x = group))+
  geom_bar(aes(fill = BestCall),
           position = position_dodge2(padding = 0,
                                      preserve = "single"))+
  theme_classic()+
  scale_fill_paletteer_d("RColorBrewer::Set2")

library(dplyr)
dat2 <- dat %>%
 group_by(group, BestCall) %>% 
 summarise(count = n())

ggplot() + geom_bar(data =dat2, aes(x = group, y = count, fill = BestCall),
 stat = "identity",
 position = "fill")+
  theme_classic()+
  scale_fill_paletteer_d("RColorBrewer::Set2")