查看原文
其他

TCGA-5.GDC数据整理-后续(含情感体验)

豆豆花花 生信星球 2022-06-06

 今天是生信星球陪你的第516天


   大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~

   就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~

   这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!

花花写于2020-01-07

本系列是我的TCGA学习记录,跟着生信技能树B站课程学的,已获得授权(嗯,真的^_^)。课程链接:https://www.bilibili.com/video/av49363776

目录:

TCGA-1.数据下载

TCGA2.GDC数据整理

TCGA3.R包TCGA-biolinks下载数据

TCGA-4.使用RTCGA包获取数据

书接第二回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语言功底和必胜的信念(手动狗头)。

插个小广告!

生信零基础入门学习小组长期报名中

GEO数据挖掘广州专场课程

再给生信技能树打个call!

全国巡讲第21站(长沙线下培训)

全球公益巡讲招学徒

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存