前言
干细胞(stem cells
, SC
)是人体内具有自我复制和多向分化潜能的原始细胞,而这种自我更新和分化为成熟细胞的能力称为干性
肿瘤干细胞 (cancer stem cells
,CSCs
) 是指拥有与正常干细胞相关特征的癌细胞,能在特定癌症样本中产生所有细胞类型。通常这类的细胞被认为有形成肿瘤、发展成癌症的潜力,特别是随着癌症的转移,能产生新型的癌症
与不具有干性的肿瘤细胞相比,CSCs
通过 DNA
损伤修复、抑制凋亡通路和产生耐药蛋白等自我保护机制,会导致肿瘤的进展、转移、耐药并增强了自我更新能力
我们要介绍的是一篇发表在 Cell
上的文章所提出的一种机器学习算法(OCLR
,One Class Linear Regression
),对肿瘤样本的干性进行量化。通过对干细胞转录组、甲基化组等多平台数据进行分析,计算出两种不同的干性指数:
-
mRNAsi
:反映干细胞的基因表达特征 -
mDNAsi
:反映干细胞的表观遗传特征
下面我们来介绍这两种干性指数的计算
mRNAsi
1. 数据处理
训练数据使用的是 Progenitor Cell Biology Consortium
(PCBC
)(https://www.synapse.org
) 提供的人类干细胞数据,数据组成包括
我们使用 synapser
包来下载对应的数据,当然也可以手动下载。
首先,安装并导入包
install.packages("synapser", repos = c("http://ran.synapse.org", "http://cran.fhcrc.org"))
library(synapser)
在 https://www.synapse.org/
网站上注册账号,并记住邮箱号密码
在 R
中,调用 synLogin
函数并传入刚才注册的邮箱和密码来登录
synLogin(email = "email@xxx.com", password = "123456")
# Welcome, !NULL
登录成功之后,使用 synGet
函数获取 RNA
表达数据
synRNA <- synGet( "syn2701943", downloadLocation = "~/Downloads/PCBC" )
读入表达数据
library(tidyverse)
exp <- read_delim(file = synRNA$path) %>%
# 去除 Ensembl ID 的后缀
separate(col = "tracking_id", sep = "\\.", into = c("Ensembl_ID", "suffix")) %>%
dplyr::select(-suffix) %>%
column_to_rownames("Ensembl_ID") %>%
as.matrix()
> exp[1:3,1:3]
H9.102.2.5 H9.102.2.6 H9.119.3.7
ENSG00000000003 0.6306521 0.6071539 0.7197784
ENSG00000000005 -1.2838794 -1.0152271 -0.8893850
ENSG00000000419 0.6509436 0.5833117 0.7618753
将 ENSEMBL
转换为 Symbol
library(org.Hs.eg.db)
unimap <- mapIds(
org.Hs.eg.db, keys = rownames(exp), keytype = "ENSEMBL",
column = "SYMBOL", multiVals = "filter"
)
data.exp <- exp[names(unimap),]
rownames(data.exp) <- unimap
使用 synTableQuery
函数来获取元数据,使用 SQL
语法选择 syn3156503
数据中的 UID
和 Diffname_short
列
synMeta <- synTableQuery("SELECT UID, Diffname_short FROM syn3156503")
处理数据
metaInfo <- synMeta$asDataFrame() %>%
dplyr::select(UID, Diffname_short) %>%
column_to_rownames("UID") %>%
filter(!is.na(Diffname_short))
X <- data.exp
y <- metaInfo[colnames(X), ]
names(y) <- colnames(X)
> head(y)
H9.102.2.5 H9.102.2.6 H9.119.3.7 H9.119.5.3 H9.144.7.7 H9EB.558.12.6
"SC" "SC" "SC" "SC" "SC" "EB"
2. 构建模型
在上一步处理完数据之后,可以使用 gelnet
包的 gelnet
函数来构建模型,该函数主要传入 4
个参数
gelnet(X, y, l1, l2)
-
X
: 行为样本,列为基因=>
transpose(X.sc)
-
y
: 这里使用单类模型=>
NULL
-
l1
:L1
-正则化惩罚系数=>
0
-
l2
:L2
-正则化惩罚系数=>
1
# 对数据进行均值中心化
X <- data.exp
m <- apply(X, 1, mean)
X <- X - m
# 将样本分为干细胞组和非干细胞组
sc <- which(y == "SC")
X.sc <- X[, sc]
X.or <- X[, -sc]
model.RNA <- gelnet(X.sc, NULL, 0, 1)
得到每个基因的权重
> head(model.RNA$w)
TSPAN6 TNMD DPM1 SCYL3 C1orf112 FUCA2
6.828043e-05 3.261994e-03 3.170824e-04 -1.860469e-05 1.114508e-03 2.916178e-04
保存模型
save(X, y, model.RNA, file = "~/Downloads/PCBC/model.rda")
3. 模型预测
构建出模型之后,我们可以直接使用该模型来预测其他样本的 mRNAsi
值。
我们首先获取 TCGA
的部分样本
library(TCGAbiolinks)
library(SummarizedExperiment)
get_expression <- function(proj, n = 20) {
query <- GDCquery(
project = proj,
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - FPKM"
)
query$results[[1]] <- query$results[[1]][1:n,]
GDCdownload(query)
data <- GDCprepare(query)
exp <- assay(data) %>% as.data.frame() %>%
rownames_to_column(var = "Ensembl_ID") %>% {
unimap <- mapIds(
org.Hs.eg.db, keys = .$Ensembl_ID, keytype = "ENSEMBL",
column = "SYMBOL", multiVals = "filter")
filter(., Ensembl_ID %in% names(unimap))
} %>%
mutate(Symbol = AnnotationDbi::select(
org.Hs.eg.db, keys = .$Ensembl_ID, keytype = "ENSEMBL",
columns = "SYMBOL")$SYMBOL, .before = Ensembl_ID) %>%
dplyr::select(!Ensembl_ID) %>%
filter(!is.na(Symbol)) %>%
group_by(Symbol) %>%
summarise_all(mean) %>%
column_to_rownames(var = "Symbol")
return(exp)
}
exp <- get_expression("TCGA-BRCA")
> exp[1:3, 1:3]
TCGA-E2-A15G-01A-11R-A12D-07 TCGA-E2-A1B5-01A-21R-A12P-07 TCGA-EW-A2FS-01A-11R-A17B-07
A1BG 0.200303082 0.078878919 0.113499829
A1BG-AS1 1.964474694 0.945018651 0.802335988
A1CF 0.004516686 0.001935602 0.005072972
导入模型
load('~/Downloads/PCBC/model.rda')
提取交叠基因的表达谱及权重
common <- intersect(names(model.RNA$w), rownames(exp))
X <- exp[common, ]
w <- model.RNA$w[common]
对于 RNA
表达数据,使用 spearman
计算权重与表达值之间的相关性来衡量样本的干性指数,并进行标准化使其落在 [0,1]
之间
score <- apply(X, 2, function(z) {cor(z, w, method="sp", use="complete.obs")})
score <- score - min(score)
score <- score / max(score)
样本的干性指数为
> head(score)
TCGA-E2-A15G-01A-11R-A12D-07 TCGA-E2-A1B5-01A-21R-A12P-07 TCGA-EW-A2FS-01A-11R-A17B-07
0.4484917 0.5286787 0.3824672
TCGA-EW-A1P7-01A-21R-A144-07 TCGA-LL-A5YO-01A-21R-A28M-07 TCGA-BH-A1FN-01A-11R-A13Q-07
0.1841952 0.8273317 0.6462018
封装成函数
predict.mRNAsi <- function(exp, modelPath='model.rda') {
load(modelPath)
common <- intersect(names(model.RNA$w), rownames(exp))
X <- exp[common, ]
w <- model.RNA$w[common]
score <- apply(X, 2, function(z) {cor(z, w, method="sp", use="complete.obs")})
score <- score - min(score)
score <- score / max(score)
}
score <- predict.mRNAsi(exp, '~/Downloads/PCBC/model.rda')
mDNAsi
其实 DNA
甲基化干性特征并不是从整个探针范围来寻找的,而是将不同的甲基化特征区域作为输入,根据输入特征的不同,共包含三种方法
-
DMPsi
:全称为differentially methylated probes-based stemness index
,使用差异甲基化探针区域(包含很多过滤条件)作为OCLR
算法的输入,来构建预测模型 -
ENHsi
:enhancer-based stemness index
,将增强子区域的甲基化探针作为OCLR
算法的输入来构建预测模型 -
EREG-mDNAsi
:使用ELMER
包从DNA
甲基化和转录组表达数据重建基因调控网络,将识别出的特征作为OCLR
算法的输入来构建预测模型,由于该包不仅会输出甲基化探针,也包含基因,这部分基因用可以用于构建预测模型,称为EREG-mRNAsi
这三种方法共包含 219
个甲基化探针,将这 219
个探针作为 mDNAsi
的特征。而 EREG-mRNAsi
共包含 103
个基因
1. 构建模型
我们使用文章提供的处理好的甲基化数据来构建模型
> pcbc.data[1:3, 1:3]
SC14_069MESO_3999442133_R01C01 SC14_070EB_3999442133_R01C02 SC14_073MESO_3999442133_R02C02
cg02927655 0.4490163 0.6460517 0.5091612
cg15948871 0.3110161 0.4863665 0.4063254
cg17676824 0.4746514 0.5468089 0.5240562
> pcbc.pd.f[1:3, 1:5]
ROW_ID ROW_VERSION UID biologicalSampleName C4_Cell_Line_ID
13 1986 30 SC11_003_A SC11-003A SC11-003
15 1988 30 SC11_004_A SC11-004A SC11-004
17 1990 30 SC11_014_B SC11-014B SC11-014
load('~/Downloads/PCBC/pcbc.data.Rda')
load('~/Downloads/PCBC/pcbc.pd.f.Rda')
m <- apply(pcbc.data, 1, mean)
pcbc.data.norm <- pcbc.data - m
SC <- pcbc.pd.f[pcbc.pd.f$Diffname_short %in% 'SC',]
X <- pcbc.data.norm[, SC$UID]
model.DNA <- gelnet(t(X), NULL, 0, 1)
save(model.DNA, model.RNA, file = "~/Downloads/PCBC/model-weight.rda")
> length(model.DNA$w)
[1] 219
> head(model.DNA$w)
cg02927655 cg15948871 cg17676824 cg25552705 cg21434327 cg06678662
-0.03251136 -0.03270579 -0.03151675 -0.04281638 -0.02921891 -0.03178537
模型共包含 219
个甲基化探针区域
2. 模型预测
获取 TCGA
甲基化数据
coad.samples <- matchedMetExp("TCGA-COAD", n = 10)
query <- GDCquery(
project = "TCGA-COAD",
data.category = "DNA methylation",
platform = "Illumina Human Methylation 450",
legacy = TRUE,
barcode = coad.samples
)
GDCdownload(query)
met <- GDCprepare(query, save = FALSE)
其实这里还涉及到补缺失值的问题,补缺失值的方法有很多,为了方便,在我直接删除带缺失值的探针
data.met <- subset(met,subset = (rowSums(is.na(assay(met))) == 0))
对于 DNA
数据,使用线性预测模型来计算样本的干性指数
predict.mDNAsi <- function(met, modelPath='model-weight.rda') {
load(modelPath)
common <- intersect(names(model.DNA$w), rownames(met))
X <- met[common, ]
w <- model.DNA$w[common]
score <- t(w) %*% X
score <- score - min(score)
score <- score / max(score)
}
score.meth <- predict.mDNAsi(assay(data.met), "~/Downloads/PCBC/model-weight.rda")
> head(t(score.meth))
[,1]
TCGA-4T-AA8H-01A-11D-A40X-05 1.0000000
TCGA-AZ-4614-01A-01D-1407-05 0.7734937
TCGA-A6-2675-01A-02D-1721-05 0.3427734
TCGA-A6-6781-01A-22D-A27A-05 0.1482315
TCGA-A6-6781-01B-06D-A27A-05 0.1319317
TCGA-A6-6781-01A-22D-1926-05 0.0000000
总结
这些模型已经全部构建了,可以从
https://github.com/dxsbiocc/learn/tree/main/R/CSCs
中下载 Stemness_index.rda
文件,获取这些模型的特征和权重值
计算其他样本的干性指数也很简单,RNA-seq
数据使用的是 spearman
相关,DNA
甲基化数据使用的是向量点乘