前面我们提到过一个bioconductor包居然发在了《cancer research》杂志，该文章 标题是：《Software for the Integration of Multiomics Experiments in Bioconductor》，文章链接是：https://cancerres.aacrjournals.org/content/77/21/e39
其bioconductor 链接是：https://bioconductor.org/packages/release/bioc/html/MultiAssayExperiment.html
主要是定义了一个数据结构（S4对象），如下所示：
MultiAssayExperiment结构
S4对象基本上是R语言分水岭了，无论是理解并且使用它还是创造它，都是一个门槛，甚至我在讲解单细胞数据分析流程的时候，把S4对象的理解作为了基本功！
安装和加载R包相信已经无需我多说了：
代码语言：javascript复制 BiocManager::install("curatedTCGAData")
BiocManager::install("TCGAutils")
library(curatedTCGAData)
library(MultiAssayExperiment)
library(TCGAutils)
首先查看TCGA数据库有哪些数据
代码如下：
代码语言：javascript复制data('diseaseCodes', package = "TCGAutils")
head(diseaseCodes)
可以看到是居然是37个数据集，但是有分子分型的并不多：
37个数据集
主要的癌症简称，以及其各个癌症内部的数据集形式，如下所示：
代码语言：javascript复制
Available Cancer codes:
ACC BLCA BRCA CESC CHOL COAD DLBC ESCA GBM HNSC KICH
KIRC KIRP LAML LGG LIHC LUAD LUSC MESO OV PAAD PCPG
PRAD READ SARC SKCM STAD TGCT THCA THYM UCEC UCS UVM
Available Data Types:
CNACGH CNACGH_CGH_hg_244a
CNACGH_CGH_hg_415k_g4124a CNASeq CNASNP
CNVSNP GISTIC_AllByGene GISTIC_Peaks
GISTIC_ThresholdedByGene Methylation
Methylation_methyl27 Methylation_methyl450
miRNAArray miRNASeqGene mRNAArray
mRNAArray_huex mRNAArray_TX_g4502a
mRNAArray_TX_g4502a_1
mRNAArray_TX_ht_hg_u133a Mutation
RNASeq2GeneNorm RNASeqGene RPPAArray
多种癌症，多种数据类型。也是需要联网下载数据，但是可以使用 dry.run 控制是否真的下载，因为如果是下载甲基化信号值矩阵或者表达量矩阵，会耗时很长。
获取突变信息和拷贝数信息
因为下载甲基化信号值矩阵或者表达量矩阵，会耗时很长，所以我们拿简单的突变信息和拷贝数信息举例。
代码语言：javascript复制(accmae <- curatedTCGAData(
"ACC", c("CN*", "Mutation"), version = "2.0.1", dry.run = FALSE
))
accmae
得到一个MultiAssayExperiment结构的变量，如下所示：
代码语言：javascript复制> accmae
A MultiAssayExperiment object of 3 listed
experiments with user-defined names and respective classes.
Containing an ExperimentList class object of length 3:
[1] ACC_CNASNP-20160128: RaggedExperiment with 79861 rows and 180 columns
[2] ACC_CNVSNP-20160128: RaggedExperiment with 21052 rows and 180 columns
[3] ACC_Mutation-20160128: RaggedExperiment with 20166 rows and 90 columns
Functionality:
experiments() - obtain the ExperimentList instance
colData() - the primary/phenotype DataFrame
sampleMap() - the sample coordination DataFrame
`$`, `[`, `[[` - extract colData columns, subset, or experiment
*Format() - convert into a long or wide DataFrame
assays() - convert ExperimentList to a SimpleList of matrices
exportClass() - save all data to files
可以看到我们前面的 参数 c("CN*", "Mutation")， 设置了可以下载一个突变信息文件，和两个拷贝数信息文件。因为在线下载，所以仍然是取决于网络速度哦。
探索临床信息
前面的MultiAssayExperiment结构的变量，可以被很多函数来处理，包括 getSubtypeMap 和 sampleTables。我们的突变信息和两个拷贝数信息，都需要有临床信息才能进行生物学意义探索。
如下所示：
代码语言：javascript复制> head(getSubtypeMap(accmae))
ACC_annotations ACC_subtype
1 Patient_ID patientID
2 histological_subtypes Histology
3 mrna_subtypes C1A/C1B
4 mrna_subtypes mRNA_K4
5 cimp MethyLevel
6 microrna_subtypes miRNA cluster
> sampleTables(accmae)
$`ACC_CNASNP-20160128`
01 10 11
90 85 5
$`ACC_CNVSNP-20160128`
01 10 11
90 85 5
$`ACC_Mutation-20160128`
01
90
可以看到一个突变信息文件里面仅仅是记录了ACC这个癌症的全部的病人，而两个拷贝数信息文件是包含了ACC癌症的病人的正常组织和肿瘤组织。如果你对 01 和 10 不了解，也可以看看备忘录。
代码语言：javascript复制data(sampleTypes, package = "TCGAutils")
head(sampleTypes)
sampleTypes[sampleTypes[["Code"]] %in% c("01", "10"), ]
可以看到：
代码语言：javascript复制> sampleTypes[sampleTypes[["Code"]] %in% c("01", "10"), ]
Code Definition Short.Letter.Code
1 01 Primary Solid Tumor TP
10 10 Blood Derived Normal NB
前面的 01 10 11 代表 是否是肿瘤样品，就可以提取肿瘤相关数据：
代码语言：javascript复制splitAssays(accmae, c("01", "10"))
tums <- TCGAsampleSelect(colnames(accmae), "01")
accmae[, tums, ]
只需要 exportClass 函数即可，把前面的 MultiAssayExperiment对象 里面的数据写出来。
代码语言：javascript复制exportClass(accmae, dir = './', fmt = "csv", ext = ".csv")
Writing about 7 files to disk...
[1] ".//accmae_META_0.csv"
[2] ".//accmae_META_1.csv"
[3] ".//accmae_ACC_CNASNP-20160128.csv"
[4] ".//accmae_ACC_CNVSNP-20160128.csv"
[5] ".//accmae_ACC_Mutation-20160128.csv"
[6] ".//accmae_colData.csv"
[7] ".//accmae_sampleMap.csv"
也可以获取临床信息；
代码语言：javascript复制> ## -----------------------------------------------------------------------------
> head(getClinicalNames("ACC"))
[1] "years_to_birth" "vital_status" "days_to_death"
[4] "days_to_last_followup" "tumor_tissue_site" "pathologic_stage"
> colData(accmae)[, getClinicalNames("ACC")][1:5, 1:3]
DataFrame with 5 rows and 3 columns
years_to_birth vital_status days_to_death
<integer> <integer> <integer>
TCGA-OR-A5J1 58 1 1355
TCGA-OR-A5J2 44 1 1677
TCGA-OR-A5J3 23 0 NA
TCGA-OR-A5J4 23 1 423
TCGA-OR-A5J5 30 1 365
> ## -----------------------------------------------------------------------------
一个实战：提取TCGA数据库的BRCA数据集的TNBC亚型的表达量矩阵
前面我们提到过，如果是下载甲基化信号值矩阵或者表达量矩阵，会耗时很长。如果是去ucsc的xena浏览器下载，是一个130M左右的压缩包文件。
代码语言：javascript复制brca <- curatedTCGAData(diseaseCode = "BRCA",
version = "2.0.1",
assays = "RNASeq2GeneNorm", FALSE)
head(getSubtypeMap(brca ))
head(getClinicalNames("BRCA"))
信息如下：
代码语言：javascript复制> head(getSubtypeMap(brca ))
BRCA_annotations BRCA_subtype
1 Patient_ID patientID
2 mrna_subtypes PAM50 mRNA
3 mrna_subtypes SigClust Unsupervised mRNA
4 mrna_subtypes SigClust Intrinsic mRNA
5 microrna_subtypes miRNA Clusters
6 methylation_subtypes methylation Clusters
> head(getClinicalNames("BRCA"))
[1] "years_to_birth" "vital_status" "days_to_death"
[4] "days_to_last_followup" "tumor_tissue_site" "pathologic_stage"
简单的搜索一下 ER.Status ， PR.Status ，HER2.Final.Status
代码语言：javascript复制 phe=as.data.frame(colData( brca ))
# Gender Age.at.Initial.Pathologic.Diagnosis ER.Status PR.Status HER2.Final.Status
> head( phe[,colnames(phe)[grep('Status', colnames(phe),ignore.case = T)][56:59] ])
ER.Status PR.Status HER2.Final.Status Vital.Status
TCGA-A1-A0SB Positive Negative Negative LIVING
TCGA-A1-A0SD Positive Positive Negative LIVING
TCGA-A1-A0SE Positive Positive Negative LIVING
TCGA-A1-A0SF Positive Positive Negative LIVING
TCGA-A1-A0SG Positive Positive Negative LIVING
TCGA-A1-A0SH Negative Positive Negative LIVING
这样的话，就可以根据 ER.Status ， PR.Status ，HER2.Final.Status 这3个信息来判断TNBC样品，进而去筛选前面的表达量矩阵啦。
整体评价：用起来还行，但也不是非他不可！
一个分类模型
如果是分类模型构建，比如我们的模型希望是可以根据转录组测序的表达量来区分肿瘤样品和正常组织，获取数据的代码如下：
代码语言：javascript复制library(curatedTCGAData)
library(TCGAutils)
params <- list(seed = 29221)
brca <- curatedTCGAData(diseaseCode = "BRCA",
version = "2.0.1",
assays = "RNASeq2GeneNorm", FALSE)
## ----data.show, warning=FALSE, error=FALSE------------------------------------
brca <- TCGAutils::splitAssays(brca, c('01','11'))
xdata.raw <- t(cbind(assay(brca[[1]]), assay(brca[[2]])))
library(dplyr)
# Get matches between survival and assay data
class.v <- TCGAbiospec(rownames(xdata.raw))$sample_definition %>% factor
names(class.v) <- rownames(xdata.raw)
# keep features with standard deviation > 0
xdata.raw <- xdata.raw %>%
{ (apply(., 2, sd) != 0) } %>%
{ xdata.raw[, .] } %>%
scale
set.seed(params$seed)
small.subset <- c('CD5', 'CSF2RB', 'HSF1', 'IRGC', 'LRRC37A6P', 'NEUROG2',
'NLRC4', 'PDE11A', 'PIK3CB', 'QARS', 'RPGRIP1L', 'SDC1',
'TMEM31', 'YME1L1', 'ZBTB11',
sample(colnames(xdata.raw), 100))
xdata <- xdata.raw[, small.subset[small.subset %in% colnames(xdata.raw)]]
ydata <- class.v
save(xdata,ydata,file = 'brca_for_lasso.rdata')
上面的代码需要在线下载。
可以看到 肿瘤样品大概是正常样品的10倍数量。
代码语言：javascript复制
> table(ydata)
ydata
Primary Solid Tumor Solid Tissue Normal
1093 112
有了上面的乳腺癌患者的正常组织以及肿瘤样品的表达量矩阵数据，接下来就可以使用函数cv.glmHub来处理xdata和ydata即可进行lasso回归。
如果ydata存的是肿瘤或者正常样品这样的属性，我们使用cv.glmHub来做分类模型构建，如果ydata存的是肿瘤生存信息，那么我们使用cv.glmHub来做生存分析的模型构建。
代码语言：javascript复制library(dplyr)
library(ggplot2)
library(survival)
library(loose.rock)
library(futile.logger)
library(curatedTCGAData)
library(TCGAutils)
library(glmSparseNet)
因为前面ydata存的是肿瘤或者正常样品这样的属性，是二分类，使用函数cv.glmHub 一句话代码即可完成分析：
代码语言：javascript复制load('brca_for_lasso.rdata')
table(ydata)
fitted <- cv.glmHub(xdata, ydata,
family = 'binomial',
network = 'correlation',
nlambda = 1000,
network.options = networkOptions(cutoff = .6,
min.degree = .2))
plot(fitted)
resp <- predict(fitted, s = 'lambda.min', newx = xdata, type = 'class')
table(resp,ydata)
关键的参数是 family = 'binomial'，这个参数告诉我们函数这个时候使用的分类模型。然后可以看到，这个分类器效果还不错：
代码语言：javascript复制> table(resp,ydata)
ydata
resp Primary Solid Tumor Solid Tissue Normal
Primary Solid Tumor 1089 7
Solid Tissue Normal 4 105
基本上没有分错的样品，当然了，还可以去提取构建好的分类器的基因，以及绘制ROC曲线等等，是另外的知识点了，这里仅仅是分享代码：
代码语言：javascript复制coefs.v <- coef(fitted, s = 'lambda.min')[,1] %>% { .[. != 0]}
coefs.v %>% {
data.frame(ensembl.id = names(.),
gene.name = geneNames(names(.))$external_gene_name,
coefficient = .,
stringsAsFactors = FALSE)
} %>%
arrange(gene.name) %>%
knitr::kable()
flog.info('Misclassified (%d)', sum(ydata != resp))
flog.info(' * False primary solid tumour: %d',
sum(resp != ydata & resp == 'Primary Solid Tumor'))
flog.info(' * False normal : %d',
sum(resp != ydata & resp == 'Solid Tissue Normal'))
## ----predict, echo=FALSE, warning=FALSE---------------------------------------
response <- predict(fitted, s = 'lambda.min', newx = xdata, type = 'response')
qplot(response, bins = 100)
## ----roc, echo=FALSE----------------------------------------------------------
roc_obj <- pROC::roc(ydata, as.vector(response))
data.frame(TPR = roc_obj$sensitivities, FPR = 1 - roc_obj$specificities) %>%
ggplot() geom_line(aes(FPR,TPR), color = 2, size = 1, alpha = 0.7)
labs(title= sprintf("ROC curve (AUC = %f)", pROC::auc(roc_obj)),
x = "False Positive Rate (1-Specificity)",
y = "True Positive Rate (Sensitivity)")
可以看到这个分类器的效果超级无敌的好：
分类器的效果超级无敌的好
上面的代码都是复制粘贴就可以运行的。
其实判断肿瘤样品和正常样品是比较容易的，因为肿瘤的特征非常明显，比如引用达到五位数的经典综述《hallmarks of cancer：the next generation》 里面的：
- sustaining proliferative signaling：维持增殖信号
- evading growth suppressors：逃避生长抑制
- resisting cell death：抵抗细胞死亡
- replicative immortality：复制无限
- Angiogenesis：血管生成
- invasion and metastasis：侵袭和转移
- genome instability：基因组不稳定
- Inflammation：炎症
- reprogramming of energy metabolism：能量代谢重编程
- evading immune destruction：逃避免疫破坏
- tumor microenvironment：肿瘤微环境
