R画图技术专题good code

TCGA 数据分析实战 —— 干性指数（mRNAsi、mDNAs

2021-08-10  本文已影响0人  名本无名

前言

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

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

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

我们要介绍的是一篇发表在 Cell 上的文章所提出的一种机器学习算法（OCLR，One Class Linear Regression），对肿瘤样本的干性进行量化。通过对干细胞转录组、甲基化组等多平台数据进行分析，计算出两种不同的干性指数：

下面我们来介绍这两种干性指数的计算

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 <- 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 甲基化干性特征并不是从整个探针范围来寻找的，而是将不同的甲基化特征区域作为输入，根据输入特征的不同，共包含三种方法

这三种方法共包含 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 甲基化数据使用的是向量点乘

上一篇 下一篇

猜你喜欢

热点阅读