免疫亚型的鉴定与可视化

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
##               5            6
## 1  5.132385e-02 5.121503e-05
## 2  1.000665e-04 2.261875e-04
## 3  5.731729e-03 1.741630e-04
## 4  1.236814e-03 7.847964e-05
## 5  1.232723e-02 5.358944e-04
## 6  1.019272e-03 2.854635e-04
## 7  3.649302e-05 6.113418e-01
## 8  8.199657e-04 1.142911e-04
## 9  2.182218e-03 1.554600e-04
## 10 3.058250e-05 6.004089e-04
## 11 4.063185e-03 6.252101e-04
## 12 1.320550e-02 2.124778e-04
## 13 1.996231e-04 1.273039e-04
## 14 2.466791e-03 5.618507e-04
## 15 4.465038e-03 2.503012e-04
## 16 4.555504e-05 6.263918e-03
## 17 3.616062e-05 1.134532e-04
## 18 2.343850e-04 2.590704e-03
## 19 2.194258e-02 9.931145e-05
## 20 1.666186e-03 1.220569e-02

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
##               5            6
## 1  5.132385e-02 5.121503e-05
## 2  1.000665e-04 2.261875e-04
## 3  5.731729e-03 1.741630e-04
## 4  1.236814e-03 7.847964e-05
## 5  1.232723e-02 5.358944e-04
## 6  1.019272e-03 2.854635e-04
## 7  3.649302e-05 6.113418e-01
## 8  8.199657e-04 1.142911e-04
## 9  2.182218e-03 1.554600e-04
## 10 3.058250e-05 6.004089e-04
## 11 4.063185e-03 6.252101e-04
## 12 1.320550e-02 2.124778e-04
## 13 1.996231e-04 1.273039e-04
## 14 2.466791e-03 5.618507e-04
## 15 4.465038e-03 2.503012e-04
## 16 4.555504e-05 6.263918e-03
## 17 3.616062e-05 1.134532e-04
## 18 2.343850e-04 2.590704e-03
## 19 2.194258e-02 9.931145e-05
## 20 1.666186e-03 1.220569e-02

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
##               5            6
## 1  0.0970802158 0.0054663994
## 2  0.0022650487 0.0541812675
## 3  0.0191998631 0.0012868469
## 4  0.0070779554 0.0004822279
## 5  0.0083340211 0.0135210226
## 6  0.0007782122 0.0006131011
## 7  0.0004781231 0.1474968009
## 8  0.0134209383 0.0001762331
## 9  0.0017751575 0.0001619612
## 10 0.0004091656 0.0014408990
## 11 0.0018340806 0.0065515421
## 12 0.0795782432 0.0040155074
## 13 0.0009391853 0.0014942663
## 14 0.0016644751 0.0035965345
## 15 0.0084981914 0.0013056903
## 16 0.0001257564 0.1152272336
## 17 0.0021467893 0.0009310142
## 18 0.0043335482 0.0023881147
## 19 0.0033517949 0.0070664277
## 20 0.0082207613 0.0109260329

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

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 195,783评论 5赞 462
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 82,360评论 2赞 373
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 142,942评论 0赞 325
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 52,507评论 1赞 267
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 61,324评论 5赞 358
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 46,299评论 1赞 273
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 36,685评论 3赞 386
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,358评论 0赞 254
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 39,652评论 1赞 293
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 34,704评论 2赞 312
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 36,465评论 1赞 326
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,318评论 3赞 313
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 37,711评论 3赞 299
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 28,991评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,265评论 1赞 251
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 41,661评论 2赞 342
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 40,864评论 2赞 335