介绍
TCGA，癌症基因组图谱，汇集许多肿瘤样本的多组学数据。
我们将讨论如何使用RTCGAToolbox 包获取的数据。说明该系统的大部分工作在于从存档中下载大量资源。工具箱包的插图描述了用于分析的高级实用程序; 在这里，我们将专注于数据协调和手动分析。
以下是直肠腺瘤中三种数据类型下载工作的说明：
library(ph525x)
firehose()
综合
我们使用了命令
library(RTCGAToolbox)
readData = getFirehoseData (dataset="READ", runDate="20150402",forceDownload = TRUE,
Clinic=TRUE, Mutation=TRUE, Methylation=TRUE, RNASeq2GeneNorm=TRUE)
这需要大约10分钟才能获得并节省良好的无线连接。对象的show方法打印
readData
## READ FirehoseData objectStandard run date: 20150402
## Analysis running date: 20160128
## Available data types:
## clinical: A data frame of phenotype data, dim: 171 x 22
## RNASeq2GeneNorm: A matrix of count or normalized data, dim: 20501 x 105
## Methylation: A list of FirehoseMethylationArray object(s), length: 2
## GISTIC: A FirehoseGISTIC for copy number data
## Mutation: A data.frame, dim: 22075 x 39
## To export data, use the 'getData' function.
并隐藏两种甲基化测定的维度。这些可以通过
> lapply(readData@Methylation, function(x) dim(x@DataMatrix))
[[1]]
[1] 27578 76
[[2]]
[1] 485577 109
临床数据的视图
对数据集的完整理解需要注意组装肿瘤“群组”的方法，包括建立疾病进展或死亡事件时间的共同时间来源，以及肿瘤取样和测定程序的细节。出于本课程的目的，我们假设我们可以有意义地组合我们检索到的所有数据。
选择临床参数
clin = getData(readData, "clinical")
names(clin)
## [1] "Composite Element REF"
## [2] "years_to_birth"
## [3] "vital_status"
## [4] "days_to_death"
## [5] "days_to_last_followup"
## [6] "primary_site_of_disease"
## [7] "neoplasm_diseasestage"
## [8] "pathology_T_stage"
## [9] "pathology_N_stage"
## [10] "pathology_M_stage"
## [11] "dcc_upload_date"
## [12] "gender"
## [13] "date_of_initial_pathologic_diagnosis"
## [14] "days_to_last_known_alive"
## [15] "radiation_therapy"
## [16] "histological_type"
## [17] "radiations_radiation_regimenindication"
## [18] "completeness_of_resection"
## [19] "number_of_lymph_nodes"
## [20] "race"
## [21] "ethnicity"
## [22] "batch_number"
疾病的严重程度将在病理学变量中指出。T分期是指肿瘤的大小和侵袭性，N分期是指各种淋巴结中存在癌细胞。
with(clin, table(pathology_T_stage, pathology_N_stage))
## pathology_N_stage
## pathology_T_stage n0 n1 n1a n1b n1c n2 n2a n2b nx
## t1 8 1 0 0 0 0 0 0 0
## t2 28 3 0 1 0 0 0 0 0
## t3 49 30 3 4 1 22 2 2 2
## t4 2 1 0 0 0 2 0 0 0
## t4a 1 1 0 0 0 1 0 5 0
## t4b 0 0 0 0 0 0 1 0 0
我们看到两种分期措施都存在差异。我们将减少T分期以避免小班级。
clin$t_stage = factor(substr(clin$pathology_T_stage,1,2))
table(clin$t_stage)
##
## t1 t2 t3 t4
## 9 32 115 14
定义生存时间
我们猜测，重要状态变量对应于最后一次跟进时的状态，并且最后一次跟进的日期是从诊断过程中的共同起源记录的。患者出现，肿瘤分级并分期，随后的日历开始。
以下Kaplan-Meier显示是粗略的健全性检查，显示肿瘤阶段1没有观察到事件，阶段2和阶3按照我们预期的前1000天左右排序，然后曲线交叉; 数据很稀疏。
library(survival)
ev = 1*(clin$vital == 1)
fut = as.numeric(clin$days_to_last_followup)
su = Surv(fut, ev)
plot(survfit(su~t_stage, data=clin), lwd=2, lty=1:4, xlim=c(0,2000))
ntab = table(clin$t_stage)
ns = paste("[n=", ntab, "]", sep="")
legend(100, .4, lty=1:4, lwd=2, legend=paste(levels(clin$t_stage), ns))
介绍突变数据
mut = getData(readData, "Mutation")
dim(mut)
## [1] 22075 39
table(mut$Variant_Classification)
##
## 3'UTR 5'Flank 5'UTR
## 11 3 14
## De_novo_Start_InFrame De_novo_Start_OutOfFrame Frame_Shift_Del
## 4 31 176
## Frame_Shift_Ins IGR In_Frame_Del
## 161 24 30
## In_Frame_Ins Intron Missense_Mutation
## 8 132 14767
## Nonsense_Mutation Nonstop_Mutation Read-through
## 1964 6 10
## RNA Silent Splice_Site
## 16 4675 43
让我们通过记录的错义或无义突变的数量来命令基因。
gt = table(mut$Hugo, mut$Variant_Classification)
mn = apply(gt[,12:13], 1, sum)
omn = order(mn, decreasing=TRUE)
gt[omn[1:20], c(12:13,17,18)]
##
## Missense_Mutation Nonsense_Mutation Silent Splice_Site
## TTN 79 8 20 0
## APC 9 65 0 1
## TP53 33 8 1 1
## KRAS 37 1 0 0
## MUC16 31 2 12 0
## SYNE1 21 2 6 0
## DNAH10 19 1 2 0
## DNAH5 18 2 2 0
## LRP1B 19 1 4 0
## ABCA13 17 0 1 0
## NEB 14 3 2 0
## ZFHX4 16 1 5 0
## AHNAK2 15 1 1 0
## FAT4 16 0 5 0
## HMCN1 14 2 1 0
## HYDIN 13 3 5 0
## PCLO 13 3 4 0
## PKHD1 16 0 0 0
## SACS 15 1 4 0
## DNAH11 12 2 4 0
KRAS和TP53在此列表中的事实给出了另一个粗略的理智检查。
多组学101：我们可以结合突变和临床数据吗？
这并不简单，因为不共享样本标识符。
clin[1:4,1:3]
## Composite Element REF years_to_birth vital_status
## tcga.af.2687 value 57 0
## tcga.af.2689 value 41 1
## tcga.af.2690 value 76 1
## tcga.af.2691 value 48 0
mut[1:4,c(1,16)]
## Hugo_Symbol Tumor_Sample_Barcode
## 1 OR5B12 TCGA-AG-3586-01A-02W-0831-10
## 2 KRT17 TCGA-AG-3586-01A-02W-0831-10
## 3 SLC10A4 TCGA-AG-3586-01A-02W-0831-10
## 4 C7orf58 TCGA-AG-3586-01A-02W-0831-10
我们猜测以下转换会为突变数据生成适当的标识符。
mid = tolower(substr(mut[,16],1,12))
mid = gsub("-", ".", mid)
mean(mid %in% rownames(clin))
## [1] 1
mut$sampid = mid
让我们简单总结每个人的总突变负担。
nmut = sapply(split(mut$sampid, mut$sampid),length)
nmut[1:4]
## tcga.af.2689 tcga.af.2691 tcga.af.2692 tcga.af.3400
## 64 73 45 33
length(nmut)
## [1] 69
dim(clin)
## [1] 171 23
我们看到并非所有具有临床数据的个体都进行了突变研究。我们将临床数据分组，并通过肿瘤阶段可视化突变计数的分布。
clinwmut = clin[names(nmut),]
clinwmut$nmut = nmut
with(clinwmut, boxplot(split(nmut, t_stage), log="y"))
表达数据
数据没有附带实验级元数据，但我们知道这是通过RSEM进行转录本丰度估计的illumina hiseq RNA测序。如何将数据转换为基因水平需要进行调查。
rnaseq = getData(readData, "RNASeq2GeneNorm")
rnaseq[1:4,1:4]
## TCGA-AF-2687-01A-02R-1736-07 TCGA-AF-2689-11A-01R-A32Z-07
## A1BG 20.1873 43.4263
## A1CF 51.0856 313.3531
## A2BP1 0.4257 18.9911
## A2LD1 90.2639 92.2611
## TCGA-AF-2690-01A-02R-1736-07 TCGA-AF-2691-11A-01R-A32Z-07
## A1BG 56.4619 35.9451
## A1CF 24.9913 218.1571
## A2BP1 0.9256 22.6758
## A2LD1 164.3365 113.6528
我们再次必须转换样本标识符字符串。
rid = tolower(substr(colnames(rnaseq),1,12))
rid = gsub("-", ".", rid)
mean(rid %in% rownames(clin))
## [1] 1
colnames(rnaseq) = rid
遗憾的是，突变和表达数据之间没有太多重叠。
intersect(rid,mid)
## [1] "tcga.af.2689" "tcga.af.2691" "tcga.af.2692" "tcga.af.3400"
## [5] "tcga.ag.3902"
我们注意到一些重复的转换标识符; 不清楚要保留哪个副本，所以我们放下第二个。
which(duplicated(colnames(rnaseq)))
## [1] 11 21 23 25 27 35 104
pairs(log2(rnaseq[101:200,c(10:11,20,21,22,23,50,55)]))
rnaseq = rnaseq[,-which(duplicated(colnames(rnaseq)))]
让我们创建一个ExpressionSet：
library(Biobase)
readES = ExpressionSet(log2(rnaseq+1))
pData(readES) = clin[sampleNames(readES),]
readES
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 20501 features, 98 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: tcga.af.2687 tcga.af.2689 ... tcga.g5.6641 (98
## total)
## varLabels: Composite Element REF years_to_birth ... t_stage (23
## total)
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:
这里有一个人缺失肿瘤阶段，我们将简单地消除这个人。
readES = readES[,-97]
将肿瘤阶段与基因表达变异联系起来
我们将使用非常粗略的分类方法，并将在练习中探索替代方案。我们将使用适度F检验来检验肿瘤分期中常见平均表达的零假设。
library(limma)
##
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
##
## plotMA
mm = model.matrix(~t_stage, data=pData(readES))
f1 = lmFit(readES, mm)
ef1 = eBayes(f1)
## Warning: Zero sample variances detected, have been offset away from zero
topTable(ef1, 2:4)
## t_staget2 t_staget3 t_staget4 AveExpr F
## LOC100128977 -0.2196981 -0.2196981 -0.2196981 0.01132465 16.060640
## CELA2B -0.4924981 -0.4856924 -0.3966889 0.04003472 14.236630
## GDEP -0.9003575 -0.8376533 -0.8484637 0.09571759 11.852702
## FOLR4 -0.6721456 -0.5992014 -0.5847564 0.09479200 10.805454
## C18orf26 -0.3168032 -0.3042073 -0.3168032 0.02516016 10.417589
## GFRA4 -0.2622521 -0.3289370 -0.3481737 0.04383371 10.269518
## TAS2R41 -0.2117273 -0.2051461 -0.2117273 0.01552745 10.213776
## SEMA5B 1.7503244 1.7161443 2.4269386 4.52996758 9.131087
## CTRB2 -1.5604625 -1.5221324 -1.3959781 0.15557817 8.983471
## TXNDC17 -0.7046034 -0.9283568 -1.4545590 9.81107126 8.216157
## P.Value adj.P.Val
## LOC100128977 1.660698e-08 0.0003404596
## CELA2B 1.010504e-07 0.0010358169
## GDEP 1.187769e-06 0.0081168181
## FOLR4 3.647543e-06 0.0186945693
## C18orf26 5.561914e-06 0.0203603949
## GFRA4 6.539887e-06 0.0203603949
## TAS2R41 6.951991e-06 0.0203603949
## SEMA5B 2.311338e-05 0.0592309205
## CTRB2 2.728557e-05 0.0621535018
## TXNDC17 6.519177e-05 0.1336496376
boxplot(split(exprs(readES)["LOC100128977",], readES$t_stage))
介绍450k甲基化数据
450k数据具有与表达数据相同的复杂条件 - 需要转换，重复和仅与临床数据部分匹配的标识符。
me450k = getData(readData, "Methylation", 2)
fanno = me450k[,1:3]
me450k = data.matrix(me450k[,-c(1:3)])
med = tolower(substr(colnames(me450k),1,12))
med = gsub("-", ".", med)
mean(med %in% rownames(clin))
## [1] 1
sum(duplicated(med))
## [1] 8
todrop = which(duplicated(med))
me450k = me450k[,-todrop]
med = med[-todrop]
colnames(me450k) = med
ok = intersect(rownames(clin), colnames(me450k))
me450kES = ExpressionSet(me450k[,ok])
pData(me450kES) = clin[ok,]
fData(me450kES) = fanno
me450kES = me450kES[,-which(is.na(me450kES$t_stage))]
与基因g相关的450k DNA甲基化指标是否与g的表达变化相关？我们需要同步两个数据集。
ok = intersect(sampleNames(me450kES), sampleNames(readES))
meMatch = me450kES[,ok]
esMatch = readES[,ok]
鉴于这些定义，我们可以编写一个辅助函数
corv = function (sym, mpick = 3)
{
mind = which(fData(meMatch)[, 1] == sym)
if (length(mind) > mpick)
mind = mind[1:mpick]
eind = which(featureNames(esMatch) == sym)
dat = cbind(t(exprs(meMatch)[mind, , drop = FALSE]), t(exprs(esMatch)[eind,
, drop = FALSE]), t_stage = jitter(as.numeric(esMatch$t_stage)))
bad = apply(dat, 2, function(x) all(is.na(x)))
if (any(bad))
dat = dat[, -which(bad)]
pairs(dat)
}
corv("ZNF300") # learned about it from firebrowse.org
结论
TCGA是支持多元分析的基础设施开发的明显候选者。
该
MultiAssayExperiment和curatedTCGAData软件包是解决此任务的最新贡献。我们已经看到了一些挑战，即使使用像RTCGAToolbox这样的精心开发的工具来获取数据也是如此：我们必须警惕不匹配的样本标识符标签，缺少数据，样本来源和测定行为的记录不足等等。人类的努力总是需要的; 数据质量标准必须超越数值准确性并解决透明度和可用性问题。curatedTCGAData包使用手动策划的表示，因此本节中演示的协调任务不是必需的。但总的来说，我们需要知道如何协调，所以这些任务在本课程中保留。
通过适当改变本文档中的代码片段，您可以访问其他模态的多元数据（包括microRNA，拷贝数变异和蛋白质组学）。当您发现解释这些措施以创建生物学洞察力的新方法时，请通过向Bioconductor添加包或工作流程将它们传达给全世界。
原创： 小伍说事 医知圈
网友评论