TCGA

Author

Anqi Dong

Published

May 30, 2026

本文档使用TCGAbiolinks包对 TCGA 数据库中 bulk mRNA 测序数据进行下游分析。TCGAbiolinks包会自动下载好 raw counts、tpm 和 fpkm 方便使用。相比于 UCSC XENA 网站的数个月的延迟，TCGAbiolinks包直接从 TCGA 官网下载最新、最全的数据。

以 TCGA-COAD 数据集作为示例。

library(TCGAbiolinks)
library(SummarizedExperiment)

# query, download, prepare ----
query <- GDCquery(
  project = "TCGA-COAD",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification",
  workflow.type = "STAR - Counts",
  sample.type = c("Primary Tumor")
)

GDCdownload(query, directory = '~/data/TCGA-COAD', files.per.chunk = 6)

se <- GDCprepare(query, save = T, directory = '~/data/TCGA-COAD', 
                   save.filename = 'TCGA-COAD_mRNA.Rdata')

# keep mRNA data
se_mrna <- se[rowData(se)$gene_type == "protein_coding",]

# extract TPM data ----
tpm <- assay(se_mrna, 'tpm_unstrand')
## unstrand = raw counts

se和se_mrna数据类型均为SummarizedExperiment，可以类比 Python 中的 AnnData 进行理解。

tpm是一个行为基因、列为样本的 matrix。

此处根据某基因的表达量高低将样本分为高中低三等分，取高表达组和低表达组使用DESeq2包计算差异基因。直接将 se_mrna 传入 DESeqDataSet() 函数中，该函数会自动将 se_mrna 中的第一个 assay （unstrand）当作 raw counts 用于计算。

loc <- which(rowData(se_mrna)$gene_name == 'GENE_A')
rowData(se_mrna)[loc,]
colData(se_mrna)$group_by_gene <- ntile(tpm[loc, ], 3) |> as.factor()

# DESeq2 ----
## pass SE object to DESeqDataSet() function
ddsSE <- DESeqDataSet(se_mrna, design = ~ group_by_gene)
## renaming the first element in assays to 'counts'
ddsSE

## filter out low-expression genes
ddsSE <- ddsSE[rowSums(counts(ddsSE)) > 10, ]

## calculate DEGs
ddsSE <- DESeq(ddsSE)

# 注意此处contrast参数的顺序，表示分组信息储存在group_by_gene中，计算3相对于1的变化，也就是 > 0 表示在3组上调。
res <- results(ddsSE, contrast = c("group_by_gene", "3", "1"))

head(res)

res_sig <- res |> as.data.frame() |> 
  tibble::rownames_to_column(var = 'gene') |> 
  dplyr::mutate(signif = case_when(
    log2FoldChange > 1 & pvalue < 0.05 ~ 'up',
    log2FoldChange < -1 & pvalue < 0.05 ~ 'down',
    .default = 'ns'
  )) |> as_tibble()
dplyr::count(res_sig, signif)