让TCGA表达矩阵和临床信息方便的匹配

0.背景

关于样本id和分组信息：TCGA的样本id里藏着分组信息

TCGA数据库的表达矩阵是按照样本来组织的，而临床信息是按照病人(个体)来组织的。

病人数量与样本数量并非一一对应

原因挺多的，我罗列了一下：

1.有的病人有tumor、normal两个样本，有的病人只有tumor一个样本。
2.有的病人会取两个及以上tumor样本，比如一个原发一个转移，或者一个冷冻样本一个石蜡包埋。
3.有的病人没有被记录临床信息（少，但真的有）
4.有的临床信息表格里的病人没有RNA-seq样本

前三个都好理解，第四个可能有点懵逼。这是因为TCGA是多组学数据库，同一个癌症都是同一帮病人，不管分析什么组学数据，临床信息都是一样的，所以那些没有RNA-seq样本但有其他组学的样本的病人，也有临床信息。

2.解决的思路

1.以病人为中心，如果找临床信息是为了生存分析用，那么可以剔除normal样本，tumor样本去重。得到的是每个病人的临床信息，样本与病人一一对应，不存在重复。

2.以样本为中心，如果不想去重，就把临床信息表格加上一列样本id，即以样本为中心，得到的是每个样本的临床信息。每个病人会有一个或多个样本。

3.代码实现

2.1示例数据

我放在tinyarray包里了，安装好了包就能使用示例数据。

if(!require(tinyarray))devtools::install_github("xjsun1221/tinyarray")
library(tinyarray)
exp_hub1[1:4,1:4]

##        GTEX-S33H-1226-SM-4AD69 GTEX-VJYA-0826-SM-4KL1M TCGA-FB-A545-01
## CXCL8                 34.00041                50.99857        1341.029
## FN1                  591.01914               960.00492      137998.164
## COL3A1              2226.94824              3682.08761      177553.205
## ISG15                125.00175               135.00350        1784.980
##        GTEX-ZF3C-2026-SM-4WWB5
## CXCL8                 388.9918
## FN1                  5218.1531
## COL3A1              15621.0916
## ISG15                 513.9898

meta1[1:4,1:4]

##               sample event    X_PATIENT time
## 216 TCGA-3A-A9IO-01A     0 TCGA-3A-A9IO 1942
## 172 TCGA-US-A774-01A     1 TCGA-US-A774  695
## 128 TCGA-HZ-A49H-01A     0 TCGA-HZ-A49H  491
## 37  TCGA-FB-A4P5-01A     1 TCGA-FB-A4P5  179

exp_hub1不只这一个用处，里面有一些gtex样本，可以去掉，其实留着也没事。

library(tidyverse)
k = str_starts(colnames(exp_hub1),"GTEX");table(k)

## k
## FALSE  TRUE 
##   183   167

exp = exp_hub1[,!k]
table(make_tcga_group(exp))

## 
## normal  tumor 
##      4    179

可以看到有四个正常样本。

我写成函数了，装最新版本的tinyarray就能用了

3.2 以病人为中心

去掉正常样本，每个病人只保留一个肿瘤样本。

match_exp_cl(exp,meta1[,2:4],"X_PATIENT",sample_centric = F)
dim(exp_matched)

## [1]   8 176

dim(cl_matched)

## [1] 176   4

table(make_tcga_group(exp_matched))

## 
## normal  tumor 
##      0    176

3.3 以样本为中心

不去重，每个病人会保留一个或多个样本。

match_exp_cl(exp,meta1[,2:4],"X_PATIENT")
dim(exp_matched)

## [1]   8 181

dim(cl_matched)

## [1] 181   4

table(make_tcga_group(exp_matched))

## 
## normal  tumor 
##      4    177

想明白为什么表达矩阵有183列，病人却只有176个了么 想明白为什么表达矩阵有183列，匹配后却只有181列了么

length(unique(str_sub(colnames(exp),1,12)))

## [1] 178

table(unique(str_sub(colnames(exp),1,12)) %in% meta1$X_PATIENT)

## 
## FALSE  TRUE 
##     2   176

想明白了就回家吃饭吧~都是R语言技巧，如果想看源代码，可以在你的R里只输入match_exp_cl，不写括号，回车。

插个小广告！

生信零基础入门学习小组

生信入门班（四周线上直播课，长期开班）

数据挖掘班（医生/医学生首选，三周线上直播课，长期开班）

一起来学单细胞吗？

生信星球答疑公告