書接第二回TCGA-2.GDC數(shù)據(jù)整理，留下了一個大坑，GDC下載的數(shù)據(jù)整理出來表達矩陣是缺少行名的。
手癢的我要動手了。
加行名并不容易，我搗鼓了好久啊。
思路
1.表達矩陣的行名應(yīng)該是TCGA的樣本id（id是行名還是列名無所謂啦，t()一個函數(shù)搞定）
2.要找到下載的mirna矩陣單個文件名和TCGA的樣本id的對應(yīng)關(guān)系（難難難哭了）
內(nèi)心掙扎思考記錄：
很堅定，這個對應(yīng)關(guān)系文件一定是存在的，否則沒有行名的表達是沒辦法分析的，gdc作為官方下載工具絕不可能沒救。其次，這個對應(yīng)關(guān)系一定是一對一，表達矩陣是567行，id也應(yīng)該是567個，否則就會丟失數(shù)據(jù)。
所以這個對應(yīng)關(guān)系不在我們整理的xml信息里，因為xml只有522個病人的信息，567個樣本是因為有的病人既有tumor樣本又有normal樣本。TCGA-biolinks也是基于gdc做了更多開發(fā)，用它下載的表達矩陣就有行名，一個樣本都不少。所以這個對應(yīng)關(guān)系一定能從gdc網(wǎng)頁上找到。
嘗試
1.在gdc網(wǎng)頁樣本選擇頁面找
我把這個json下載下來研究了半天，也沒發(fā)現(xiàn)任何id的影子。
2.在其他地方各種搜啊找，折騰半天找不到。
我為什么一定要找到它？
此時我的心情已經(jīng)非常暴躁了，中午12點半了，都不想吃飯。雖然TCGA-biolinks能夠更方便的實現(xiàn)，但它不是官方工具（萬一哪天。。。嗯），而且gdc我已經(jīng)學了，發(fā)兩篇筆記了，它又是官方玩法，干到一半不能撂挑子不管啊，我就是這么執(zhí)著的人。
3.藏在購物車里了
估計只有開上帝視角才能直接找到吧，全靠執(zhí)念。
把選出來的樣本全都加進購物車，點開購物車就會發(fā)現(xiàn)，這里有一排下載按鈕。挨個下載了一下，發(fā)現(xiàn)
metadata下載的就是一個json文件（"metadata.cart.2020-01-08.json"，自己從購物車下載即可），讀進R里就是567行！
開始進入R語言的表演
接下來的代碼是上一篇TCGA-2.GDC數(shù)據(jù)整理中的代碼為基礎(chǔ)繼續(xù)寫的，運行完上一篇的代碼接著運行這一篇的。
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函數(shù)了。詳見：http://www.itdecent.cn/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語言功底和必勝的信念（手動狗頭）。
微信公眾號生信星球同步更新我的文章，歡迎大家掃碼關(guān)注！
我們有為生信初學者準備的學習小組,點擊查看??
想要參加我的線上線下課程，也可加好友咨詢??
如果需要提問，請先看生信星球答疑公告
