TCGA 数据分析实战 —— 干性指数（mRNAsi、mDNAsi）
前言
干细胞(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))
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: 这里使用单类模型=>NULLl1:L1-正则化惩罚系数=>0l2: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(cancer) {
query <- GDCquery(
project = cancer,
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "STAR - Counts",
sample.type = c("Primary Tumor"),
)
# 选择 20 个样本
query$results[[1]] <- query$results[[1]][1:20,]
GDCdownload(query)
# 获取 read count
prep <- GDCprepare(
query,
summarizedExperiment = TRUE,
)
tpm <- TCGAanalyze_Preprocessing(
object = prep,
cor.cut = 0.6,
datatype = "tpm_unstrand"
)
exp <- as.data.frame(tpm) %>%
tibble::rownames_to_column("gene_id") %>%
inner_join(dplyr::select(TCGAbiolinks::get.GRCh.bioMart("hg38"), gene_id, gene_name)) %>%
dplyr::select(-gene_id) %>%
group_by(gene_name) %>%
summarise_all(mean) %>%
tibble::column_to_rownames(var = "gene_name") %>%
rename_with(function(x) substr(x, 1, 12))
return(exp)
}
exp <- get_expression("TCGA-BRCA")
exp[1:3,1:3]
# TCGA-3C-AALJ TCGA-A2-A0CO TCGA-A2-A0YM
# 5S_rRNA 0.1292222 0.2757778 0.1316
# 5_8S_rRNA 0.0000000 0.0000000 0.0000
# 7SK 0.0000000 0.0000000 0.0000
导入模型
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(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 甲基化数据使用的是向量点乘。
对于 TCGA 样本的干性指数可以从 TCGA PanCanStemness-2018 网站上下载
TCGAbiolinks 还提供了函数，用于计算各种干性指标。例如
Stemness_score <- TCGAanalyze_Stemness(
stemSig = SC_PCBC_stemSig,
dataGE = exp,
colname.score = "SC_PCBC_stem_score"
)
ECTO_score <- TCGAanalyze_Stemness(
stemSig = ECTO_PCBC_stemSig,
dataGE = exp,
colname.score = "ECTO_PCBC_stem_score"
)
MESO_score <- TCGAanalyze_Stemness(
stemSig = MESO_PCBC_stemSig,
dataGE = exp,
colname.score = "MESO_PCBC_stem_score"
)
扫一扫
748
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
