TCGA 数据分析实战 —— 干性指数(mRNAsi、mDNAsi)

前言

干细胞(stem cells, SC)是人体内具有自我复制和多向分化潜能的原始细胞,而这种自我更新和分化为成熟细胞的能力称为干性

肿瘤干细胞 (cancer stem cells,CSCs) 是指拥有与正常干细胞相关特征的癌细胞,能在特定癌症样本中产生所有细胞类型。通常这类的细胞被认为有形成肿瘤、发展成癌症的潜力,特别是随着癌症的转移,能产生新型的癌症

与不具有干性的肿瘤细胞相比,CSCs 通过 DNA 损伤修复、抑制凋亡通路和产生耐药蛋白等自我保护机制,会导致肿瘤的进展、转移、耐药并增强了自我更新能力

我们要介绍的是一篇发表在 Cell 上的文章所提出的一种机器学习算法(OCLROne 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 数据中的 UIDDiffname_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 算法的输入,来构建预测模型
  • ENHsienhancer-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 甲基化数据使用的是向量点乘

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,222评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,455评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,720评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,568评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,696评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,879评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,028评论 3 409
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,773评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,220评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,550评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,697评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,360评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,002评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,782评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,010评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,433评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,587评论 2 350

推荐阅读更多精彩内容