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包
代码语言：javascript复制rm(list = ls())
options(stringsAsFactors = F)
# 作者：DoubleHelix，微信公众号:MedBioInfoCloud
library(rjson) #没有安装事先安装
library(dplyr) #没有安装事先安装
library(tidyr)
library(limma) #没有安装事先安装
library(stringr) #没有安装事先安装
library(miRBaseVersions.db)处理一下json文件。
代码语言：javascript复制###---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数据下载与整理】中一样的。
这里读入其中一个文件：
代码语言：javascript复制filepath <- dir(path ="./raw",pattern=".isoforms.quantification.txt$",full.names = TRUE,recursive = T)
oneSampExp <- read.table(filepath[1],header = T,check.names = F)
选择我们需要的列，并重新命名。
代码语言：javascript复制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多个。
代码语言：javascript复制acc <- oneSampExp$ACCESSION
acc <- acc[!duplicated(acc)]
length(acc)
我们可以通过miRBaseVersions.db包进行转换。
代码语言：javascript复制ano <- select(miRBaseVersions.db,
keys = acc,
keytype = "VW-MIMAT-22.0",
columns = c("ACCESSION","NAME"))
head(ano)
与前面的数据融合。
代码语言：javascript复制data <- merge(ano,oneSampExp,by= "ACCESSION")
head(data)
我们只提取hsa-let-7a-1的RPM数据。求和是3670。
代码语言：javascript复制temp <- data[data$miRNA_ID=="hsa-let-7a-1",]
sum(temp$RPM)
这里只是一个病人的数据，该病人的ID是：TCGA-AA-3531-01A-01T-0822-13
代码语言：javascript复制tempPath <- unlist(strsplit(filepath[1],"/"))
filename <- tempPath[length(tempPath)]
filesNameToBarcode[filename,2]
利用这个ID，我们加载前面通过文章【TCGA数据库：miRNA数据下载与整理】得到的数据。提取该病人的RPM数据。当然，个别miRNA可能会有点出入，值没多大区别。
代码语言：javascript复制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数据下载与整理】介绍的差不多，那样我们会得到如下这样的数据。
代码在文章【TCGA数据库：miRNA数据下载与整理】，以前付费过的，你回复文章中的关键词，重新获取。
