以前都是直接拿前體數(shù)據(jù)去處理肋殴，今天偶然在群里看到，想用isform做一下試試坦弟。我非常喜歡gdc-client這個(gè)官方下載工具,它的數(shù)據(jù)以樣本的方式組織护锤，優(yōu)點(diǎn)是可以隨便組和，缺點(diǎn)是需要后期批量讀取和整理酿傍，代碼難度大烙懦。無(wú)所謂啦，習(xí)慣了都一樣赤炒。

日常小記：我8月22號(hào)講完了線上課氯析，趁現(xiàn)在疫情已經(jīng)穩(wěn)定了亏较，我趕緊一張機(jī)票把自己鏟回了老家。目前處于爸媽愛(ài)不釋手期（應(yīng)該是過(guò)幾天才會(huì)開始被嫌棄）掩缓。這一篇推文我寫了不到一個(gè)半小時(shí)雪情，被投喂了四次，一次是二姨家的玉米你辣，一次是鄰居家的大桃子巡通，一次是家門口中的葡萄，就在寫這段話的時(shí)候麻麻又去商店拎回來(lái)一箱牛奶哈哈哈哈

正文開始舍哄！

1.從網(wǎng)頁(yè)選擇數(shù)據(jù)宴凉，下載manifest文件

參考(本來(lái)這里應(yīng)該有個(gè)鏈接，但是被封了表悬，沒(méi)了弥锄，去國(guó)民聊天軟件搜生信星球看吧)，選擇file的時(shí)候找到miRNAseq的isform數(shù)據(jù)即可签孔。

2.使用gdc-client工具下載

注意：

將gdc-client(mac)或gdc-client.exe(windows)放在工作目錄下叉讥；

將manifest文件放在工作目錄下。

library(stringr)
if(!dir.exists("expdata"))dir.create("expdata")
dir()
#> [1] "1_gdc-client.html"            
#> [2] "1_gdc-client.Rmd"             
#> [3] "clinical"                     
#> [4] "expdata"                      
#> [5] "gdc-client.exe"               
#> [6] "gdc_manifest.2021-08-22.txt"  
#> [7] "metadata.cart.2021-08-24.json"
#> [8] "tcga_1.Rproj"
## 這句下載的代碼要在terminal運(yùn)行
#./gdc-client.exe download -m gdc_manifest.2021-08-22.txt -d expdata

length(dir("./expdata/"))
#> [1] 45

可以看到饥追，下載的文件是按樣本存放的图仓，我們需要得到的是表格，需要將他們批量讀入R語(yǔ)言并整理但绕。

3.整理表達(dá)矩陣

探索數(shù)據(jù)：先任選兩個(gè)counts文件讀取救崔，并觀察geneid的順序是否一致。

library(tidyverse)
x = read.table("expdata/0cbb843b-c6d4-4ad8-8678-6df63074d1a0/40a39e83-23c6-4240-a1a0-ddf221f23103.mirbase21.isoforms.quantification.txt",sep = "\t",header = T) 
head(x)
#>       miRNA_ID                isoform_coords
#> 1 hsa-let-7a-1 hg38:chr9:94175961-94175982:+
#> 2 hsa-let-7a-1 hg38:chr9:94175961-94175983:+
#> 3 hsa-let-7a-1 hg38:chr9:94175961-94175984:+
#> 4 hsa-let-7a-1 hg38:chr9:94175962-94175981:+
#> 5 hsa-let-7a-1 hg38:chr9:94175962-94175982:+
#> 6 hsa-let-7a-1 hg38:chr9:94175962-94175983:+
#>   read_count reads_per_million_miRNA_mapped
#> 1          5                       0.756391
#> 2          7                       1.058948
#> 3         12                       1.815339
#> 4        302                      45.686022
#> 5      12255                    1853.914558
#> 6      15830                    2394.734187
#>   cross.mapped        miRNA_region
#> 1            N mature,MIMAT0000062
#> 2            N mature,MIMAT0000062
#> 3            N mature,MIMAT0000062
#> 4            N mature,MIMAT0000062
#> 5            N mature,MIMAT0000062
#> 6            N mature,MIMAT0000062

可見(jiàn)這個(gè)數(shù)據(jù)里不只有count數(shù)據(jù)捏顺，還有前體和成熟體的id六孵，而且還是分了區(qū)間。接下來(lái)需要：

• 只要成熟體id和count兩列
• 按照成熟體分組求和幅骄，得出每個(gè)成熟體的count之和

一個(gè)文件整理成功后劫窒，對(duì)所有的樣本文件批量做相同的處理

x2 = x %>%
  dplyr::select(c(6,3)) %>% 
  group_by(miRNA_region) %>% 
  summarise(read_count = sum(read_count))

count_files = dir("expdata/",pattern = "*isoforms.quantification.txt$",recursive = T)

exp = list()
for(i in 1:length(count_files)){
  exp[[i]] <- read.table(paste0("expdata/",count_files[[i]]),sep = "\t",header = T) %>%
  dplyr::select(c(6,3)) %>% 
  group_by(miRNA_region) %>% 
  summarise(read_count = sum(read_count)) 
}
sapply(exp, nrow)
#>  [1] 650 828 799 791 805 725 778 783 751 781 703 780
#> [13] 757 767 731 800 714 781 722 749 777 795 735 828
#> [25] 832 760 805 911 629 719 765 733 801 825 839 924
#> [37] 730 815 710 740 776 788 711 908 832

檢查發(fā)現(xiàn)，每個(gè)文件都進(jìn)來(lái)的數(shù)據(jù)拆座，它的行數(shù)（也就是miRNA成熟體的個(gè)數(shù)）是不一樣多的主巍。

所以就不能cbind了，Resuce配merge挪凑，空著的地方會(huì)填充上NA孕索，給他改成０即可。

另躏碳，最后三行需要去掉搞旭，不是表達(dá)矩陣矩陣正式內(nèi)容。

m = Reduce(function(x, y) merge(x, y, by= 'miRNA_region',all = T), exp)
m[is.na(m)]=0
exp <- column_to_rownames(m,var = "miRNA_region")
exp = exp[-((nrow(exp)-2):nrow(exp)),]
dim(exp)
#> [1] 1753   45
exp[1:4,1:4]
#>                     read_count.x read_count.y
#> mature,MIMAT0000062        27647       167773
#> mature,MIMAT0000063        15016       122311
#> mature,MIMAT0000064         1196        10586
#> mature,MIMAT0000065          194         1521
#>                     read_count.x.1 read_count.y.1
#> mature,MIMAT0000062         183894         152801
#> mature,MIMAT0000063         206980         106415
#> mature,MIMAT0000064          13853          10759
#> mature,MIMAT0000065            912            978

4.行名轉(zhuǎn)換

行名里的前綴＂mature＂可以去掉。MIM開頭的是mirBase數(shù)據(jù)庫(kù)里的id肄渗，需要轉(zhuǎn)換為大多以hsa-miR開頭的成熟體id镇眷。

GDC數(shù)據(jù)庫(kù)使用的mirBasev21版本的id，我們轉(zhuǎn)換時(shí)需要使用相同的版本,使用miRBaseVersions.db包更絲滑

table(str_detect(rownames(exp),"mature,"))
#> 
#> TRUE 
#> 1753
rownames(exp) = str_remove(rownames(exp),"mature,")
library(miRBaseVersions.db)
mh <- select(miRBaseVersions.db,
               keys = rownames(exp),
               keytype = "MIMAT",
               columns = c("ACCESSION","NAME","VERSION"))
head(mh)
#>      ACCESSION          NAME VERSION
#> 1 MIMAT0000062 hsa-let-7a-5p      22
#> 2 MIMAT0000063 hsa-let-7b-5p      22
#> 3 MIMAT0000064 hsa-let-7c-5p      22
#> 4 MIMAT0000065 hsa-let-7d-5p      22
#> 5 MIMAT0000066 hsa-let-7e-5p      22
#> 6 MIMAT0000067 hsa-let-7f-5p      22
mh = mh[mh$VERSION=="21",]
mh = mh[match(rownames(exp),mh$ACCESSION),]
identical(rownames(exp),mh$ACCESSION)
#> [1] TRUE
table(!duplicated(mh$NAME))
#> 
#> TRUE 
#> 1753
rownames(exp) = mh$NAME
exp[1:4,1:4]
#>               read_count.x read_count.y
#> hsa-let-7a-5p        27647       167773
#> hsa-let-7b-5p        15016       122311
#> hsa-let-7c-5p         1196        10586
#> hsa-let-7d-5p          194         1521
#>               read_count.x.1 read_count.y.1
#> hsa-let-7a-5p         183894         152801
#> hsa-let-7b-5p         206980         106415
#> hsa-let-7c-5p          13853          10759
#> hsa-let-7d-5p            912            978

經(jīng)過(guò)這一圈折騰恳啥，終于得到一個(gè)正常行名的表達(dá)矩陣?yán)财印Ｎ冶容^了一下成熟體和前體id的差別和聯(lián)系，如下：

x = read.table("expdata/0cbb843b-c6d4-4ad8-8678-6df63074d1a0/40a39e83-23c6-4240-a1a0-ddf221f23103.mirbase21.isoforms.quantification.txt",sep = "\t",header = T) %>%
  dplyr::select(c(6,1)) %>% 
  distinct()
x$miRNA_region = str_remove(x$miRNA_region,"mature,")
x = merge(x,mh,by.x = "miRNA_region",by.y = "ACCESSION")  
head(x)
#>   miRNA_region     miRNA_ID          NAME VERSION
#> 1 MIMAT0000062 hsa-let-7a-3 hsa-let-7a-5p      21
#> 2 MIMAT0000062 hsa-let-7a-2 hsa-let-7a-5p      21
#> 3 MIMAT0000062 hsa-let-7a-1 hsa-let-7a-5p      21
#> 4 MIMAT0000063   hsa-let-7b hsa-let-7b-5p      21
#> 5 MIMAT0000064   hsa-let-7c hsa-let-7c-5p      21
#> 6 MIMAT0000065   hsa-let-7d hsa-let-7d-5p      21

5.列名轉(zhuǎn)換

這個(gè)仍然是參考之前的鏈接（http://www.reibang.com/p/af9b257c8a20）钝的，難得啊沒(méi)給我封了翁垂。

meta <- jsonlite::fromJSON("metadata.cart.2021-08-24.json")
ID = sapply(meta$associated_entities,
            function(x){x$entity_submitter_id})
file2id = data.frame(file_name = meta$file_name,
                     ID = ID)
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
exp[1:4,1:4]
#>               TCGA-W5-AA2O-01A-11R-A41D-13
#> hsa-let-7a-5p                        27647
#> hsa-let-7b-5p                        15016
#> hsa-let-7c-5p                         1196
#> hsa-let-7d-5p                          194
#>               TCGA-4G-AAZT-01A-11R-A41D-13
#> hsa-let-7a-5p                       167773
#> hsa-let-7b-5p                       122311
#> hsa-let-7c-5p                        10586
#> hsa-let-7d-5p                         1521
#>               TCGA-ZH-A8Y1-01A-11R-A41D-13
#> hsa-let-7a-5p                       183894
#> hsa-let-7b-5p                       206980
#> hsa-let-7c-5p                        13853
#> hsa-let-7d-5p                          912
#>               TCGA-YR-A95A-01A-12R-A41D-13
#> hsa-let-7a-5p                       152801
#> hsa-let-7b-5p                       106415
#> hsa-let-7c-5p                        10759
#> hsa-let-7d-5p                          978

6.表達(dá)矩陣的過(guò)濾

表達(dá)矩陣整理完成,需要過(guò)濾一下那些在很多樣本里表達(dá)量都為0的基因。過(guò)濾標(biāo)準(zhǔn)不唯一硝桩。

dim(exp)
#> [1] 1753   45
#exp = exp[rowSums(exp)>0,]
k = apply(exp, 1, function(x) sum(x > 0) > 9);table(k)
#> k
#> FALSE  TRUE 
#>   802   951
exp = exp[k,]
dim(exp)
#> [1] 951  45
exp[1:4,1:4]
#>               TCGA-W5-AA2O-01A-11R-A41D-13
#> hsa-let-7a-5p                        27647
#> hsa-let-7b-5p                        15016
#> hsa-let-7c-5p                         1196
#> hsa-let-7d-5p                          194
#>               TCGA-4G-AAZT-01A-11R-A41D-13
#> hsa-let-7a-5p                       167773
#> hsa-let-7b-5p                       122311
#> hsa-let-7c-5p                        10586
#> hsa-let-7d-5p                         1521
#>               TCGA-ZH-A8Y1-01A-11R-A41D-13
#> hsa-let-7a-5p                       183894
#> hsa-let-7b-5p                       206980
#> hsa-let-7c-5p                        13853
#> hsa-let-7d-5p                          912
#>               TCGA-YR-A95A-01A-12R-A41D-13
#> hsa-let-7a-5p                       152801
#> hsa-let-7b-5p                       106415
#> hsa-let-7c-5p                        10759
#> hsa-let-7d-5p                          978
exp = as.matrix(exp)

大功告成咯~