TCGA数据下载就易用性来说，RTCGA包应该更好用，且由于是已经下载好的数据，使用比较稳定。但是也由于是下载好的数据，不能保证数据都是全新的。TCGAbiolinks包是实时调用GDC的API，所以可以获取最新的数据。
数据下载三部曲
数据下载三部曲GDCquery、GDCdownload、GDCprepare。
GDCquery用于查询GDC数据库，里面获取所有需要下载的TCGA数据的各项记录。
GDCdownload根据GDCquery的检索结果进行文件下载。
下载完成后，GDCprepare同样根据GDCquery的文件结果可以将下载数据规整为summarizedExperiment对象或者是返回一个data.frame。
目前有两大类TCGA数据可供下载，一个是Legacy，主要是一些使用 GRCh37 (hg19) 和GRCh36 (hg18)的数据，另一个是harmonized数据，统一使用GRCh38 (hg38)作为参考序列。
There are two available sources to download GDC data using TCGAbiolinks:
- GDC Legacy Archive : provides access to an unmodified copy of data that was previously stored in CGHub and in the TCGA Data Portal hosted by the TCGA Data Coordinating Center (DCC), in which uses as references GRCh37 (hg19) and GRCh36 (hg18).
- GDC harmonized database: data available was harmonized against GRCh38 (hg38) using GDC Bioinformatics Pipelines which provides methods to the standardization of biospecimen and clinical data.
这两种的GDCquery的参数会有少许不同，这里主要以harmonized数据为主，下载TCGA-READ和TCGA-COAD项目的RNA-seq数据。
代码语言：javascript复制library(TCGAbiolinks)
library(tidyverse)
# 获取全部的TCGA项目信息
AllProj <- getGDCprojects()
# 每个项目的摘要
proj_READ <- TCGAbiolinks:::getProjectSummary("TCGA-READ")
proj_COAD <- TCGAbiolinks:::getProjectSummary("TCGA-COAD")
根据vignette的内容（http://www.bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/query.html）中的Harmonized data options (legacy = FALSE)，转录组数据的Data.category是Transcriptome Profiling，Data.type是Gene Expression Quantification，WorkflowType有四种：STAR - Counts、HTSeq - FPKM-UQ、HTSeq - FPKM和HTSeq - Counts。
这里选择下载HTSeq - Counts，也就是RawCounts，不使用FPKM Normalization数据，后面的Normalization使用DESeq2来做。
代码语言：javascript复制# query
query_READ <- GDCquery(
project = "TCGA-READ",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts", # "STAR - Counts"
)
query_COAD <- GDCquery(
project = "TCGA-COAD",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts" # "STAR - Counts"
)
下载使用GDCdownload，由于TCGA的下载不是特别稳定，所以可以使用files.per.chunk定为一个值，几个文件打包为一个压缩文件来下载。另外如果这里的下载如果不成功的话，需要重复运行几次，直至完全下载成功。
代码语言：javascript复制# download
GDCdownload(
query = query_READ,
method = "api",
files.per.chunk = 5
)
GDCdownload(
query = query_COAD,
method = "api",
files.per.chunk = 5
)
GDCprepare后即可拿到TCGA的数据，默认是会返回一个summarizedExperiment对象。这个过程中，GDCprepare还会将生存数据自动合并到summarizedExperiment对象的colData中。
代码语言：javascript复制summarizedExperiment对象和ExpressionSet等对象类型类似，核心组件就是三大件：表达量、列注释和行注释。 表达量：一个表达量矩阵，行是基因或者相关特征，列是样本或相关特征; 列注释：样本相关的注释，比如病人信息、生存数据等等； 行注释：基因相关的注释，比如基因名称、长度、位置、ID等等。
# prepare
dat_READ <- GDCprepare(
query = query_READ,
save = TRUE,
save.filename = "data_READ.rda",
remove.files.prepared = FALSE
)
dat_COAD <- GDCprepare(
query = query_COAD,
save = TRUE,
save.filename = "data_COAD.rda",
remove.files.prepared = FALSE
)
生存分析
TCGAbiolinks也有自定义的分析函数，这里使用survival和survminer包手动分析生成数据。
数据合并并使用DEseq2的Normalization方法。
代码语言：javascript复制library(DESeq2)
# 合并数据
dat_couts <-
list(dat_counts_READ,
dat_counts_COAD) %>%
map(~rownames_to_column(.x, "Gene")) %>%
do.call(left_join, .) %>%
column_to_rownames("Gene")
dds <- DESeq2::DESeqDataSetFromMatrix(
dat_couts,
colData = data.frame(row.names = colnames(dat_couts)),
design = ~1
)
dds <- DESeq2::estimateSizeFactors(dds)
dat_couts_norm <- counts(dds, normalized = TRUE)
RNA-seq的Normalization方法有很多种，但是最不建议的就是RPKM或者FPKM。
Normalization需要控制的三个不均衡因素是文库大小、基因长度及文库组成：
- 文库大小：比如样本A是样本B的测序量的两倍，那么在同等表达水平下，样本A的基因的Counts值就是B的两倍；
- 基因长度：如果需要进行Gene间的比较，那么需要对基因长度做Normalization，否则的话，基因越长，获得的基因Counts也就越多；
- 文库组成：比如样本A是敲除样本，而样本B是过表达样本，那么这个时候我们就不能假定两个样本的转录组总表达量是一致的，样本A很有可能是要比样本B拥有更多的Counts数量的，这个情况下就会对TPM、CPM及FPKM/RPKM等数据提出挑战。
DESeq2的Normalization方法，已经有很多资料了，这里只说它的效果就是可以校正文库大小和文库组成，也就是说可以进行样本间比较，无法进行基因间比较，大多数情况下，我们都是不需要基因间比较的。
对生存数据做数据清洗，关键信息是vital_status, days_to_death, days_to_last_follow_up，根据vital_status选择生存数据，但是有些病例的生存数据比较奇怪：正常情况下如果vital_status是Dead，那么生存时间是days_to_death，否则生存时间是days_to_last_follow_up，但是也有days_to_death和days_to_last_follow_up同时为NA的情况，这种情况下直接丢弃这些样本。
代码语言：javascript复制dat_surv_READ <-
data_READ@colData %>%
as.data.frame() %>%
dplyr::select(
barcode,
sample_type,
vital_status,
days_to_death,
days_to_last_follow_up
) %>%
dplyr::filter(sample_type == "Primary Tumor") %>%
mutate(status = ifelse(vital_status == 'Dead', "1", "0")) %>%
mutate(
OS = case_when(
status == 1 ~ days_to_death,
status == 0 ~ days_to_last_follow_up,
TRUE ~ NA_integer_
)
) %>%
dplyr::filter(!is.na(OS))
dat_surv_COAD <-
data_COAD@colData %>%
as.data.frame() %>%
dplyr::select(
barcode,
sample_type,
vital_status,
days_to_death,
days_to_last_follow_up
) %>%
dplyr::filter(sample_type == "Primary Tumor") %>%
mutate(status = ifelse(vital_status == 'Dead', "1", "0")) %>%
mutate(
OS = case_when(
status == 1 ~ days_to_death,
status == 0 ~ days_to_last_follow_up,
TRUE ~ NA_integer_
)
) %>%
dplyr::filter(!is.na(OS))
# 合并READ COAD
dat_surv <- rbind(dat_surv_READ, dat_surv_COAD)
# barcode的前15位是病人ID，根据barcode去重
dat_surv <- dat_surv %>% mutate(Patient_ID = str_sub(barcode, 15))
dat_surv_uniq <- dat_surv %>% distinct(Patient_ID, .keep_all = TRUE)
随机选择100个基因的Normalization Data做后续分析，合并Normalization Data和生成分析数据。
代码语言：javascript复制### 随机选100个基因
filter_gene <- sample(rownames(dat_couts_norm), 100)
filter_dat_norm <-
dat_couts_norm[filter_gene, ] %>% t %>%
as.data.frame %>%
rownames_to_column("barcode")
filter_dat_norm_full <- left_join(dat_surv_uniq, filter_dat_norm, by = "barcode")
使用survival进行生存分析，使用survminer进行可视化。
生存分析时根据基因的中位数将其分为High和Low，使用log-rank检验显著性，也可以使用cox回归。log-rank和cox回归的区别在于是cox是半参数检验，需要对数据有一些先验假设，另外cox回归并不不局限于拟合数据是分类变量，也可以是连续变量。
代码语言：javascript复制library(survival)
library(survminer)
res_surv_full <-
map(filter_gene,
function(x){
# 由于是随机选的基因，不少基因几乎没有表达量
# 这里做一下判断，如果一个基因的均值低于5，就不做生存分析了
if(mean(filter_dat_norm_full[[x]]) < 5){
return(NULL)
}
# 根据基因表达的中位数标记为High和Low
dat <- filter_dat_norm_full %>%
mutate('{x}' := ifelse(.data[[x]] > median(.data[[x]]), "High", "Low"))
# 生存信息status保证是数值类型
dat$status <- as.numeric(dat$status)
# 手动构造生存分析的拟合公式
pattern <- str_glue("Surv(OS, status) ~ {x}")
formu <- pattern %>% as.formula()
# 拟合
fit <- survfit(formu, data=dat)
fit$call$formula <- formu
# 差异分析
diff <- survdiff(formu, data=dat)
diff$call$formula <- formu
p.value <- 1 - pchisq(diff$chisq, length(diff$n) -1)
# 返回结果
list(fit = fit,
diff = diff,
p.value = p.value,
dat = dat)
})
# 未做生存分析的基因直接筛除
res_surv_full <- res_surv_full %>% .[!map_lgl(., is.null)]
length(res_surv_full)
# [1] 42
# 此时的res_surv_full的每个元素都是一个四元素列表
res_surv_full[[1]] %>% str(max.level = 1)
#List of 4
# $ fit :List of 17
# ..- attr(*, "class")= chr "survfit"
# $ diff :List of 6
# ..- attr(*, "class")= chr "survdiff"
# $ p.value: num 0.294
# $ dat :'data.frame': 94 obs. of 108 variables:
可视化，ggsurvplot函数返回的是一个ggsurvplot对象的列表，选1个看一下结果，感兴趣的可以做后续的拼图、定制等后续美化。
代码语言：javascript复制res_surv_plots <-
map(res_surv_full,
function(x){
ggsurvplot(x$fit,
data = x$dat,
pval = TRUE,
risk.table = TRUE)
})
# 查看第一个结果
res_surv_plots[[1]]
ggsurvplot对象其实就一个基于列表的S3对象，里面的plot就是实际的ggplto2对象，如果有添加risk.table的话，那么里面的table元素就是实际的ggplto2对象。 可以自己提取元素plot和table，然后使用patchwork或者cowplot合并，则可以将ggsurvplot转为ggplot2对象，然后就可以自由的拼合多个生成图形了。
