TCGA数据库| 如何将表达矩阵与样本临床数据进行合并？

生信技能树 · 公众号 · · 2025-01-07 20:42

正文

前面我们已经给大家介绍过TCGA数据库中样本barcode的详细组成： TCGA样本barcode详细介绍，现在我们来看看如何将基因表达矩阵与样本临床信息进行合并，方便后续做比如生存分析，基因在不同样本分期、性别、年龄分组等中的差异表达情况。

首先我们去TGCA下载如乳腺癌的基因表达矩阵

这里使用R包 TCGAbiolinks 去TCGA官网下载数据。

1、加载包：

## download tcga data
## update: 2024-02-22
## Author: zhang juan

rm(list=ls())
# 当然，需要先去安装这个包，如果已安装就可以跳过：
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("TCGAbiolinks")

## load packages
library(TCGAbiolinks)
library(SummarizedExperiment)
suppressPackageStartupMessages(library(tidyverse))

2、癌症类型选择：

# 癌症类型，用 getGDCprojects()$project_id 查看所有
getGDCprojects()$project_id

# [1] "TARGET-AML"                "MATCH-Z1I"                 "HCMI-CMDC"                 "MATCH-W"         
# [5] "MATCH-Z1D"                 "MATCH-Z1A"                 "MATCH-U"                   "MATCH-Q"         
# [9] "TCGA-PCPG"                 "MATCH-H"                   "MATCH-C1"                  "TCGA-THYM"       
# [13] "MATCH-I"                   "MATCH-S1"                  "MATCH-P"                   "MATCH-R"      
# [17] "MATCH-Z1B"                 "TCGA-PAAD"                 "TCGA-STAD"                 "TCGA-TGCT"   
# [21] "MATCH-S2"                  "TCGA-SARC"                 "TCGA-PRAD"                 "TCGA-READ"     
# [25] "TCGA-UCS"                  "TCGA-UVM"                  "TRIO-CRU"                  "VAREPOP-APOLLO"
# [29] "WCDT-MCRPC"                "TARGET-ALL-P1"             "REBC-THYR"                 "TARGET-ALL-P2"   
# ...

不同缩写代表的含义可取这个地址查看： https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/bcr-batch-codes

本次乳腺癌缩写为： BRCA

3、下载：

# 设置query 
query   project = "TCGA-BRCA",  # 癌症类型，用 getGDCprojects()$project_id 查看所有
  data.category="Transcriptome Profiling",     # 数据类别, 用getProjectSummary(project)查看所有类别
  data.type ="Gene Expression Quantification", # 数据类型
  workflow.type="STAR - Counts"                # 工作流类型
  )

## 下载数据
GDCdownload(query=query, files.per.chunk= 50, directory = "./")

下来后的数据为一个样本一个tsv文件：如 8d1641ea-7552-4d23-9298-094e0056386a.rna_seq.augmented_star_gene_counts.tsv

4、整合成一个表达矩阵：

## 整理数据并存储为 R对象
GDCprepare(query,save=T,save.filename="TCGA-BRCA.transcriptome.Rdata", directory = "./")

## 再次加载
load("TCGA-BRCA.transcriptome.Rdata")
ls()
names(assays(data))
rowdata

5、提取mRNA的SummarizedExperiment对象，根据 gene_type 取子集，太简单了！

table(rowdata$gene_type)

tcga_mrna $gene_type == "protein_coding",]

tcga_mrna_count "unstranded")   # mRNA的counts矩阵
tcga_mrna_tpm "tpm_unstrand")  # mRNA的tpm矩阵
tcga_mrna_fpkm "fpkm_unstrand") # mRNA的fpkm矩阵


# 添加gene_symbol, 先提取gene_name
symbol_mrna $gene_name
head(symbol_mrna)

####################################################### count值
tcga_mrna_count_symbol 
# 去重复保留最大的那个
tcga_mrna_count_symbol1 % 
  as_tibble() %>% # tibble不支持row name，我竟然才发现！
  mutate(meanrow = rowMeans(.[,-1]), .before=2) %>% 
  arrange(desc(meanrow)) %>% 
  distinct(symbol_mrna,.keep_all=T) %>% 
  select(-meanrow)

saveRDS(tcga_mrna_count_symbol1, file = "tcga_mrna_count_symbol.rds")
write.table(tcga_mrna_count_symbol1, file ="tcga_mrna_count_symbol.xls",row.names = F,sep = "\t",quote = F)

####################################################### tpm值
tcga_mrna_tpm_symbol 
# 去重复保留最大的那个
tcga_mrna_tpm_symbol1 % 
  as_tibble() %>% # tibble不支持row name，我竟然才发现！
  mutate(meanrow = rowMeans(.[,-1]), .before=2) %>% 
  arrange(desc(meanrow)) %>% 
  distinct(symbol_mrna,.keep_all=T) %>% 
  select(-meanrow)

saveRDS(tcga_mrna_tpm_symbol1, file = "tcga_mrna_tpm_symbol.rds")
write.table(tcga_mrna_tpm_symbol1, file = "tcga_mrna_tpm_symbol.xls",row.names = F,sep = "\t",quote = F)

####################################################### fpkm值
tcga_mrna_fpkm_symbol 
# 去重复保留最大的那个
tcga_mrna_fpkm_symbol1 % 
  as_tibble() %>% # tibble不支持row name，我竟然才发现！
  mutate(meanrow = rowMeans(.[,-1]), .before=2) %>% 
  arrange(desc(meanrow)) %>% 
  distinct(symbol_mrna,.keep_all=T) %>% 
  select(-meanrow)

saveRDS(tcga_mrna_fpkm_symbol1, file = "tcga_mrna_fpkm_symbol.rds")
write.table(tcga_mrna_fpkm_symbol1, file = "tcga_mrna_fpkm_symbol.xls",row.names = F,sep = "\t",quote = F)

接着下载样本临床信息

1、同样首先需要联网进行 query：

##############################################################################
########################## 3.批量下载临床数据 ###################################
##############################################################################
# ref: https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/clinical.html

query   project = "TCGA-BRCA", 
  data.category = "Clinical",
  data.format = "bcr xml"
  )
save(query, file = "TCGA-BRCA.clinic.query.rdata")

# 下载到当前目录
GDCdownload(query, files.per.chunk= 50, directory = "./")

2、对下载的数据进行整理：

clinical "patient", directory = "./")
clinical.follow_up "follow_up", directory = "./")
clinical.stage_event "stage_event", directory = "./")
clinical.drug "drug", directory = "./")
clinical.radiation "radiation", directory = "./")
clinical.admin "admin", directory = "./")

# 保存
saveRDS(clinical, file = "TCGA-BRCA.clinical_patient.rds")
saveRDS(clinical.admin, file = "TCGA-BRCA.clinical_admin.rds")
saveRDS(clinical.drug, file = "TCGA-BRCA.clinical_drug.rds")
saveRDS(clinical.follow_up, file = "TCGA-BRCA.clinical_follow_up.rds")
saveRDS(clinical.radiation, file = "TCGA-BRCA.clinical_radiation.rds")
saveRDS(clinical.stage_event, file = "TCGA-BRCA.clinical_stage_event.rds")

现在将基因表达矩阵与临床信息整合在一起

先看看各自的样本ID名，根据前面的介绍《 TCGA样本barcode详细介绍》，可以看到表达矩阵里面的是样本ID ，临床信息中是patient ID ，一个病人可能会取多个样本，比如同时存在正常样本与肿瘤样本，也可能同时具有好几个肿瘤样本：

# 表达矩阵 样本名
mrna_fpkm "tcga_mrna_fpkm_symbol.rds")
head(colnames(mrna_fpkm))

# [1] "symbol_mrna"                  "TCGA-5L-AAT0-01A-12R-A41B-07" "TCGA-A2-A04U-01A-11R-A115-07" "TCGA-AN-A04A-01A-21R-A034-07"
# [5] "TCGA-A7-A13D-01A-13R-A12P-07" "TCGA-BH-A201-01A-11R-A14M-07"

# 临床信息
clinical "TCGA-BRCA.clinical_patient.rds")
colnames(clinical)
head(clinical[,1:6])
# 我们后面相比较不同病理分期间某个基因表达差异，这里过滤一下样本
clinical "bcr_patient_barcode", "stage_event_pathologic_stage")]
colnames(clinical) "bcr_patient_barcode", "pathologic_stage")
str(clinical)
table(clinical$pathologic_stage)
clinical$pathologic_stage $pathologic_stage)
clinical $pathologic_stage!="",]
clinical head(clinical)

# 变成 stage I 、II、III、IV、
clinical$stage $pathologic_stage
clinical$stage[grepl("Stage I$|Stage IA$|Stage IB$",clinical$pathologic_stage)] "Stage I"
clinical$stage[grepl("Stage II$|Stage IIA$|Stage IIB$",clinical$pathologic_stage)] "Stage II"
clinical$stage[grepl("Stage III$|Stage IIIA$|Stage IIIB$|Stage IIIC$",clinical$pathologic_stage)] "Stage III"
table(clinical$stage)
table(clinical$pathologic_stage,clinical$stage)
clinical$stage $stage, levels = c("Stage I","Stage II","Stage III","Stage IV"))

那么，这里对应的时候，一般可以先将样本分为肿瘤样本与正常样本，看看肿瘤样本中某个基因表达的高低分组生存曲线KM差异：

肿瘤样本的编号一般为样本 barcode中的第14-15位编码字符：

01-09为肿瘤样本，10以及10以上的为对照样本。肿瘤样本里面又有很多细小的分类：

https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes

我们这里直接提取 01A类的实体瘤样本：

# 提取 01A类的实体瘤样本