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，不写括号，回车。
网友评论