Workshop-TCGA data analysis

rpubs.com/tiagochst/TCGA...

1 工作描述

1.1 目标

本次研讨会将:

1. 在 NCI 的基因组数据共享(GDC)中可用的 TCGA(癌症基因组图谱)数据。
2. 演示使用 R/Bioconductor 包 tcgabiollinks 访问和导入 TCGA 数据。
3. 提供使用 R/Bioconductor 软件包 maftools 可视化拷贝数改变和突变数据的指导。

1.2 可用的链接

• the NCI's Genomic Data Commons (GDC): https://gdc.cancer.gov/.
• TCGAbiolinks package: http://bioconductor.org/packages/TCGAbiolinks/
• maftools package: http://bioconductor.org/packages/maftools/

1.3 参考文献

• TCGAbiolinks:

  • Colaprico, Antonio, et al. “TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data.” Nucleic acids research 44.8 (2015): e71-e71. Silva, Tiago C., et al. “TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages.” F1000Research 5 (2016). (https://f1000research.com/articles/5-1542/v2)
  • Mounir, Mohamed, et al. “New functionalities in the TCGAbiolinks package for the study and integration of cancer data from GDC and GTEx.” PLoS computational biology 15.3 (2019): e1006701. (https://doi.org/10.1371/journal.pcbi.1006701)
  • Silva TC, Colaprico A, Olsen C et al.TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages [version 2; peer review: 1 approved, 2 approved with reservations]. F1000Research 2016, 5:1542 (https://doi.org/10.12688/f1000research.8923.2).

• The NCI’s Genomic Data Commons (GDC):

  • Gao, Galen F., et al. “Before and After: Comparison of Legacy and Harmonized TCGA Genomic Data Commons’ Data.” Cell systems 9.1 (2019): 24-34. (https://doi.org/10.1016/j.cels.2019.06.006)

• Maftools:

  • Mayakonda A, Lin D, Assenov Y, Plass C, Koeffler PH (2018). “Maftools: efficient and comprehensive analysis of somatic variants in cancer.” Genome Research. doi: 10.1101/gr.239244.118.

• Complexheatmap:

  • Gu Z, Eils R, Schlesner M (2016). “Complex heatmaps reveal patterns and correlations in multidimensional genomic data.” Bioinformatics.

1.4 使用 R/Bioconductor 包

• TCGAbiolinks
• MultiAssayExperiment
• maftools
• ComplexHeatmap

1.5 前提条件

• Basic knowledge of R syntax
• Understand the pipe operator (“%>%”) (help material https://r4ds.had.co.nz/pipes.html)
• Understand the SummarizedExperiment data structure (help material http://bioconductor.org/packages/SummarizedExperiment/)

2 引言

2.1 了解 GDC(基因组学数据公共数据库)

在本次研讨会中，我们将访问 NCI 基因组数据共享(GDC)中可用的(癌症基因组图谱)数据(portal.gdc.cancer.gov/)。在我们开始之前，重要的是要知道数据存储在两个不同的数据库中:

​​

• **==The legacy archive==**(https://portal.gdc.cancer.gov/legacy-archive/search/f)，其中包含来自早于 GDC(例如 chub)的存储库的未协调遗留数据。此外，遗留数据没有被 GDC 积极维护、处理或协调(来源:https://docs.gdc.cancer.gov/Data_Portal/Users_Guide/Legacy_Archive/)
• ==Harmonized database==(https://portal.gdc.cancer.gov/)，其中包含使用标准化管道和参考基因组 GRCh38 (hg38)处理的许多癌症项目的数据。这提供了分析多种癌症类型或跨多个项目分析同一癌症类型的优势。您可以在 https://gdc.cancer.gov/about-data/gdc-data-harmonization 和 https://docs.gdc.cancer.gov/Encyclopedia/pages/Harmonized_Data/找到有关支持数据协调的管道的更多信息。

3.2 理解 TCGA 数据

3.2.1 数据访问

GDC 提供了两种访问级别的数据:

​​

• Open: 包括不能单独识别的高层次基因组数据，以及大多数临床和所有生物标本数据元素。
• Controlled: 包括个体可识别的数据，如低水平基因组测序数据、种系变异、SNP6 基因型数据和某些临床数据元素
• 详细教程链接：gdc.cancer.gov/access-da...

3.2.2 TCGA 条码描述

每个 TCGA 样品都有一个唯一的标识符，称为 TCGA 条形码，它包含了每个样品的重要信息。条形码的描述如下所示。(详细内容: https://docs.gdc.cancer.gov/Encyclopedia/pages/TCGA_Barcode/) TCGA 样本命名规则详解 | 快速获取样本类型、组织来源等分组信息¹

​​

​​

样品类型编号

​​

3.2.3 数据结构

为了筛选 GDC 中可用的数据，可以使用项目(TCGA, TARGET 等)，数据类别(转录组分析 ，​DNA 甲基化 ，​临床等)，数据类型(基因表达量化[Gene Expression Quantification]，==Isoform 表达量化==[Isoform Expression Quantification]，**甲基化 β 值**[Methylation Beta Value]，实验策略(miRNA-Seq, RNA-Seq 等)，工作流类型，平台，访问类型等。

​​

就数据粒度而言，一个项目有几个类别的数据，每个类别包含几个数据类型，这些数据类型可能是由不同的工作流、实验策略和平台产生的。这样，如果您选择数据类型“基因表达量化[Gene Expression Quantification]”，数据类别将是转录组分析。

​​

3.3 SummarizedExperiment 数据结构

在我们开始之前，重要的是要知道 R/Bioconductor 环境提供了一个名为 SummarizedExperiment 的数据结构，该数据结构用于处理同一对象中的样本元数据(年龄，性别等)，基因组学数据(即 DNA 甲基化 β 值)和基因组学元数据信息(chr，开始，结束，基因符号)。您可以使用 colData 函数访问样本元数据，使用 assays 访问基因组数据，使用 rowRanges 访问基因组元数据。

​​

3.4 加载需要的 R 包

library(TCGAbiolinks)
library(MultiAssayExperiment)
library(maftools)
library(dplyr)
library(ComplexHeatmap)

如果这些软件包不可用，您可以使用 BiocManager 安装这些软件包。

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!requireNamespace("ComplexHeatmap", quietly = TRUE))
  BiocManager::install("ComplexHeatmap") # from bioconductor

if (!requireNamespace("TCGAbiolinks", quietly = TRUE))
  BiocManager::install("BioinformaticsFMRP/TCGAbiolinks") # from github

4 处理 TCGA 数据

在本次研讨会期间，我们将使用来自统一 GDC 数据库(hg38 版本)的开放获取数据，该数据库可在 https://portal.gdc.cancer.gov/上访问。

通过接下来的部分，我们将:

4.1 临床数据

在 GDC 数据库中，临床数据可以从不同的来源检索:

• 索引临床:使用 XML 文件创建的精细临床数据。
• XML 文件:数据的原始来源
• BCR Biotab:从 XML 文件解析的 tsv 文件

索引临床文件和 XML 文件之间有两个主要区别:

• XML 有更多的信息:辐射、药物信息、随访、生物标本等。因此，索引文件只是 XML 文件的一个子集
• 索引的数据包含带有后续信息的更新数据。例如:如果患者在收集临床数据时还活着，然后在下一次随访中他已经死亡，则索引数据将显示死亡。XML 将有两个字段，第一个字段表示他还活着(在临床部分)，第二个字段表示他已经死亡。

其他有用的临床资料包括:

• 组织切片图像
• 病理报告-幻灯片图像

您可以在 tcgabiollinks 的插图中查看示例 https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/clinical.html.

​​

# Access indexed clinical data
clinical <- GDCquery_clinic("TCGA-COAD")
head(clinical)

​​

临床索引数据与您在 GDC 数据门户中看到的相同。观察:GDC API 提供以天为单位的 age_at_diagnosis。

# Same case as figure above
clinical %>% 
  dplyr::filter(submitter_id == "TCGA-AA-3562") %>% 
  t %>% 
  as.data.frame

​​

您还可以访问从 XML 文件创建的 TSV 文件的 Biotab 文件(仅适用于 tcgabiollinks 2.13.5 以上版本)。

query <- GDCquery(project = "TCGA-ACC", 
                  data.category = "Clinical",
                  data.type = "Clinical Supplement", 
                  data.format = "BCR Biotab")
GDCdownload(query)
clinical.BCRtab.all <- GDCprepare(query)
names(clinical.BCRtab.all)

# [1] "clinical_omf_v4.0_acc"           "clinical_follow_up_v4.0_nte_acc"
## [3] "clinical_radiation_acc"          "clinical_patient_acc"         
## [5] "clinical_nte_acc"                "clinical_follow_up_v4.0_acc"  
## [7] "clinical_drug_acc"

clinical.BCRtab.all$clinical_drug_acc  %>% 
  head  %>% 
  as.data.frame

​​

4.2 RNA 测序数据

RNA-Seq 流水线产生原始计数，FPKM 和 FPKM- uq 定量，并进行了描述 3. 来自 GDC 官网 mRNA 分析流程²

https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/.

​​

​​

以下选项用于使用 tcgabiollinks 搜索 mRNA 结果:

• data.category: “Transcriptome Profiling”
• data.type: “Gene Expression Quantification”
• workflow.type: “HTSeq - Counts”, “HTSeq - FPKM”, “HTSeq - FPKM-UQ”

Here is the example to download the raw counts, which can be used with DESeq2 (http://bioconductor.org/packages/DESeq2/) for differential expression analysis.

query.exp.hg38 <- GDCquery(project = "TCGA-GBM", 
                  data.category = "Transcriptome Profiling", 
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts",
                  barcode =  c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"))
GDCdownload(query.exp.hg38)
raw.counts <- GDCprepare(query = query.exp.hg38, summarizedExperiment = FALSE)
head(raw.count)

​​

下面是下载 FPKM-UQ 计数的示例。

query.exp.hg38 <- GDCquery(project = "TCGA-GBM", 
                  data.category = "Transcriptome Profiling", 
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - FPKM-UQ",
                  barcode =  c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"))
GDCdownload(query.exp.hg38)
fpkm.uq.counts <- GDCprepare(query = query.exp.hg38, summarizedExperiment = FALSE)
head(fpkm.uq.counts)

​​

4.3 突变分析

tcgabiollinks 提供了一些从 GDC 下载突变数据的功能。下载数据有两种方式:

1. 使用 GDCquery_Maf 将下载与 hg38 对齐的 MAF。

2. 您可以使用 tcgabiollinks 的 GDCquery_Maf 函数下载数据。

# GDCquery_Maf download the data from GDC
maf <- GDCquery_Maf("COAD", pipelines = "muse")
maf %>% head %>% as.data.frame

​​

# using maftools for data summary 
maftools.input <- maf %>% read.maf

# Check summary
plotmafSummary(maf = maftools.input, 
               rmOutlier = TRUE, 
               addStat = 'median', 
               dashboard = TRUE)

​​

oncoplot(maf = maftools.input, 
         top = 10, 
         removeNonMutated = TRUE)

​​

# classifies Single Nucleotide Variants into Transitions and Transversions
titv = titv(maf = maftools.input, 
            plot = FALSE, 
            useSyn = TRUE)
plotTiTv(res = titv)

​​

getSampleSummary(maftools.input)

​​

4.4 拷贝数变异数据

拷贝数变异分析管道在 https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/上描述了数字焦点级拷贝数变异(CNV)值，该值是在项目级别上使用 GISTIC2 从肿瘤同质组中使用“屏蔽拷贝数片段”文件生成的。仅保留蛋白质编码基因，其数字 CNV 值通过 0.3 的噪声截断进一步阈值化:

• Genes with focal CNV values smaller than -0.3 are categorized as a “loss” (-1)
• Genes with focal CNV values larger than 0.3 are categorized as a “gain” (+1)
• Genes with focal CNV values between and including -0.3 and 0.3 are categorized as “neutral” (0).

query <- GDCquery(project = "TCGA-GBM",
             data.category = "Copy Number Variation",
             data.type = "Gene Level Copy Number Scores",            
             access = "open")
GDCdownload(query)
scores <- GDCprepare(query)
scores[1:5,1:5]

​​

scores.matrix <- scores %>% 
  dplyr::select(-c(1:3)) %>%  # Removes metadata from the first 3 columns
  as.matrix

rownames(scores.matrix) <- paste0(scores$`Gene Symbol`,"_",scores$Cytoband)

# gain in more than 100 samples
gain.more.than.hundred.samples <- which(rowSums(scores.matrix == 1) > 100)

# loss in more than 100 samples
loss.more.than.hundred.samples <- which(rowSums(scores.matrix == -1) > 100)

lines.selected <- c(gain.more.than.hundred.samples,loss.more.than.hundred.samples)

Heatmap(scores.matrix[lines.selected,],
        show_column_names = FALSE, 
        show_row_names = TRUE,
         row_names_gp = gpar(fontsize = 8),
        col = circlize::colorRamp2(c(-1,0,1), colors = c("red","white","blue")))

​​

4.5 DNA 甲基化数据

处理后的 DNA 甲基化数据测量已知 CpG 位点的甲基化水平为 β 值，从阵列强度(2 级数据)计算为 β = M/(M+U)。

(Zhou, Laird, and Shen 2017)，范围从 0 到 1 完全甲基化。https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Methylation_LO_Pipeline/.

我们将下载两个胶质母细胞瘤(GBM)作为总结实验对象。

query_met.hg38 <- GDCquery(project = "TCGA-GBM", 
                           data.category = "DNA Methylation", 
                           platform = "Illumina Human Methylation 27", 
                           barcode = c("TCGA-02-0116-01A","TCGA-14-3477-01A-01D"))
GDCdownload(query_met.hg38)
data.hg38 <- GDCprepare(query_met.hg38,summarizedExperiment = TRUE)
data.hg38

## class: RangedSummarizedExperiment 
## dim: 27578 2 
## metadata(1): data_release
## assays(1): ''
## rownames(27578): cg00000292 cg00002426 ... cg27662877 cg27665659
## rowData names(7): Composite.Element.REF Gene_Symbol ...
##   CGI_Coordinate Feature_Type
## colnames(2): TCGA-02-0116-01A-01D-0199-05
##   TCGA-14-3477-01A-01D-0915-05
## colData names(113): sample patient ...
##   subtype_Telomere.length.estimate.in.blood.normal..Kb.
##   subtype_Telomere.length.estimate.in.tumor..Kb.

使用 rowrange 访问探针信息

data.hg38 %>% rowRanges %>% as.data.frame

​​

使用 colData 访问示例元数据。

data.hg38 %>% colData %>% as.data.frame

​​

通过 assay 获得 DNA 甲基化水平。

data.hg38 %>% assay %>% head %>% as.data.frame

​​

# plot 10 most variable probes
data.hg38 %>% 
  assay %>% 
  rowVars %>% 
  order(decreasing = TRUE) %>% 
  head(10) -> idx

pal_methylation <- colorRampPalette(c("#000436",
                                      "#021EA9",
                                      "#1632FB",
                                      "#6E34FC",
                                      "#C732D5",
                                      "#FD619D",
                                      "#FF9965",
                                      "#FFD32B",
                                      "#FFFC5A"))(100)

Heatmap(assay(data.hg38)[idx,],
        show_column_names = TRUE, 
        show_row_names = TRUE,
         name = "Methylation Beta-value", 
         row_names_gp = gpar(fontsize = 8),
         column_names_gp = gpar(fontsize = 8),
        col = circlize::colorRamp2(seq(0, 1, by = 1/99), pal_methylation))

​​

4.6 ATAC 数据

rpubs.com/tiagochst/atac...