TCGA-3.GDC数据整理-后续（含情感体验）

2020-01-08 本文已影响0人小洁忘了怎么分身

书接第二回TCGA-2.GDC数据整理，留下了一个大坑，GDC下载的数据整理出来表达矩阵是缺少行名的。

手痒的我要动手了。

加行名并不容易，我捣鼓了好久啊。

思路

1.表达矩阵的行名应该是TCGA的样本id（id是行名还是列名无所谓啦，t()一个函数搞定）
2.要找到下载的mirna矩阵单个文件名和TCGA的样本id的对应关系（难难难哭了）

内心挣扎思考记录：
很坚定，这个对应关系文件一定是存在的，否则没有行名的表达是没办法分析的，gdc作为官方下载工具绝不可能没救。

其次，这个对应关系一定是一对一，表达矩阵是567行，id也应该是567个，否则就会丢失数据。
所以这个对应关系不在我们整理的xml信息里，因为xml只有522个病人的信息，567个样本是因为有的病人既有tumor样本又有normal样本。

TCGA-biolinks也是基于gdc做了更多开发，用它下载的表达矩阵就有行名，一个样本都不少。所以这个对应关系一定能从gdc网页上找到。

尝试

1.在gdc网页样本选择页面找

我把这个json下载下来研究了半天，也没发现任何id的影子。

2.在其他地方各种搜啊找，折腾半天找不到。

我为什么一定要找到它？
此时我的心情已经非常暴躁了，中午12点半了，都不想吃饭。虽然TCGA-biolinks能够更方便的实现，但它不是官方工具（万一哪天。。。嗯），而且gdc我已经学了，发两篇笔记了，它又是官方玩法，干到一半不能撂挑子不管啊，我就是这么执着的人。

3.藏在购物车里了

估计只有开上帝视角才能直接找到吧，全靠执念。

把选出来的样本全都加进购物车，点开购物车就会发现，这里有一排下载按钮。挨个下载了一下，发现metadata下载的就是一个json文件（"metadata.cart.2020-01-08.json"，自己从购物车下载即可），读进R里就是567行！

开始进入R语言的表演

接下来的代码是上一篇TCGA-2.GDC数据整理中的代码为基础继续写的，运行完上一篇的代码接着运行这一篇的。

1.将json读进R里并简化

meta <- jsonlite::fromJSON("metadata.cart.2020-01-08.json")
colnames(meta)
# [1] "md5sum"                "data_type"             "file_name"             "file_size"            
# [5] "data_format"           "analysis"              "submitter_id"          "access"               
# [9] "state"                 "file_id"               "data_category"         "associated_entities"  
# [13] "experimental_strategy" "annotations" 
#找了大半天的TCGA-ID就在associated_entities里，说出来你可能不信，data.frame的一个格子里竟然可以装列表，我今天也是第一次见。
entity <- meta$associated_entities
meta$associated_entities[[1]]
#                             entity_id                              case_id          entity_submitter_id
#1 81b91c93-e57c-4aa1-a19d-6b004ff539e7 75a0bb0b-6528-4fb8-a9e6-254905d21df4 TCGA-#MN-A4N1-01A-11H-A24S-13
#  entity_type
#1     aliquot
class(entity)
#[1] "list"
#可以看到这个列表的第三个元素就是tcgaid
jh = function(x){
  as.character(x[3])
}
jh(entity[[1]])
#[1] "TCGA-MN-A4N1-01A-11H-A24S-13"
ID = sapply(entity,jh)
options(stringsAsFactors = F)
file2id = data.frame(file_name = meta$file_name,
                     ID = ID)

上一篇用到的变量mis就是读取文件的顺序，和行名顺序一致。

head(mis)
# [1] "00537534-e079-41dd-8cb9-ae1b0a27ad24/49d56cba-70c8-47a1-a890-9784ccfbb156.mirbase21.mirnas.quantification.txt"
# [2] "00c46e8b-f303-4d20-bd6d-d650c36895f5/af29b644-b3a8-455b-9f07-b956d41f6ec4.mirbase21.mirnas.quantification.txt"
# [3] "00fe6d87-6f9f-4818-8e9d-9672af400474/853e5acf-2710-4d13-989f-d977f5d3db51.mirbase21.mirnas.quantification.txt"
# [4] "0104295c-9b03-45c0-a598-36d9b788b450/bbb0a08f-8a93-420b-b96c-58bb3e07c9e6.mirbase21.mirnas.quantification.txt"
# [5] "0178f4cc-26db-47ce-8025-6fec60523991/33568f1e-8704-4ea5-948b-564057e5a36b.mirbase21.mirnas.quantification.txt"
# [6] "023b0bf9-3954-439b-8230-689e11ebd39a/9fb0324a-992d-4d1e-9f04-76d307ccb754.mirbase21.mirnas.quantification.txt"

mis2 = stringr::str_split(mis,"/",simplify = T)[,2]
mis2[1] %in% file2id$file_name
##[1] TRUE

下面涉及到一个比较绕的知识点，就是向量重排序，需要深入理解match函数了。详见：https://www.jianshu.com/p/335dcb2cd27a。

head(match(mis2,file2id$file_name))
## [1] 143 348 119 500 479  45
#排序可以验证一下，match返回值的第一个是143，意思mis2的第一个元素是file2id$file_name的第143个元素。
mis2[1] == file2id$file_name[143]
##[1] TRUE
row_tcga = file2id[match(mis2,file2id$file_name),]
rownames(mi_df) = row_tcga$ID
mi_df[1:4,1:4]
##                             hsa-let-7a-1 hsa-let-7a-2 hsa-let-7a-3 hsa-let-7b
## TCGA-MN-A4N4-01A-12H-A24S-13        43563        44020        43680      65476
## TCGA-55-7816-01A-11H-2169-13        75848        75451        75769      54045
## TCGA-78-7160-01A-11H-2038-13        56189        56527        56437      44829
## TCGA-05-4382-01A-01T-1207-13         7956         7951         7834      15317

就只是加个行名，背后需要强大的R语言功底和必胜的信念（手动狗头）。