获取TCGA中miRNA成熟体表达

前言：这个帖子本应该在2018年就要写完的，有个网友甚至当时加了我微信给了50块钱我都没有给他一个代码版。在这里希望这个网友能看到，对得起你的50块。（内心充满愧疚，一工作后就被工作的各种事情占满了精力然后每天下班回家都只想在床上安静的躺着什么都不想干）

以前还会使用TCGA官网工具gdc-client去下载数据，自从后来发现可以用TCGAbiolinks这个R包之后，尤其是与后续的分析可以在R中实现无缝衔接，果断放弃了前面的工具。下面就用这个包来写一个代码版的数据下载吧。

废话有点多还啰嗦，不说了，切入正题，以下载肺腺癌的miRNA成熟体为主示例。

# Author: ScorpioGirl
# Date: 2019-12-07
# Aim: Download miRNA isoform expression from TCGA database</pre>

加载所需要的包以及环境设置

library(pacman)
p_load(TCGAbiolinks)
p_load(stringi)
packageVersion("TCGAbiolinks")

rm(list = ls())
options(stringsAsFactors = F)
setwd("/media/bmk/Elements SE/2.work/DZH/3.TCGA")</pre>

数据查询与下载

# 可以通过TCGAbiolinks:::getGDCprojects()$project_id得到每个项目的ID
TCGAbiolinks:::getGDCprojects()$project_id
project_id <- "TCGA-LUAD"

# 可以通过TCGAbiolinks:::getProjectSummary(project_id)得到可以下载的数据类型
TCGAbiolinks:::getProjectSummary(project_id)
data_category <- "Transcriptome Profiling"

# data.type
data_type <- "Isoform Expression Quantification"

# 查询数据
query <- GDCquery(project = project_id,
data.category = data_category,
data.type = data_type)

# 下载
GDCdownload(query)

# 整合，这里得到的dataAssy是一个数据框，并且每一个样本的每一个miRNA一行，后面需要处理一下
dataAssy <- GDCprepare(query = query,
summarizedExperiment=F)</pre>

下载的数据如下图，第一列根据命名习惯我们可以知道是前体miRNA

img

得到有多少个样本

总共有567个样本，并且每个病人中miRNA表达的数目不一样

sample_id <- as.data.frame(table(dataAssy$barcode))</pre>

img

得到有多少个成熟体miRNA

总共有2216行，其中有三个是precursor，stemloop以及没有注释上的unannotated

miRNA_id <- as.data.frame(table(dataAssy$miRNA_region))</pre>

img

生成一个有567个样本和2216个成熟体miRNA的矩阵

miRNA_matrue_RPM <- matrix(NA,ncol = nrow(sample_id),nrow = nrow(miRNA_id))
colnames(miRNA_matrue_RPM) <- sample_id$Var1
rownames(miRNA_matrue_RPM) <- as.character(miRNA_id$Var1)</pre>

处理数据得到RPM值

未表达的miRNA这里处理成了0，后续根据需要可以进行过滤，如去掉在80%的样本中表达值都为0的miRNA

for(i in 1:nrow(sample_id)){
 temp1 <- dataAssy[which(dataAssy$barcode==as.character(sample_id[i,1])),]

 for(j in 1:nrow(miRNA_id)){
     loc <- which(temp1$miRNA_region==as.character(miRNA_id[j,1])
     if(length(loc)>0){
       miRNA_matrue_RPM[j,i] <- sum(temp1[loc,4])
     }else{
           miRNA_matrue_RPM[j,i] <- 0
           }
     }

 print(i)
}

miRNA_matrue_RPM1 <- miRNA_matrue_RPM[1:2213,]</pre>

得到的表达谱如下：

img

对名字进行处理

#name <- unlist(strsplit(rownames(miRNA_matrue_RPM1),split="mature,",fixed=T))[seq(2,2213*2,by=2)]
name <- substr(rownames(miRNA_matrue_RPM1),8,nchar(rownames(miRNA_matrue_RPM1)))

# 使用R包miRBaseVersions.db进行ID转换
p_load(miRBaseVersions.db)
items <- select(miRBaseVersions.db,
keys = name,
keytype = "MIMAT",
columns = c("ACCESSION","NAME","VERSION"))
id_name <- items[items$VERSION == 21.0, c("ACCESSION","NAME")]
miRNA_matrue_RPM2 <- cbind(id_name,miRNA_matrue_RPM1)</pre>

最后的成熟体表达谱

数据是这个样子的，mir变成了miR

img

愉快的保存

# 我喜欢RData格式，省空间又方便下次使用加载。
save(miRNA_matrue_RPM2,file = paste(project_id,"_miRNA_matrue_RPM.RData",sep=""))

# 当然也可以是你们习惯的格式
write.table(miRNA_matrue_RPM2,file = paste(project_id,"_miRNA_matrue_RPM.xls",sep=""),sep="\t",row.names = F,quote=F)</pre>

！！！划重点总结！！！

1.一个前体可以生成两个成熟体或者其他的
如：hsa-mir-10a，有两个成熟体MIMAT0000253，MIMAT0004555和precursor。

2.一个成熟体可以由多个前体生成
如：MIMAT0000062由hsa-let-7a-1，hsa-let-7a-2和hsa-let-7a-3三个前体生成。

3.ID转换的时候，慎重选择miRBase数据库版本
如：VERSION ==22.0，最后会发现只剩下了2197个ID，如果选择21版本，发现还是2213个ID。Because：GDC Data Portal 中的miRNA 数据是基于 GRCh38 和 mirbase21的。也不能冒然使用最新版的数