TCGA-5.GDC数据整理-后续(含情感体验)
大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~
就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~
这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!
花花写于2020-01-07
本系列是我的TCGA
学习记录,跟着生信技能树B站课程学的,已获得授权(嗯,真的^_^)。课程链接:https://www.bilibili.com/video/av49363776
目录:
书接第二回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里并简化
1meta <- jsonlite::fromJSON("metadata.cart.2020-01-08.json")
2colnames(meta)
3# [1] "md5sum" "data_type" "file_name" "file_size"
4# [5] "data_format" "analysis" "submitter_id" "access"
5# [9] "state" "file_id" "data_category" "associated_entities"
6# [13] "experimental_strategy" "annotations"
7#找了大半天的TCGA-ID就在associated_entities里,说出来你可能不信,data.frame的一个格子里竟然可以装列表,我今天也是第一次见。
8entity <- meta$associated_entities
9meta$associated_entities[[1]]
10# entity_id case_id entity_submitter_id
11#1 81b91c93-e57c-4aa1-a19d-6b004ff539e7 75a0bb0b-6528-4fb8-a9e6-254905d21df4 TCGA-#MN-A4N1-01A-11H-A24S-13
12# entity_type
13#1 aliquot
14class(entity)
15#[1] "list"
16#可以看到这个列表的第三个元素就是tcgaid
17jh = function(x){
18 as.character(x[3])
19}
20jh(entity[[1]])
21#[1] "TCGA-MN-A4N1-01A-11H-A24S-13"
22ID = sapply(entity,jh)
23options(stringsAsFactors = F)
24file2id = data.frame(file_name = meta$file_name,
25 ID = ID)
上一篇用到的变量mis就是读取文件的顺序,和行名顺序一致。
1head(mis)
2# [1] "00537534-e079-41dd-8cb9-ae1b0a27ad24/49d56cba-70c8-47a1-a890-9784ccfbb156.mirbase21.mirnas.quantification.txt"
3# [2] "00c46e8b-f303-4d20-bd6d-d650c36895f5/af29b644-b3a8-455b-9f07-b956d41f6ec4.mirbase21.mirnas.quantification.txt"
4# [3] "00fe6d87-6f9f-4818-8e9d-9672af400474/853e5acf-2710-4d13-989f-d977f5d3db51.mirbase21.mirnas.quantification.txt"
5# [4] "0104295c-9b03-45c0-a598-36d9b788b450/bbb0a08f-8a93-420b-b96c-58bb3e07c9e6.mirbase21.mirnas.quantification.txt"
6# [5] "0178f4cc-26db-47ce-8025-6fec60523991/33568f1e-8704-4ea5-948b-564057e5a36b.mirbase21.mirnas.quantification.txt"
7# [6] "023b0bf9-3954-439b-8230-689e11ebd39a/9fb0324a-992d-4d1e-9f04-76d307ccb754.mirbase21.mirnas.quantification.txt"
8
9mis2 = stringr::str_split(mis,"/",simplify = T)[,2]
10mis2[1] %in% file2id$file_name
11##[1] TRUE
下面涉及到一个比较绕的知识点,就是向量重排序,需要深入理解match函数了。详见:向量排序:match,rank,order,sort,rev
1head(match(mis2,file2id$file_name))
2## [1] 143 348 119 500 479 45
3#排序可以验证一下,match返回值的第一个是143,意思mis2的第一个元素是file2id$file_name的第143个元素。
4mis2[1] == file2id$file_name[143]
5##[1] TRUE
6row_tcga = file2id[match(mis2,file2id$file_name),]
7rownames(mi_df) = row_tcga$ID
8mi_df[1:4,1:4]
9## hsa-let-7a-1 hsa-let-7a-2 hsa-let-7a-3 hsa-let-7b
10## TCGA-MN-A4N4-01A-12H-A24S-13 43563 44020 43680 65476
11## TCGA-55-7816-01A-11H-2169-13 75848 75451 75769 54045
12## TCGA-78-7160-01A-11H-2038-13 56189 56527 56437 44829
13## TCGA-05-4382-01A-01T-1207-13 7956 7951 7834 15317
就只是加个行名,背后需要强大的R语言功底和必胜的信念(手动狗头)。
插个小广告!
再给生信技能树打个call!