TCGA数据库：miRNA数据下载与整理(2)

DoubleHelix

发布于 2021-03-16 16:00:26

5.7K0

发布于 2021-03-16 16:00:26

文章被收录于专栏：生物信息云生物信息云

我前面有一篇介绍TCGA数据库中miRN

miRNA Expression Quantification数据。

我整理得到的数据结果是这样的。

这里我们可以发现，miRNA的前体可能对应多个成熟的miRNA，比如hsa-let-7a-1，有两个对应的成熟体，MIMAT0000062(hsa-let-7a-5p)和MIMAT0004481(hsa-let-7a-3p)。这里的值是对所有成熟体miRNA求和的结果。

如果我们想要得到成熟体miRNA的数据，我们需要下载Isoform Expression Quantification的数据。

下载后的数据处理方式和前面文章差不多。下载后的数据，我们打开单个文件看一下，比miRNA Expression Quantification数据多了2列。在miRNA_region列中就有是否成熟的信息。我们可以参考前面的文章【miRNA注释包:miRBaseVersions.db】进行转换为我们常用的名称，比如：hsa-let-7a-5p。

下面我们来看一下，这些数据的差异。

加载一些需要用到的R包

rm(list = ls())
options(stringsAsFactors = F)
# 作者：DoubleHelix，微信公众号:MedBioInfoCloud
library(rjson) #没有安装事先安装
library(dplyr) #没有安装事先安装
library(tidyr)
library(limma) #没有安装事先安装
library(stringr) #没有安装事先安装
library(miRBaseVersions.db)

处理一下json文件。

###---1.处理json文件----------------------
jsonFile <- fromJSON(file="./metadata.cart.2021-03-10.json")
filesNameToBarcode <- data.frame(filesName = c(),TCGA_Barcode = c())
for(i in 1:length(jsonFile)){
  TCGA_Barcode <- jsonFile[[i]][["associated_entities"]][[1]][["entity_submitter_id"]]
  file_name <- jsonFile[[i]][["file_name"]]
  filesNameToBarcode <- rbind(filesNameToBarcode,data.frame(filesName = file_name,TCGA_Barcode = TCGA_Barcode))
}
rownames(filesNameToBarcode) <- filesNameToBarcode[,1]

这个和前文【TCGA数据库：miRNA数据下载与整理】中一样的。

这里读入其中一个文件：

filepath <- dir(path ="./raw",pattern=".isoforms.quantification.txt$",full.names = TRUE,recursive = T)
oneSampExp <- read.table(filepath[1],header = T,check.names = F)

选择我们需要的列，并重新命名。

oneSampExp <- oneSampExp[grep("mature",oneSampExp$miRNA_region),]
oneSampExp <- separate(oneSampExp,into=c("State","ACCESSION"),col="miRNA_region",sep=",")
oneSampExp <- oneSampExp[c("miRNA_ID","ACCESSION","read_count","reads_per_million_miRNA_mapped")]
colnames(oneSampExp)<- c("miRNA_ID","ACCESSION","count","RPM")

在众多前体下，这里获得536个成熟体。据我自己所知，目前已知的miRNA成熟体应该有2000多，这里一个样本只检测到500多，是因为不是所有miRNA在一个样本中都能检测到，你把所有的样本加起来统计大概也就2000多个。

acc <- oneSampExp$ACCESSION
acc <- acc[!duplicated(acc)]
length(acc)

我们可以通过miRBaseVersions.db包进行转换。

ano <- select(miRBaseVersions.db, 
                 keys = acc, 
                 keytype = "VW-MIMAT-22.0", 
                 columns = c("ACCESSION","NAME"))
head(ano)

与前面的数据融合。

data <- merge(ano,oneSampExp,by= "ACCESSION")
head(data)

我们只提取hsa-let-7a-1的RPM数据。求和是3670。

temp <- data[data$miRNA_ID=="hsa-let-7a-1",]
sum(temp$RPM)

这里只是一个病人的数据，该病人的ID是：TCGA-AA-3531-01A-01T-0822-13

tempPath <- unlist(strsplit(filepath[1],"/"))
filename <- tempPath[length(tempPath)]
filesNameToBarcode[filename,2]

利用这个ID，我们加载前面通过文章【TCGA数据库：miRNA数据下载与整理】得到的数据。提取该病人的RPM数据。当然，个别miRNA可能会有点出入，值没多大区别。

load("miRNAdata.Rdata")
rpm <- miRNAdata[["miRNA_RPM_Expdata"]][["Exp"]]
rpm["hsa-let-7a-1",filesNameToBarcode[filename,2]]%>%as.numeric()

这里的值就是一样的。也就是说，我们通过miRNA Expression Quantification数据得到的是miRNA前体的数据。一个miRNA前体可能有多个成熟体，一个成熟体可能来自多个前体。

当然，还有其他处理方式，比如利用TCGAbiolinks包，关于TCGAbiolinks包我在很久前的一篇文章【TCGA数据挖掘（一）：TCGAbiolinks包介绍】中有介绍，同时，也有差异分析的教程【TCGAbiolinks包下载TCGA数据进行表达差异分析-乳腺癌案例】，但是现在这个包下载数据是真的感人，完全看GDC的心情，而且下载很慢，所以我是自己下载数据处理，我也发了自己下载提取数据进行差异分析的文章【一文就会TCGA数据库基因表达差异分析】，当然还有GDCRNATools包【TCGA数据库：GDCRNATools包下载数据、处理数据以及差异分析】。这个包下载数据的时候有时候也会挂，开VPN有时候可以解决。但我们可以自己下载数据后用该包的函数提取数据。如果你想偷懒，那么你可以在一些数据库中直接下载，比如：GDAC Firehose数据库，以及UCSC数据库【UCSC数据库下载TCGA数据需要注意的细节】，但从UCSC下载的数据和我们前面处理的一样，是miRNA Expression Quantification的数据，而不是Isoform Expression Quantification。当然还是TCGAbiolinks包好用。

最后说明，上面我只是读入一个样本的数据来进行讲解，你可以写一个循环，把所有文件遍历读入，处理合并就行，和前文【TCGA数据库：miRNA数据下载与整理】介绍的差不多，那样我们会得到如下这样的数据。