TCGA 數(shù)據(jù)分析實戰(zhàn) —— 干性指數(shù)（mRNAsi、mDNAsi）

前言

干細(xì)胞(stem cells, SC)是人體內(nèi)具有自我復(fù)制和多向分化潛能的原始細(xì)胞，而這種自我更新和分化為成熟細(xì)胞的能力稱為干性

腫瘤干細(xì)胞 (cancer stem cells,CSCs) 是指擁有與正常干細(xì)胞相關(guān)特征的癌細(xì)胞，能在特定癌癥樣本中產(chǎn)生所有細(xì)胞類型。通常這類的細(xì)胞被認(rèn)為有形成腫瘤、發(fā)展成癌癥的潛力，特別是隨著癌癥的轉(zhuǎn)移，能產(chǎn)生新型的癌癥

與不具有干性的腫瘤細(xì)胞相比，CSCs 通過 DNA 損傷修復(fù)、抑制凋亡通路和產(chǎn)生耐藥蛋白等自我保護機制，會導(dǎo)致腫瘤的進展、轉(zhuǎn)移、耐藥并增強了自我更新能力

我們要介紹的是一篇發(fā)表在 Cell 上的文章所提出的一種機器學(xué)習(xí)算法（OCLR，One Class Linear Regression），對腫瘤樣本的干性進行量化。通過對干細(xì)胞轉(zhuǎn)錄組、甲基化組等多平臺數(shù)據(jù)進行分析，計算出兩種不同的干性指數(shù)：

  • mRNAsi：反映干細(xì)胞的基因表達特征
  • mDNAsi：反映干細(xì)胞的表觀遺傳特征

下面我們來介紹這兩種干性指數(shù)的計算

mRNAsi

1. 數(shù)據(jù)處理

訓(xùn)練數(shù)據(jù)使用的是 Progenitor Cell Biology Consortium (PCBC)(https://www.synapse.org) 提供的人類干細(xì)胞數(shù)據(jù)，數(shù)據(jù)組成包括

我們使用 synapser 包來下載對應(yīng)的數(shù)據(jù)，當(dāng)然也可以手動下載。

首先，安裝并導(dǎo)入包

install.packages("synapser", repos = c("http://ran.synapse.org", "http://cran.fhcrc.org"))

library(synapser)

在 https://www.synapse.org/ 網(wǎng)站上注冊賬號，并記住郵箱號密碼

在 R 中，調(diào)用 synLogin 函數(shù)并傳入剛才注冊的郵箱和密碼來登錄

synLogin(email = "email@xxx.com", password = "123456")
# Welcome,  !NULL

登錄成功之后，使用 synGet 函數(shù)獲取 RNA 表達數(shù)據(jù)

synRNA <- synGet( "syn2701943", downloadLocation = "~/Downloads/PCBC" )

讀入表達數(shù)據(jù)

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 轉(zhuǎn)換為 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 函數(shù)來獲取元數(shù)據(jù)，使用 SQL 語法選擇 syn3156503 數(shù)據(jù)中的 UID 和 Diffname_short 列

synMeta <- synTableQuery("SELECT UID, Diffname_short FROM syn3156503")

處理數(shù)據(jù)

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. 構(gòu)建模型

在上一步處理完數(shù)據(jù)之后，可以使用 gelnet 包的 gelnet 函數(shù)來構(gòu)建模型，該函數(shù)主要傳入 4 個參數(shù)

gelnet(X, y, l1, l2)
  • X: 行為樣本，列為基因 => transpose(X.sc)
  • y: 這里使用單類模型 => NULL
  • l1: L1-正則化懲罰系數(shù) => 0
  • l2: L2-正則化懲罰系數(shù) => 1 
# 對數(shù)據(jù)進行均值中心化
X <- data.exp
m <- apply(X, 1, mean)
X <- X - m
# 將樣本分為干細(xì)胞組和非干細(xì)胞組
sc <- which(y == "SC")
X.sc <- X[, sc]
X.or <- X[, -sc]

model.RNA <- gelnet(X.sc, NULL, 0, 1)

得到每個基因的權(quán)重

> 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. 模型預(yù)測

構(gòu)建出模型之后，我們可以直接使用該模型來預(yù)測其他樣本的 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

導(dǎo)入模型

load('~/Downloads/PCBC/model.rda')

提取交疊基因的表達譜及權(quán)重

common <- intersect(names(model.RNA$w), rownames(exp))
X <- exp[common, ]
w <- model.RNA$w[common]

對于 RNA 表達數(shù)據(jù)，使用 spearman 計算權(quán)重與表達值之間的相關(guān)性來衡量樣本的干性指數(shù)，并進行標(biāo)準(zhǔn)化使其落在 [0,1] 之間

score <- apply(X, 2, function(z) {cor(z, w, method="sp", use="complete.obs")})
score <- score - min(score)
score <- score / max(score)

樣本的干性指數(shù)為

> 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

封裝成函數(shù)

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 甲基化干性特征并不是從整個探針范圍來尋找的，而是將不同的甲基化特征區(qū)域作為輸入，根據(jù)輸入特征的不同，共包含三種方法

  • DMPsi：全稱為 differentially methylated probes-based stemness index，使用差異甲基化探針區(qū)域（包含很多過濾條件）作為 OCLR 算法的輸入，來構(gòu)建預(yù)測模型
  • ENHsi：enhancer-based stemness index，將增強子區(qū)域的甲基化探針作為 OCLR 算法的輸入來構(gòu)建預(yù)測模型
  • EREG-mDNAsi：使用 ELMER 包從 DNA 甲基化和轉(zhuǎn)錄組表達數(shù)據(jù)重建基因調(diào)控網(wǎng)絡(luò)，將識別出的特征作為 OCLR 算法的輸入來構(gòu)建預(yù)測模型，由于該包不僅會輸出甲基化探針，也包含基因，這部分基因用可以用于構(gòu)建預(yù)測模型，稱為 EREG-mRNAsi

這三種方法共包含 219 個甲基化探針，將這 219 個探針作為 mDNAsi 的特征。而 EREG-mRNAsi 共包含 103 個基因

1. 構(gòu)建模型

我們使用文章提供的處理好的甲基化數(shù)據(jù)來構(gòu)建模型

> 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 個甲基化探針區(qū)域

2. 模型預(yù)測

獲取 TCGA 甲基化數(shù)據(jù)

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 數(shù)據(jù)，使用線性預(yù)測模型來計算樣本的干性指數(shù)

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

總結(jié)

這些模型已經(jīng)全部構(gòu)建了，可以從
https://github.com/dxsbiocc/learn/tree/main/R/CSCs 中下載 Stemness_index.rda 文件，獲取這些模型的特征和權(quán)重值

計算其他樣本的干性指數(shù)也很簡單，RNA-seq 數(shù)據(jù)使用的是 spearman 相關(guān)，DNA 甲基化數(shù)據(jù)使用的是向量點乘

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容