ixxmu
commented
1 year ago 看完还不会来揍/找我 | TCGA 与 GTEx 数据库联合分析 | 附完整代码 + 注释 by 生信技能树
看完还不会来揍/找我 | TCGA 与 GTEx 数据库联合分析 | 附完整代码 + 注释
TCGA（The Cancer Genome Atlas）数据库提供了丰富的癌症组织样本数据，但由于其主要关注癌症类型，其中的正常组织样本数量有限。比如说乳腺癌，1200 多的转录组数据，其中 1100 左右都是肿瘤组织的测序数据，正常对照只有可怜巴巴一丢丢。这个时候我们就需要想办法加大正常组织样本量啦，既然 TCGA 数据库没有，我们就从其他数据库借一点过来！强烈推荐 GTEx（Genotype-Tissue Expression）数据库，因为 GTEx 涵盖了多个器官和组织中的大量正常组织样本数据。
两个数据库联合分析，怎么搞嘞！能随随便便拼接吗！那肯定不行！今天，咱们就一起来看看要怎么把它俩合并！
本文所用到的数据和代码，我已经上传到了 GitHub，有需要的话，大家可以在 https://github.com/hzyao/merge_TCGA_GTEx 进行查看下载（点击阅读原文即可直达）。
如果有朋友想知道如何上传项目到 GitHub 的话，也可以看这里：使用 git 命令行上传项目到 GitHub（以 R 包为例）
TCGA 与 GTEx 数据库联合分析
数据库介绍
TCGA（The Cancer Genome Atlas）数据库： TCGA 数据库主要关注癌症组织，它收集了来自不同癌症类型的肿瘤组织样本的基因组、转录组、蛋白质组等多维度数据。这些数据有助于研究人员深入了解不同癌症类型的分子机制、致病基因、突变频率等，从而有助于癌症的诊断、治疗和预防。
GTEx（Genotype-Tissue Expression）数据库： GTEx 数据库则关注正常组织中基因的表达模式，它收集了来自多个器官和组织的样本，分析了这些样本中基因的表达水平。通过 GTEx，研究人员可以了解不同组织中基因表达的变化，探索基因在正常生理状态下的功能和调控机制。
想详细了解 GTEx 数据库的朋友们可以看这里：GTEx 数据库 - 再也不用担心没有正常组织样本啦
数据准备
UCSC Xena 整合了来自不同研究项目、数据库和平台的生物医学数据，包括基因组学、转录组学、蛋白质组学等数据，这些数据来自于多个资源，如 TCGA、GTEx、ICGC 等。我们在这里下载后续整理也会很方便，所以咱们今天两个数据库的数据统一就在 UCSC Xena 下载！
UCSC Xena 官网：http://xena.ucsc.edu/
那我们开始下载！
首先点击Lauch Xena进入主界面（https://xenabrowser.net/）。
点击DATA SETS进入数据集界面（https://xenabrowser.net/datapages/）。
我们可以看到这个页面包含很多可下载的数据集，左侧是来自多个项目、数据库和平台的数据集列表，右侧展示了这些数据集的的来源。
TCGA 数据下载
我们今天以弥漫性大 B 细胞淋巴瘤（Lymphoid Neoplasm Diffuse Large B-cell Lymphoma）为例，它在 TCGA 中的缩写为 DLBC，所以我们点击DC TCGA Large B-cell Lymphoma (DLBC) (14 datasets)即可进入该肿瘤数据集下载页面。
里面长这个样子！
每个数据集点开都差不多，可能有的数据集包含的数据类型多一点有的少一点，如果在这里找不到你想要的，可能就需要辛苦点去原数据库官网去下载啦！
我们选择点击 gene expression RNAseq 中的HTSeq - FPKM (n=48) GDC Hub。选择 FPKM 是为了和接下来我们要下载的 GTEx 数据库中的数据计算方式保持一致（因为 GTEx 里面恰好有这种，一会儿你就知道啦！）
如果你不太清楚"Counts"、"FPKM" 和 "FPKM-UQ"有什么区别，请看这里！（该跳跳！）
"Counts"、"FPKM" 和 "FPKM-UQ" 是在基因表达数据分析中常见的三种不同的表达量计算方法，用于衡量基因的表达水平。
Counts：Counts 是一种基本的基因表达量计算方法，它直接对原始的 RNA-seq 测序数据进行处理，将每个基因在每个样本中的读数数目（即基因表达计数值）进行统计。它是通过对测序数据进行基因级别的计数来得到的，不考虑基因长度和测序深度等因素的影响。
FPKM (Fragments Per Kilobase Million)：FPKM 是一种对基因表达水平进行归一化的计算方法。它考虑了基因长度和测序深度，并将基因表达计数值调整为每千碱基对的百万个读数。FPKM 是在基因长度归一化后的表达量，使得不同长度的基因之间可以进行比较。然而，FPKM 存在对测序深度敏感的问题，在低测序深度下容易产生高估的表达量。
FPKM-UQ (FPKM-Upper Quartile normalization)：FPKM-UQ 是对 FPKM 进行进一步归一化的方法，旨在消除FPKM中的测序深度偏差。它使用上四分位数（Upper Quartile）对每个样本的FPKM值进行调整，将其调整为同一分布的中值。这样做的目的是减少测序深度偏差对 FPKM 值的影响，使得 FPKM-UQ 更适合用于样本间的比较分析。
FPKM 和 FPKM-UQ 都是早期常用的表达量计算方法，随着技术的进步和对数据的更深入理解，已经出现了更先进的表达量计算方法，比如如 TPM（Transcripts Per Million）和 TPM-UQ等等，它们对于基因表达的定量更准确和可靠。
点击HTSeq - FPKM (n=48) GDC Hub就会进入下面这个页面，点击download和ID/Gene Mapping的链接把它们都下载下来！
download：基因表达矩阵（像上图绿框框中那样）
ID/Gene Mapping：用于进行基因 ID 的转换
TCGA 中包含癌组织和癌旁组织，我们需要对它们进行区分（当然靠样本名我们也能进行区分，但是这不麻烦嘛）！所以我们下载 phenotype 中的Phenotype (n=52) GDC Hub，里面是临床信息，我们到时候用一下！顺便把生存数据（survival data (n=51) GDC Hub）也下载下来！
接下来，我们下载 GTEX 数据库中正常组织的 gene expression RNAseq 数据。
理论上 TPM 更先进，但 TCGA 只有 FPKM，为了保持数据一致性，所以我们统一下载 FPKM 数据类型。
同样也要下载它的基因表达矩阵和基因 ID 转换文件。
细心的朋友们可能已经注意到了，GTEx 的表达矩阵单位是 log2(fpkm+0.001)，和之前 TCGA 中的不一样（TCGA 中是 log2(fpkm+1)），在后续分析中，我们要记得转换成一致的！
当然，还要下载表型数据。
好！到这里，我们整合所需要的所有数据就下载完毕啦！开始合体！
数据整合
解释部分我都放在注释里啦！
本文所用到的数据和代码，我已经上传到了 GitHub，有需要的话，大家可以在 https://github.com/hzyao/merge_TCGA_GTEx 进行查看下载（点击阅读原文即可直达）。
除了这位兄台：
gtex_RSEM_gene_fpkm，它是在是太大了，大家自己去我们上面介绍的 UCSC Xena 下载好不好，或者直接去这里（https://xenabrowser.net/datapages/?dataset=gtex_RSEM_gene_fpkm&host=https%3A%2F%2Ftoil.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443），就进去它的下载页面啦！我可真贴心！）。
如果有朋友想知道如何上传项目到 GitHub 的话，也可以看这里：使用 git 命令行上传项目到 GitHub（以 R 包为例）
############################### TCGA + GTEx ####################################
# 加载要用到的包，大家如果没有可以先自己安装一下
library(data.table) # 用于高效处理大数据集的库
library(dplyr) # 数据操作和转换的库
library(tidyverse) # 数据处理和可视化的综合库
#################################### TCGA ######################################
# 读取表达矩阵
dlbc.fpkm <- fread("./download_data/TCGA-DLBC.htseq_fpkm.tsv", header = T, sep = '\t', data.table = F)
# 查看表达矩阵，我们发现第一列是ensembl id，后面的列是样本名
head(dlbc.fpkm)[1:5, 1:5]
# Ensembl_ID TCGA-RQ-A6JB-01A TCGA-FF-8046-01A TCGA-FF-A7CW-01A TCGA-RQ-A68N-01A
# 1 ENSG00000242268.2 0.00000000 0.000000 0.00000000 0.000000
# 2 ENSG00000270112.3 0.04396508 0.000000 0.01294356 0.000000
# 3 ENSG00000167578.15 3.15192692 3.399737 3.43036482 3.332282
# 4 ENSG00000273842.1 0.00000000 0.000000 0.00000000 0.000000
# 5 ENSG00000078237.5 3.03661720 3.183184 3.04236324 2.787429
# 读取基因ID转换信息，目的是为了将ensembl id转换为gene symbol
dlbc.pro <- fread("./download_data/gencode.v22.annotation.gene.probeMap", header = T, sep = '\t', data.table = F)
# 查看基因ID转换信息，我们发现第一列是ensembl id，第二列是gene symbol
head(dlbc.pro)
# id gene chrom chromStart chromEnd strand
# 1 ENSG00000223972.5 DDX11L1 chr1 11869 14409 +
# 2 ENSG00000227232.5 WASH7P chr1 14404 29570 -
# 3 ENSG00000278267.1 MIR6859-3 chr1 17369 17436 -
# 4 ENSG00000243485.3 RP11-34P13.3 chr1 29554 31109 +
# 5 ENSG00000274890.1 MIR1302-9 chr1 30366 30503 +
# 6 ENSG00000237613.2 FAM138A chr1 34554 36081 -
# 提取前两列用于进行转换
dlbc.pro <- dlbc.pro[ , c(1, 2)]
head(dlbc.pro)
# id gene
# 1 ENSG00000223972.5 DDX11L1
# 2 ENSG00000227232.5 WASH7P
# 3 ENSG00000278267.1 MIR6859-3
# 4 ENSG00000243485.3 RP11-34P13.3
# 5 ENSG00000274890.1 MIR1302-9
# 6 ENSG00000237613.2 FAM138A
# 基因ID转换信息和表达矩阵合并
dlbc.fpkm.pro <- merge(dlbc.pro, dlbc.fpkm, by.y = "Ensembl_ID", by.x = "id" )
head(dlbc.fpkm.pro)[1:5, 1:5]
# id gene TCGA-RQ-A6JB-01A TCGA-FF-8046-01A TCGA-FF-A7CW-01A
# 1 ENSG00000000003.13 TSPAN6 1.026771 0.76890645 0.35573821
# 2 ENSG00000000005.5 TNMD 0.000000 0.07154948 0.02051499
# 3 ENSG00000000419.11 DPM1 5.134474 5.01494933 5.28879599
# 4 ENSG00000000457.12 SCYL3 1.534159 1.18431712 1.91343387
# 5 ENSG00000000460.15 C1orf112 1.894523 1.72454346 2.24504933
# 去重
dlbc.fpkm.pro <- distinct(dlbc.fpkm.pro,gene, .keep_all = T)
# 把基因名转换为行名
rownames(dlbc.fpkm.pro) <- dlbc.fpkm.pro$gene
dlbc.fpkm.pro <- dlbc.fpkm.pro[ , -c(1,2)]
dim(dlbc.fpkm.pro) # 58387 48
head(dlbc.fpkm.pro)[1:5, 1:5]
# TCGA-RQ-A6JB-01A TCGA-FF-8046-01A TCGA-FF-A7CW-01A TCGA-RQ-A68N-01A TCGA-FM-8000-01A
# TSPAN6 1.026771 0.76890645 0.35573821 0.7864696 1.0575194
# TNMD 0.000000 0.07154948 0.02051499 0.0000000 0.1362624
# DPM1 5.134474 5.01494933 5.28879599 5.1281748 4.6654726
# SCYL3 1.534159 1.18431712 1.91343387 1.2004038 1.5156969
# C1orf112 1.894523 1.72454346 2.24504933 1.2145814 1.7478664
# 现在表达矩阵就构建完毕啦！
# 但是TCGA中不止有癌组织，还有癌旁组织，所以我们可以根据临床信息把它们分开
# 读取淋巴癌的临床信息
dlbc.phe <- fread("./download_data/TCGA-DLBC.GDC_phenotype.tsv", header = T, sep = '\t', data.table = F)
dlbc.phe$submitter_id.samples[1:5]
# [1] "TCGA-FA-A6HN-01A" "TCGA-GR-A4D4-01A" "TCGA-GS-A9TT-01A" "TCGA-FF-A7CQ-01A" "TCGA-FF-8046-01A"
# 把行名改成样本名
rownames(dlbc.phe) <- dlbc.phe$submitter_id.samples
# 查看临床信息的sample_type.samples列
table(dlbc.phe$sample_type.samples)
# Bone Marrow Normal Primary Tumor
# 4 48
# Primary Tumor为癌组织，我们将癌组织提取出来
dlbc.phe.t <- filter(dlbc.phe, sample_type.samples == "Primary Tumor")
# 临床信息与表达矩阵取交集
merge_phe_fpkm <- intersect(rownames(dlbc.phe.t), colnames(dlbc.fpkm.pro))
# 提取癌组织的表达矩阵
dlbc.exp <- dlbc.fpkm.pro[ , merge_phe_fpkm]
dim(dlbc.exp)
# [1] 58387 48
# TCGA DLBC 癌组织样本表达矩阵构建完成，保存一下！
saveRDS(dlbc.exp, file='./generated_data/dlbc.exp.rds')
#################################### GTEx ######################################
# 读取gtex的表达矩阵，解压和不解压其实都可以读取（如果超级慢或者出现error，不要慌，据说可能因为内存小）
gtex.exp <- fread("./download_data/gtex_RSEM_gene_fpkm", header = T, sep = '\t', data.table = F)
gtex.exp[1:5, 1:4]
# sample GTEX-S4Q7-0003-SM-3NM8M GTEX-QV31-1626-SM-2S1QC GTEX-13QIC-0011-R1a-SM-5O9CJ
# 1 ENSG00000242268.2 -3.8160 -9.9658 -9.9658
# 2 ENSG00000259041.1 -9.9658 -9.9658 -9.9658
# 3 ENSG00000270112.3 -4.0350 -2.9324 -2.1779
# 4 ENSG00000167578.16 4.2958 3.9204 6.1528
# 5 ENSG00000278814.1 -9.9658 -9.9658 -9.9658
dim(gtex.exp)
# [1] 60498 7863
# 注意这里包含GTEx中所有正常组织来源的样本，我们保存一下，后面继续提取自己需要的组织样本
# saveRDS(gtex.exp, file = './generated_data/gtex.exp.rds')
# 读取gtex的基因ID转换信息，将ensembl id转换为gene symbol
gtex.pro <- fread("./download_data/probeMap_gencode.v23.annotation.gene.probemap", header = T, sep = '\t', data.table = F)
head(gtex.pro)
# id gene chrom chromStart chromEnd strand
# 1 ENSG00000223972.5 DDX11L1 chr1 11869 14409 +
# 2 ENSG00000227232.5 WASH7P chr1 14404 29570 -
# 3 ENSG00000278267.1 MIR6859-1 chr1 17369 17436 -
# 4 ENSG00000243485.3 RP11-34P13.3 chr1 29554 31109 +
# 5 ENSG00000274890.1 MIR1302-2 chr1 30366 30503 +
# 6 ENSG00000237613.2 FAM138A chr1 34554 36081 -
dim(dlbc.pro) # 60483 2
dim(gtex.pro) # 60498 6
# 我们顺便发现了v23和v22的差异有多少，60498 vs. 60483
# 提取前两列用于基因ID转换
gtex.pro <- gtex.pro[, c(1,2)]
# 基因ID转换信息和表达矩阵合并
gtex.fpkm.pro <- merge(gtex.pro, gtex.exp, by.y ="sample", by.x = "id" )
# 取交集，为了决定之后gtex和dlbc合并是按照gene symbol还是ensembl id进行
length(intersect(gtex.pro$id, dlbc.pro$id)) # 42566
length(intersect(rownames(dlbc.exp), gtex.fpkm.pro$gene)) # 57793
# 按照gene symbol合并会有57793个，按照ensembl id合并有42566个，所以这次我们选择按照gene symbol进行合并
# 接下来，提取正常淋巴组织的表达矩阵，根据gtex的临床信息来匹配淋巴组织的样本
# 读取gtex的临床信息
gtex.phe <- fread("./download_data/GTEX_phenotype", header = T, sep = '\t', data.table = F)
rownames(gtex.phe) <- gtex.phe$Sample
# 给它改一下列名，好看点
colnames(gtex.phe) <- c("Sample", "body_site_detail (SMTSD)", "primary_site", "gender", "patient", "cohort")
table(gtex.phe$primary_site)
# <not provided> Adipose Tissue Adrenal Gland Bladder Blood Blood Vessel Bone Marrow
# 5 621 161 13 595 753 102
# Brain Breast Cervix Uteri Colon Esophagus Fallopian Tube Heart
# 1426 221 11 384 805 7 493
# Kidney Liver Lung Muscle Nerve Ovary Pancreas
# 38 141 381 478 335 112 203
# Pituitary Prostate Salivary Gland Skin Small Intestine Spleen Stomach
# 126 122 71 977 106 121 209
# Testis Thyroid Uterus Vagina
# 208 366 93 99
# 上面展示的是GTEx中所有的组织来源，今天我们分析的是淋巴癌
# TCGA中数据集名为DLBC，在GTEx中对应的组织来源应该为blood，上面显示有595个
# 在本文最后，我附上了TCGA 与 GTEx 癌症类型与组织来源对应表，大家有需要可以去看哟！
# 筛选出blood组织来源的样本
gtex.phe.s <- filter(gtex.phe, primary_site == "Blood")
# 临床信息与表达矩阵取交集
merge_phe_fpkm_gtex <- intersect(rownames(gtex.phe.s), colnames(gtex.fpkm.pro)) # 444
gtex.s <- gtex.fpkm.pro[ , c("gene", merge_phe_fpkm_gtex)]
# 去重
gtex.s <- distinct(gtex.s, gene, .keep_all = T)
rownames(gtex.s) <- gtex.s$gene
gtex.s <- gtex.s[ , -1]
dim(gtex.s)
# [1] 58581 444
# 我们发现GTEx中有444个blood样本
# gtex的表达矩阵是按照log2(fpkm+0.001)处理的，dlbc是按照log2(fpkm+1)
# 所以合并之前，我们需要把它们的处理方式调整为相同的
gtex.s2 <- 2^gtex.s
gtex.s3 <- log2(gtex.s2-0.001+1)
# 现在数据的处理方式都相同，就有可比性啦！
# 正式合并开始！
# 把gtex的blood组织数据与tcga的dlbc数据合并
all.data <- merge(gtex.s3, dlbc.exp, by = 0)
all.data <- column_to_rownames(all.data, "Row.names")
dim(all.data)
# [1] 57793 492
head(all.data)[1:5, 1:4]
# GTEX-111YS-0006-SM-5NQBE GTEX-1122O-0003-SM-5Q5DL GTEX-113IC-0006-SM-5NQ9C GTEX-113JC-0006-SM-5O997
# 5_8S_rRNA -1.571525e-08 -1.571525e-08 -1.571525e-08 -1.571525e-08
# 5S_rRNA -1.571525e-08 -1.571525e-08 -1.571525e-08 -1.571525e-08
# 7SK -1.571525e-08 -1.571525e-08 -1.571525e-08 -1.571525e-08
# A1BG 2.641511e+00 3.407365e+00 3.005374e+00 1.545996e+00
# A1BG-AS1 1.580150e+00 2.618198e+00 2.657654e+00 6.690156e-01
# 前面444列是来自GTEx的正常blood组织样本，后面48列是来自TCGA的DLBC癌组织样本
# 浅浅保存一下
saveRDS(all.data, file = './generated_data/all_data_tcga_gtex.rds')
############################### 整合完毕 #######################################
# 去除批次效应
library(limma)
nromalized.data <- normalizeBetweenArrays(all.data)
nromalized.data <- as.data.frame(nromalized.data)
# 接下来，就可以进行我们的后续分析啦！
TCGA 与 GTEx 癌症类型与组织来源对应表
| TCGA | Cancer | Tumor | Normal | GTEx | Numer |
|---|---|---|---|---|---|
| ACC | Adrenocortical carcinoma | 77 | - | Adrenal Gland | 128 |
| BLCA | Bladder Urothelial Carcinoma | 404 | 19 | Bladder | 9 |
| BRCA | Breast invasive carcinoma | 1085 | 112 | Breast | 179 |
| CESC | Cervical squamous cell carcinoma and endocervical adenocarcinoma | 306 | 3 | Cervix Uteri | 10 |
| CHOL | Cholangio carcinoma | 36 | 9 | - | - |
| COAD | Colon adenocarcinoma | 275 | 41 | Colon | 308 |
| DLBC | Lymphoid Neoplasm Diffuse Large B-cell Lymphoma | 47 | - | Blood | 337 |
| ESCA | Esophageal carcinoma | 182 | 13 | Esophagus | 273 |
| GBM | Glioblastoma multiforme | 163 | - | Brain | 207 |
| HNSC | Head and Neck squamous cell carcinoma | 519 | 44 | - | - |
| KICH | Kidney Chromophobe | 66 | 25 | Kidney | 28 |
| KIRC | Kidney renal clear cell carcinoma | 523 | 72 | Kidney | 28 |
| KIRP | Kidney renal papillary cell carcinoma | 286 | 32 | Kidney | 28 |
| LAML | Acute Myeloid Leukemia | 173 | - | Bone Marrow | 70 |
| LGG | Brain Lower Grade Glioma | 518 | - | Brain | 207 |
| LIHC | Liver hepatocellular carcinoma | 369 | 50 | Liver | 110 |
| LUAD | Lung adenocarcinoma | 483 | 59 | Lung | 288 |
| LUSC | Lung squamous cell carcinoma | 486 | 50 | Lung | 288 |
| MESO | Mesothelioma | 87 | - | - | - |
| OV | Ovarian serous cystadenocarcinoma | 426 | - | Ovary | 88 |
| PAAD | Pancreatic adenocarcinoma | 179 | 4 | Pancreas | 167 |
| PCPG | Pheochromocytoma and Paraganglioma | 182 | 3 | - | - |
| PRAD | Prostate adenocarcinoma | 492 | 52 | Prostate | 100 |
| READ | Rectum adenocarcinoma | 92 | 10 | Colon | 308 |
| SARC | Sarcoma | 262 | 2 | - | - |
| SKCM | Skin Cutaneous Melanoma | 461 | 1 | Skin | 557 |
| STAD | Stomach adenocarcinoma | 408 | 36 | Stomach | 175 |
| TGCT | Testicular Germ Cell Tumors | 137 | - | Testis | 165 |
| THCA | Thyroid carcinoma | 512 | 59 | Thyroid | 278 |
| THYM | Thymoma | 118 | 2 | Blood | 337 |
| UCEC | Uterine Corpus Endometrial Carcinoma | 174 | 13 | Uterus | 78 |
| UCS | Uterine Carcinosarcoma | 57 | - | Uterus | 78 |
| UVM | Uveal Melanoma | 79 | - | - | - |
参考资料
https://www.jianshu.com/p/dd5806d4678a https://zhuanlan.zhihu.com/p/83140956 https://www.jingege.wang/2022/03/16/gtex联合tcga数据库差异分析/
Sunyaxin666
https://mp.weixin.qq.com/s/_xVOBTJ0DkpxiIj03NY6TQ