听说你想要TCGA的miRNA表达数据

作者: 小洁忘了怎么分身 | 来源:发表于2021-08-24 15:04 被阅读0次

#本月回顾#基因组分析方法教程汇总
听说你想要TCGA的miRNA表达数据
#TCGA系列#TCGA基因/miRNA表达谱数据整合
#TCGA系列#TCGA基因/miRNA表达谱数据整合(二)
零代码合并TCGA表达谱数据
TCGA找差异表达
【R语言】获取最新的人的所有miRNA的ID号
2020-01-07-TCGA之RTCGA
计算MiRNA–mRNA表达相关性
一文解决如何下载前体或者成熟体miRNA表达矩阵（TCGA数据库

以前都是直接拿前体数据去处理，今天偶然在群里看到，想用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)

大功告成咯~

#本月回顾#基因组分析方法教程汇总
TCGA数据分析: TCGA基因/miRNA表达谱及临床数据下载TCGA基因/miRNA表达谱数据整合TCGA基因...
听说你想要TCGA的miRNA表达数据
以前都是直接拿前体数据去处理，今天偶然在群里看到，想用isform做一下试试。我非常喜欢gdc-client这个官...
#TCGA系列#TCGA基因/miRNA表达谱数据整合
上期(#TCGA系列#TCGA基因/miRNA表达谱及临床数据下载)介绍了使用TCGA 的API下载肿瘤表达谱及...
#TCGA系列#TCGA基因/miRNA表达谱数据整合(二)
上期(#TCGA系列#TCGA基因/miRNA表达谱数据整合)使用shell 对多样本表达谱文件整合,实现方式是...
零代码合并TCGA表达谱数据
前面给大家介绍了☞如何使用R代码来合并TCGA数据库下载得到的RNAseq和miRNA seq的表达谱数据☜[ht...
TCGA找差异表达
firebrowse下载TCGA数据 ##用R处理数据 dat = read.table("mirna.txt",...
【R语言】获取最新的人的所有miRNA的ID号
前面我们在介绍TCGA数据库数据挖掘的时候，课程中使用了人了所有miRNA的ID号。 TCGA数据库介绍及数据挖掘...
2020-01-07-TCGA之RTCGA
TCGA数据的另一个重要的包：RTCGA参考学习资料：TCGA之miRNA预测临床表现RTCGA（1）数据概况与...
计算MiRNA–mRNA表达相关性
计算MiRNA–mRNA表达相关性看到一篇文章提到了这个分析，作者下载了TCGA数据库的OV的Affymetri...
一文解决如何下载前体或者成熟体miRNA表达矩阵（TCGA数据库
本文目的：解决纯R语言代码版本下载最新前体或者成熟体miRNA表达矩阵（TCGA数据库）解决零代码如何下载最新...

听说你想要TCGA的miRNA表达数据

正文开始！

1.从网页选择数据，下载manifest文件

2.使用gdc-client工具下载

3.整理表达矩阵

4.行名转换

5.列名转换

6.表达矩阵的过滤

相关文章

#本月回顾#基因组分析方法教程汇总

听说你想要TCGA的miRNA表达数据

#TCGA系列#TCGA基因/miRNA表达谱数据整合

#TCGA系列#TCGA基因/miRNA表达谱数据整合(二)

零代码合并TCGA表达谱数据

TCGA找差异表达

【R语言】获取最新的人的所有miRNA的ID号

2020-01-07-TCGA之RTCGA

计算MiRNA–mRNA表达相关性

一文解决如何下载前体或者成熟体miRNA表达矩阵（TCGA数据库

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

生信