玩转 TCGA 数据库 - 转录组分析（二）

生信菜鸟团

发布于 2025-05-13 01:22:26

5000

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

转录组数据分析

TCGA 数据库因有各种的癌症转录组数据，所以最基础的就是下载转录组数据，然后对不同样本根据表型做转录组分析。这里我们做一下实战。

转录组表达矩阵下载

1. 加载包

library(tidyverse)
library(TCGAbiolinks)
library(SummarizedExperiment)

2. 查询 project 信息

我们以 PAAD 这个癌症种类为例子。

project <- "TCGA-PAAD" #⭐️ 关注 project
TCGAbiolinks:::getProjectSummary(project)

可以看到有这么多数据类型。

其实我们知道我们想要转录组数据，其实就是Transcriptome Profiling, 完全可以不要这一步，但忘了怎么拼写可以看一下。

$file_count
[1] 12956

$data_categories
   file_count case_count                data_category
1        4100        185  Simple Nucleotide Variation
2        1537        185             Sequencing Reads
3        1032        185                  Biospecimen
4         396        185                     Clinical
5        2700        185        Copy Number Variation
6         732        178      Transcriptome Profiling
7         585        184              DNA Methylation
8         120        120           Proteome Profiling
9         858        173 Somatic Structural Variation
10        896        181         Structural Variation

$case_count
[1] 185

$file_size
[1] 1.326276e+14

3. 查询 & 下载 & 准备数据

这一步使用 TCGAbiolinks 包可以直接流程化完成。如果想要下载这个癌种的转录组数据，就是这个代码，不用改任何信息。想要其他癌症的，在 project这个变量名就可以。

这一步出错，网路原因占据大多数，不过我没出现问题。

query <- GDCquery(project = project,
                   data.category = "Transcriptome Profiling",
                   data.type = "Gene Expression Quantification",
                   workflow.type = "STAR - Counts"
                  )
GDCdownload(query)
data <- GDCprepare(query = query)

我们可以看一下 data这个数据结构，他使用的是 SummarizedExperiment 对象。https://bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html

SummarizedExperiment 对象包含三个可通过 SummarizedExperiment 包访问的主要矩阵：

样本矩阵信息通过 colData(data) 访问：存储样本信息。TCGAbiolinks 会添加来自 TCGA 标志性论文的索引临床数据和亚型信息。类似于表达矩阵的列名，代表样本名，不过对这个样本名进行了多个信息记录。
检测矩阵信息通过 assay(data) 访问：存储分子数据。没有行名和列名的矩阵。
基因矩阵信息（基因信息）通过 rowRanges(data) 访问：存储特征的元数据，包括其基因组范围。类似于表达矩阵的行名，基因 ID 一样，不过对这个 ID 进行了多个信息记录。

说的通俗一点，这个结构相当于一个表达矩阵，不过我分成了三部分储存，一部分是行名（基因信息），一部分是列名（样本信息），一部分是矩阵（counts, TPM等）。

基因信息

我们分别谈每一部分，先看基因信息。

# 可以直接查看
data@rowRanges

# 也可以使用 SummarizedExperiment 对象专有提取命名。这里的rowRanges可以一个函数命令哦。
rowRanges(data)

可以看到这个其实是一个表格，第一列就是基因ID，其余列是对这个基因更详细的标注。

附：我们可以对这些基因进行过滤，比如我就想研究蛋白编码基因，或者你觉得那些非编码基因会影响到我们的统计学检验，那就可以设置条件去过滤不想要的基因。这一步看自己想法，怎样都行（可选）。

data1 <- data[which(mcols(rowRanges(data))$gene_type == "protein_coding"),]

样本信息

接着看样本信息。这个信息就很多了，我就转换为数据框查看了。

sample <- colData(data) %>% as.data.frame()

可以看到有135个样本信息统计呀，这里使用的样本ID是叫做barcode，类似于我们的身份证证号，不同的字段也代表着不同的意思，详细内容自行探索吧

https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables。TCGAbiolinks包自动整合了样本的统计信息和临床信息，所以信息不可谓不丰富。后续样本分组也主要是靠这一步。比如我们想看一下癌症和非癌症之间的转录表达谱差异，那我们首先就是得知道哪些样本是癌症的哪些不是。那我们就看哪个字段是我们想要的，然后做统计识别。

# 先样本分布统计
> sample$tumor_descriptor %>% table()
.
    Metastatic Not Applicable        Primary 
             1              4            178

那就是一共178 个原发，四个正常，一个转移的。因为样本的信息非常多，不知道哪个字段是什么意思的，用人工智能查一下也是很准确的。

那都做的这里了，我们就做一个样本分组文件吧。一列是样本 ID，一列是关注的表型，这里我们拿肿瘤和非肿瘤分组。

这里要特别注意自己是否分组正确了，根据自己的需求进行分组。

design <- data.frame(
  sample = sample$barcode,
  group = ifelse(sample$tumor_descriptor == "Not Applicable", "Normal", "Tumor")) %>% 
  arrange(group, "Normal", "Tumor")
design$group <- factor(design$group, levels = c("Normal", "Tumor"))

附：也可以在 data对象里面过滤掉不想要的样本，比如我只想要原发的肿瘤样本 (可选)：

data2 <- data[, data$tumor_descriptor == "Primary"]

表达矩阵

表达矩阵储存在 data@assays@data@listData 这个对象结构下，并且也有多个表达矩阵

> data@assays@data@listData %>% names()
[1] "unstranded"       "stranded_first"   "stranded_second"  "tpm_unstrand"     "fpkm_unstrand"   
[6] "fpkm_uq_unstrand"

unstranded：表示 RNA-Seq 数据在处理时没有考虑链的方向性。
stranded_first：表示 RNA-Seq 数据是链特异性的，并且第一条链（通常是正义链）被考虑。
stranded_second：表示 RNA-Seq 数据是链特异性的，并且第二条链（通常是反义链）被考虑。
tpm_unstrand：表示未区分链的 TPM（Transcripts Per Million）值。这是另一种标准化基因表达的方法，考虑了转录本长度和测序深度。
fpkm_unstrand：表示未区分链的 FPKM（Fragments Per Kilobase of transcript per Million mapped reads）值，用于标准化基因表达。
fpkm_uq_unstrand：表示未区分链的 Upper Quartile-normalized FPKM 值，这是一种通过上四分位数标准化的 FPKM，用于减少样本间的技术变异。

这一步可以使用assay函数去提取。最好不要直接data@assays@data@listData[["unstranded"]]，虽然都是一样的矩阵，但assay函数会自动的给你加上样本barcode和基因ID。我们常用的矩阵也就这三种。

counts <- assay(data, "unstranded")
fpkm <- assay(data, "tpm_unstrand")
tpm <- assay(data, "tpm_unstrand")

差异分析

现在有了表达矩阵和分组文件，就可以进行最常见的差异分析了。这里以DESeq2为例。

这里的表达矩阵是我上面构建的 counts对象，分组矩阵是构建的design 对象。

# 确保counts矩阵和设计矩阵行列一致
counts <- counts[, design$sample]

# 创建 DESeq2 数据集
library(DESeq2)
dds = DESeqDataSetFromMatrix(
  countData = counts,
  colData = design,
  design = ~ group
)

# 过滤：计算至少在一半以上样本中表达的基因 (可选)
keep <- rowSums(counts(dds) >= 0) >= (ncol(counts(dds)) / 2)
dds <- dds[keep, ]

# 运行 DESeq2
dds <-DESeq(dds)

# 得到归一化的 counts 矩阵
normed <- counts(dds, normalized=TRUE)
normed <- round(normed, 1)

# 提取差异基因
res = results(dds) %>% as.data.frame()