听说你想要TCGA的miRNA表达数据
2021-08-24 本文已影响0人
小洁忘了怎么分身
以前都是直接拿前体数据去处理，今天偶然在群里看到，想用isform做一下试试。我非常喜欢gdc-client这个官方下载工具,它的数据以样本的方式组织，优点是可以随便组和，缺点是需要后期批量读取和整理，代码难度大。无所谓啦，习惯了都一样。
日常小记：我8月22号讲完了线上课，趁现在疫情已经稳定了，我赶紧一张机票把自己铲回了老家。目前处于爸妈爱不释手期（应该是过几天才会开始被嫌弃）。这一篇推文我写了不到一个半小时，被投喂了四次，一次是二姨家的玉米，一次是邻居家的大桃子，一次是家门口中的葡萄，就在写这段话的时候麻麻又去商店拎回来一箱牛奶哈哈哈哈
正文开始！
1.从网页选择数据，下载manifest文件
参考(本来这里应该有个链接，但是被封了，没了，去国民聊天软件搜生信星球看吧)，选择file的时候找到miRNAseq的isform数据即可。
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运行
#./gdc-client.exe download -m gdc_manifest.2021-08-22.txt -d expdata
length(dir("./expdata/"))
#> [1] 45
可以看到，下载的文件是按样本存放的，我们需要得到的是表格，需要将他们批量读入R语言并整理。
3.整理表达矩阵
探索数据：先任选两个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
可见这个数据里不只有count数据，还有前体和成熟体的id，而且还是分了区间。接下来需要：
- 只要成熟体id和count两列
- 按照成熟体分组求和，得出每个成熟体的count之和
一个文件整理成功后，对所有的样本文件批量做相同的处理
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
检查发现，每个文件都进来的数据，它的行数（也就是miRNA成熟体的个数）是不一样多的。
所以就不能cbind了，Resuce配merge，空着的地方会填充上NA，给他改成０即可。
另，最后三行需要去掉，不是表达矩阵矩阵正式内容。
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.行名转换
行名里的前缀＂mature＂可以去掉。MIM开头的是mirBase数据库里的id，需要转换为大多以hsa-miR开头的成熟体id。
GDC数据库使用的mirBasev21版本的id，我们转换时需要使用相同的版本,使用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
经过这一圈折腾，终于得到一个正常行名的表达矩阵啦。我比较了一下成熟体和前体id的差别和联系，如下：
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.列名转换
这个仍然是参考之前的链接（https://www.jianshu.com/p/af9b257c8a20），难得啊没给我封了。
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.表达矩阵的过滤
表达矩阵整理完成,需要过滤一下那些在很多样本里表达量都为0的基因。过滤标准不唯一。
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)
大功告成咯~
