前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >基于支持向量机模型的TNBC的分子亚型预测

基于支持向量机模型的TNBC的分子亚型预测

作者头像
生信技能树
发布2022-06-27 21:01:58
7280
发布2022-06-27 21:01:58
举报
文章被收录于专栏:生信技能树

连续两次求贤令:曾经我给你带来了十万用户,但现在祝你倒闭,以及 生信技能树知识整理实习生招募,让我走大运结识了几位优秀小伙伴!有做ngs实战整理的,也有做临床数据挖掘算法工具介绍的。今天分享的是复旦大学和西北民族大学小伙伴合作的笔记

下面是笔记原文

1. 引入

最近刚好毕业,闲来无事,偶然看见曾老师的知识整理的邀请,故应邀加入知识整理的队伍。TNBC.CMS [1] R包是曾老师安排的第一个整理任务。

TNBC.CMS: prediction of TNBC consensus molecular subtype TNBC共识分子亚型的预测

[1] Kim J, Yu D, Kwon Y, et al. Genomic Characteristics of Triple-Negative Breast Cancer Nominate Molecular Subtypes That Predict Chemotherapy Response. Mol Cancer Res. 2020;18(2):253-263. doi:10.1158/1541-7786.MCR-19-0453

2. TNBC.CMS 开发背景

TNBC.CMS是一个用于三阴性乳腺癌(TNBC)分子亚型分类的包。虽然已有了很多分类策略,但缺乏准确的亚型分类方法仍是患者诊断和TNBC研究的局限。作者的基于机器学习的分类器模型使用957名TNBC患者的基因表达谱。

TNBC.CMS包将患者分为四种常见的共识分子亚型(consensus molecular subtypes,CMS):mesenchymal-like(MSL)、immunomodulatory(IM)、luminal AR(LAR)和stem-like(SL)。它还提供了基因组和临床特征的概要,包括生存率、风险比、通路活性和药物反应

后续原作者在2021年基于前述的TNBC四种亚型,探索了根据 SMARCA4活性调节的通路及Cisplatin Resistance Score顺铂耐药评分,发表于Cancers[2]。(这个蛮有意思的,值得读)

[2] Kim J, Jang G, Sim SH, Park IH, Kim K, Park C. SMARCA4 Depletion Induces Cisplatin Resistance by Activating YAP1-Mediated Epithelial-to-Mesenchymal Transition in Triple-Negative Breast Cancer. Cancers (Basel). 2021;13(21):5474. Published 2021 Oct 30. doi:10.3390/cancers13215474

在初步了解TNBC.CMS的开发背景后,我也对TNBC.CMS产生了一些兴趣,我好奇TNBC.CMS的分类方法是否具有普及性,能不能够其他肿瘤中进行推广,完成对其他肿瘤的分型。

下面将尝试重复其分类步骤,并尝试加以解读。

3. TNBC.CMS 安装

作者将TNBC.CMS存储在Biocondutor, 因此我们可以使用以下的方式安装:

代码语言:javascript
复制
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("TNBC.CMS")

如想直接查看作者给定的说明文档,可以参考以下代码:

代码语言:javascript
复制
browseVignettes("TNBC.CMS")

4. TNBC.CMS使用

下面根据作者提供的说明文档,进行试用

4.1 加载R包

代码语言:javascript
复制
library(TNBC.CMS)
## 如不想看到过多加载信息 可使用
## suppressMessages(library(TNBC.CMS))

4.2 加载数据

代码语言:javascript
复制
data("GSE25055")  #获取R包中提供的示例数据
dim(assays(GSE25055)[[1]]) #查看数据维度 dim == Dimension 
## [1] 4746 73  ##表明 有4746个genes 和76个 samples
assays(GSE25055)[[1]][1:5, 1:5]  ## 查看数据内容
##        615397   615396   615394   615393   615392
## HBA1  15.42644 16.03556 16.01673 15.33820 12.66187
## HBB   14.91473 15.70969 15.72610 14.80730 12.30998
## B2M   12.91446 13.98487 13.79634 12.62554 11.59076
## FTH1  11.98103 12.57083 12.08568 12.98890 12.16884
## HLA-B 12.49575 14.57368 14.37568 12.32062 10.64294
class(GSE25055)    ## 查看数据的类型
## [1] "SummarizedExperiment"
attr(,"package")
## [1] "SummarizedExperiment"
View(GSE25055)

可以看到TNBC.CMS 提供的是一个 “SummarizedExperiment” 格式的数据类型,在assays(GSE25055)[1]中,行名是gene symbol,列名是样本名,数据为表达值。

4.3 示例1:CMS分类(基于共识分子亚型)

`precdictCMS`函数根据输入矩阵或`SummarizedExperiment`对象将共识分子亚型分配给TNBC样本。如果输入是一个`SummarizedExperiment`对象,分析列表中的第一个元素应该是基因表达矩阵。在任何情况下,基因表达谱都不应该被缩放(原文:Scaled )或对数转换 (原文:log-transformed )。可以通过访问`probabilities`属性来检索分类概率。

代码语言:javascript
复制
predictions <- predictCMS(GSE25055)
## Warning message:
## In class(object) <- "environment" :
##   把class(x)设成"environment"集空属性; 其结果不再是S4目标对象
table(predictions)
## predictions
## MSL  IM LAR  SL 
## 12  14  23  24 
head(attr(predictions, "probabilities"))
##                MSL          IM         LAR         SL
## 615397 0.004965739 0.003456084 0.006844380 0.98473380
## 615396 0.026430368 0.014144690 0.939662186 0.01976276
## 615394 0.056331795 0.020242679 0.020387143 0.90303838
## 615393 0.001031452 0.001748296 0.956579694 0.04064056
## 615392 0.002081453 0.038822011 0.791781322 0.16731521
## 615390 0.029205717 0.008909284 0.004956605 0.95692839

这里四个分类即是前述的几个亚群:

间充质样干细胞样亚群(mesenchymal-like stem-like,MSL)

免疫调节亚群( immunomodulatory,IM)

管腔雄激素受体亚群(luminal androgen receptor,LAR)

干细胞样亚群(stem-like,SL)

实际上到这里根本看不明白这四个群是怎么来的,只知道,我丢了一个表达矩阵进去,然后这个函数就自动将各个样本分别属于各个群的概率算了出来。按上述结果描述,编号“615397”属于SL的概率为 ‘0.98473380’,因此判定为 SL。

因此为了更好的了解这一步干了什么,有必要查看一下 “predictCMS” 函数的源代码

代码语言:javascript
复制
predictCMS
function (expr) 
{
    if (is(expr, "SummarizedExperiment")) {
        exp.mat <- assays(expr)[[1]]
    }
    else {
        exp.mat <- expr
    }
    missings <- genelist[!genelist %in% rownames(exp.mat)]
    nmissing <- length(missings)
    if (nmissing > 0) {
        missing.mat <- matrix(NA, nrow = nmissing, ncol = ncol(exp.mat))
        rownames(missing.mat) <- missings
        colnames(missing.mat) <- colnames(exp.mat)
        input <- rbind(exp.mat, missing.mat)
        input <- impute(input, what = "median")
    }
    else {
        input <- exp.mat
    }
    input <- log2(input + 1)
    input <- input - apply(input, 1, median)
    input <- data.frame(t(input[genelist, ]), check.names = FALSE)
    pred <- predict(SVM.model, input, probability = TRUE)
    prob <- attr(pred, "probabilities")
    attr(pred, "probabilities") <- prob[, c("MSL", 
        "IM", "LAR", "SL")]
    return(pred)
}
<bytecode: 0x000001b9473dc4a8>
<environment: namespace:TNBC.CMS>

可以看到,出现“input”变量前面的步骤均是在对数据进行清洗处理,后续的关键的分类步骤是使用SVM.model完成的,然后就直接出现了"MSL", "IM", "LAR", "SL" 几个我们关注的因素, 而前面的 genelist 和 SVM.model 都不是在我们环境下输入的参数。

代码语言:javascript
复制
genelist
## 错误: 找不到对象'genelist'
SVM.model
## 错误: 找不到对象'SVM.model'

因此,合理的推测是 作者先将给定genelist 在这个表达矩阵中筛选出来,然后通过已构建好的SVM.model将各个样本进行预测,而这个genelist和SVM.model均被封装在这个R包,在运行这个函数的时候调用了给定的genelist 和 SVM.model。

那下一步我想试试能不能在R包给的文档里面找到这个genelist 和SVM.model。

但……我没找到,先这样吧,就不浪费过多时间了,希望有同学会的可以指导下,想学习这个genelist。

(小编寄语:其实超级简单, 使用3个冒号即可,代替之前的两个冒号的操作符。)

其实就是实习生的编程基础知识不过关,具备基础的计算机知识非常重要,我把它粗略的分成基于R语言的统计可视化,以及基于Linux的NGS数据处理

但大概意思估摸就是各种给定每组类型匹配一个 genelist,然后根据给定基因的表达量在SVM.model下去预测这个样本分别属于哪个类型的。

4.4 示例2:基因组和临床特征的概要

`TNBC.CMS`包有几个函数用于研究共识分子亚型的基因组和临床特征。在本节中,我们将这些函数应用于GSE25055基因表达和临床特征的数据集。

4.4.1 ComputeGES

`ComputeGES`函数计算以下7个基因表达特征的特征分数:EMT(epithelial-mesencymal transition上皮-间充质转化)、stromal基质、免疫、微环境、stemness干性、hormone激素和CIN(chro- mosomal instability染色体不稳定性)。有关基因表达特征的更多详细信息,请参见`ComputeGES`函数手册页面。

代码语言:javascript
复制
resultGES <- computeGES(expr = GSE25055, pred = predictions,
+                         rnaseq = FALSE)
## Warning message:
## In .gsva(expr, mapped.gset.idx.list, method, kcdf, rnaseq, abs.ranking,  :
##  Some gene sets have size one. Consider setting 'min.sz > 1'.
resultGES[,1:4]
##                       615397       615396       615394       615393
## EMT                 6.0372043    6.1889247    6.3529614    5.7296889
## Stromal          -605.4905027 -642.1497283 -552.4131551 -929.4231951
## Immune            167.0953013  282.4093345  293.8429060 -141.3963687
## Microenvironment    0.3019017    0.2577112    0.3392774    0.2170629
## Stemness            0.5490508    0.3003425    0.3994718    0.6053706
## Hormone             4.9628294    4.9485363    5.5935250    5.4526575
## CIN                 7.3907780    6.7571045    7.2265374    7.5212038

在运行的过程中,同时生成了如下所示图形

如图1所示,此函数绘制了特征分数的箱线图,并用亚型之间的p值进行比较。MSL亚型、IM亚型、LAR亚型和SL亚型的Stromal基质、immune免疫、hormone激素和stemness干性得分分别显著高于其他亚型。

照例,我们也来看一下这个函数进行了什么操作

代码语言:javascript
复制
computeGES
function (expr, pred, rnaseq = FALSE) 
{
    if (is(expr, "SummarizedExperiment")) {
        exp.mat <- assays(expr)[[1]]
    }
    else {
        exp.mat <- expr
    }
    CMS.palette <- c(MSL = "brown2", IM = "gold2", 
        LAR = "yellowgreen", SL = "midnightblue")
    pred <- factor(pred, levels = c("MSL", "IM", 
        "LAR", "SL"))
    EMT.score <- colMeans(exp.mat[rownames(exp.mat) %in% EMT.geneset, 
        ])
    CIN.score <- colMeans(exp.mat[rownames(exp.mat) %in% CIN.geneset, 
        ])
    hormone.score <- colMeans(exp.mat[rownames(exp.mat) %in% 
        c("AR", "ESR1", "ERBB2", "PGR"), 
        ])
    X <- exp.mat[rownames(exp.mat) %in% names(Stemness.geneset), 
        ]
    Stemness.geneset = Stemness.geneset[rownames(X)]
    s <- apply(X, 2, function(z) {
        cor(z, Stemness.geneset, method = "sp", use = "complete.obs")
    })
    s <- s - min(s)
    stemness.score <- s/max(s)
    ESTIMATE.score <- computeESTIMATEscore(exp.mat)
    stromal.score <- unlist(ESTIMATE.score["StromalSignature", 
        ])
    immune.score <- unlist(ESTIMATE.score["ImmuneSignature", 
        ])
    microenvironment.score <- computexCellScore(exp.mat, rnaseq)
    sig.mat <- rbind(EMT.score, stromal.score, immune.score, 
        microenvironment.score, stemness.score, hormone.score, 
        CIN.score)
    rownames(sig.mat) <- c("EMT", "Stromal", "Immune", 
        "Microenvironment", "Stemness", "Hormone", 
        "CIN")
 ……
  ## 后续计算P值及画图部分省略

可以发现,这个函数的评分方式很简单粗暴,比如EMT.score评分,实际上就是选取EMT相关基因,然后在各个样本中取ColMeans 即作为他的评分,个人感觉还不如GSVA评分或者用AUCell进行排序,单细胞的addmodule 感觉也比这个靠谱。免疫基质相关评分采用的是 Estimate 的比较经典的方式。但这里有个问题是,在函数运行的时候 出现一个gsva函数的warning,但是函数代码没有看到gsva,也可能是我没有理解到位,持保留意见。

代码语言:javascript
复制
EMT.score <- colMeans(exp.mat[rownames(exp.mat) %in% EMT.geneset,])
4.4.2 performGSVA

`performGSVA`函数对基因集进行基因集差异分析,并生成代表GSVA富集分数的热图。如果基因组没有给出,则使用标记通路基因组。用户还可以通过设置`gsva.kcdf`参数来选择用于估计表达式值的累积分布函数的要点,该参数默认设置为"`Gaussian`"。如果表达值是整数,建议使用"`Poisson`"。

代码语言:javascript
复制
resultGSVA <- performGSVA(expr = GSE25055, pred = predictions,
gene.set = NULL)
head(resultGSVA[,1:4])

这部分就不展开了,就是一个常见的GSVA的富集分析。

4.4.3 plotKM 和 plotHR

`TNBC.CMS`包提供了两个生存分析函数:`plotKM`和`plotHR`。在这里,我们使用来自GSE25055数据集的生存数据来研究总生存率和共识分子亚型之间的关系。生存数据也包含在`SummarizedExperiment`对象中,可以使用`colData`函数进行访问。

代码语言:javascript
复制
time <- colData(GSE25055)$DFS.month
event <- colData(GSE25055)$DFS.status
plotKM(pred = predictions, time = time, event = event)

`plotKM`函数为每个共识分子亚型生成Kaplan-Meier曲线,如图3所示。SL组的预后最差,这与作者前述研究是一致的。

代码语言:javascript
复制
> library(survival)
> #Test for difference of survival between low and high expression groups
> surv <- Surv(time, event)
> GSE25055.exprs <- assays(GSE25055)[[1]]
> chisq <- apply(GSE25055.exprs, 1, function(x) survdiff(surv ~ (x > median(x)))$chisq)
> pval <- 1 - pchisq(chisq, 1)
> #Select 10 genes with lowest p-values for the log-rank test
> gs <- names(sort(pval)[1:10])
> gs
 [1] "RECK"   "RELN"   "EHD4"   "PRRX2"  "FOLR1"  "UGCG"   "GOT2"   "PRPF39" "OPHN1"  "CPE"   
> ## [1] "RECK" "RELN" "EHD4" "PRRX2" "FOLR1" "UGCG" "GOT2" "PRPF39"
> ## [9] "OPHN1" "CPE"
> plotHR(expr = GSE25055, gene.symbol = gs, pred = predictions, time = time,
+        event = event, by.subtype = FALSE)

也可以设置 “by.subtype= TRUE“ 按亚群进行计算HR

代码语言:javascript
复制
plotHR(expr = GSE25055, gene.symbol = gs[1:4], pred = predictions,
time = time, event = event, by.subtype = TRUE)

4.5 示例3:药物反应研究

`TNBC.CMS` 包还提供了预测药物反应的函数。`computeDS` 函数计算 `MSigDB CGP`(chemical and genetic perturbations化学和遗传扰动)集合中相应基因集的药物特征分数,并绘制特征分数的热图。药物特征分数计算为基因集平均表达值之间的差值,这些基因集与药物反应和耐药性相关的。分数越高,患者越有可能引起反应。用户可以通过`gene.set` 参数提供他们自己的基因集。请注意,基因组的名称必须遵循 `[DRUG NAME]_[RESPONSE/RESISTANCE]_[UP/DN]` 的格式(例如 `CISPLATIN_RESISTANCE_UP`)。

图6显示了每个样本药物特征分数的热图。MSL和SL亚型似乎分别对dasatinib 达沙替尼 和doxorubicin阿霉素耐药。此外,IM和LAR亚型对 androgen agonist雄激素激动剂和SB216763(GSK3B的抑制剂)的得分高于其他亚型。

代码语言:javascript
复制
resultDS <- computeDS(expr = GSE25055, pred = predictions)
head(resultDS[,1:4])
##                 615397     615396     615394     615393
##  APLIDIN      -0.8876518 -0.4680677 -0.8102091 -0.8549382
##  CISPLATIN     0.9589450  0.6949340  1.2613627  0.2520613
##  DASATINIB     0.5460543  0.7005523 -0.1239162  0.7189546
##  FORSKOLIN     6.4843470  6.7434440  6.4589256  6.6777797
##  IMATINIB     -6.3512657 -5.0914031 -5.4834761 -5.2843840
##  TROGLITAZONE  6.8838127  6.7156754  6.7170530  6.6795939

5. 结果保存

以上是TNBC.CMS提供的几个函数的主要用法,最后输出的结果作者也提供了一个保存方式,以供后续使用。

代码语言:javascript
复制
dfCMS <- data.frame(row.names = colnames(GSE25055.exprs), CMS = predictions, t(resultGES),
t(resultDS), stringsAsFactors = FALSE)
head(dfCMS)
##       CMS      EMT    Stromal     Immune Microenvironment  Stemness  Hormone      CIN    APLIDIN CISPLATIN  DASATINIB FORSKOLIN  IMATINIB
## 615397  SL 6.037204  -605.4905  167.09530        0.3019017 0.5490508 4.962829 7.390778 -0.8876518 0.9589450  0.5460543  6.484347 -6.351266
## 615396 LAR 6.188925  -642.1497  282.40933        0.2577112 0.3003425 4.948536 6.757104 -0.4680677 0.6949340  0.7005523  6.743444 -5.091403
## 615394  SL 6.352961  -552.4132  293.84291        0.3392774 0.3994718 5.593525 7.226537 -0.8102091 1.2613627 -0.1239162  6.458926 -5.483476
## 615393 LAR 5.729689  -929.4232 -141.39637        0.2170629 0.6053706 5.452658 7.521204 -0.8549382 0.2520613  0.7189546  6.677780 -5.284384
## 615392 LAR 5.576038 -1518.6248 -560.49739        0.1650646 1.0000000 4.224790 8.239136 -1.3667877 0.3644356  0.8224894  6.703497 -4.930696
## 615390  SL 6.233911  -607.0961   31.71176        0.2690540 0.5903492 4.602222 7.540700 -0.8707087 0.7678301  0.1427109  6.795645 -5.007474
##        TROGLITAZONE      TZD SB216763 ANDROGEN TAMOXIFEN BEXAROTENE DOXORUBICIN
## 615397     6.883813 7.214105 6.142319 3.949582 -6.029159  0.3264693   -6.875380
## 615396     6.715675 5.341876 6.268575 4.039061 -6.283494  0.1351034   -6.320782
## 615394     6.717053 6.388377 5.994192 3.848465 -6.292139 -1.0281784   -6.709809
## 615393     6.679594 6.290939 5.844824 3.895471 -5.917298  0.1482893   -6.984418
## 615392     7.042860 5.494772 5.731897 4.047578 -6.097761 -0.1288242   -7.581800
## 615390     6.856523 7.490999 6.045182 4.173174 -6.205849 -0.8771287   -6.787671
代码语言:javascript
复制
6. TNBC.CMS试用总结

在最开始拿到这个TNBC.CMS的时候,我以为他是一个无监督聚类生成四个分组,然后作者根据各个分组的基因表达模式分为四个亚群。但是在逐步探索后发现,实际上这个R包还是一个根据既定分组,既定亚型相关基因做出的区分,就相当于,我先划好4个区域,再把这些样本根据不同区域的特征分选出来,分选用的方法就是SVM。因此,虽然R包命名为TNBC.CMS,但实际上和TNBC没啥关系,换个肿瘤,套入这个模式,大概率也能分成四组。

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2022-06-01,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1. 引入
  • 2. TNBC.CMS 开发背景
  • 3. TNBC.CMS 安装
  • 4. TNBC.CMS使用
    • 4.1 加载R包
      • 4.2 加载数据
        • 4.3 示例1:CMS分类(基于共识分子亚型)
          • 4.4 示例2:基因组和临床特征的概要
            • 4.4.1 ComputeGES
            • 4.4.2 performGSVA
            • 4.4.3 plotKM 和 plotHR
          • 4.5 示例3:药物反应研究
          • 5. 结果保存
          领券
          问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档