前言：這個(gè)帖子本應(yīng)該在2018年就要寫完的誊辉，有個(gè)網(wǎng)友甚至當(dāng)時(shí)加了我微信給了50塊錢我都沒(méi)有給他一個(gè)代碼版职辅。在這里希望這個(gè)網(wǎng)友能看到传透，對(duì)得起你的50塊担败。（內(nèi)心充滿愧疚刮便，一工作后就被工作的各種事情占滿了精力然后每天下班回家都只想在床上安靜的躺著什么都不想干）
以前還會(huì)使用TCGA官網(wǎng)工具gdc-client去下載數(shù)據(jù)脯爪，自從后來(lái)發(fā)現(xiàn)可以用TCGAbiolinks這個(gè)R包之后启上，尤其是與后續(xù)的分析可以在R中實(shí)現(xiàn)無(wú)縫銜接碍侦，果斷放棄了前面的工具粱坤。下面就用這個(gè)包來(lái)寫一個(gè)代碼版的數(shù)據(jù)下載吧。
廢話有點(diǎn)多還啰嗦瓷产，不說(shuō)了站玄，切入正題，以下載肺腺癌的miRNA成熟體為主示例濒旦。
# Author: ScorpioGirl
# Date: 2019-12-07
# Aim: Download miRNA isoform expression from TCGA database</pre>
加載所需要的包以及環(huán)境設(shè)置
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>
數(shù)據(jù)查詢與下載
# 可以通過(guò)TCGAbiolinks:::getGDCprojects()$project_id得到每個(gè)項(xiàng)目的ID
TCGAbiolinks:::getGDCprojects()$project_id
project_id <- "TCGA-LUAD"
?
# 可以通過(guò)TCGAbiolinks:::getProjectSummary(project_id)得到可以下載的數(shù)據(jù)類型
TCGAbiolinks:::getProjectSummary(project_id)
data_category <- "Transcriptome Profiling"
?
# data.type
data_type <- "Isoform Expression Quantification"
?
# 查詢數(shù)據(jù)
query <- GDCquery(project = project_id,
data.category = data_category,
data.type = data_type)
?
# 下載
GDCdownload(query)
?
# 整合株旷，這里得到的dataAssy是一個(gè)數(shù)據(jù)框，并且每一個(gè)樣本的每一個(gè)miRNA一行尔邓，后面需要處理一下
dataAssy <- GDCprepare(query = query,
summarizedExperiment=F)</pre>
下載的數(shù)據(jù)如下圖晾剖，第一列根據(jù)命名習(xí)慣我們可以知道是前體miRNA
得到有多少個(gè)樣本
總共有567個(gè)樣本，并且每個(gè)病人中miRNA表達(dá) 的數(shù)目不一樣
sample_id <- as.data.frame(table(dataAssy$barcode))</pre>
得到有多少個(gè)成熟體miRNA
總共有2216行梯嗽，其中有三個(gè)是precursor齿尽，stemloop以及沒(méi)有注釋上的unannotated
miRNA_id <- as.data.frame(table(dataAssy$miRNA_region))</pre>
生成一個(gè)有567個(gè)樣本和2216個(gè)成熟體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>
處理數(shù)據(jù)得到RPM值
未表達(dá)的miRNA這里處理成了0，后續(xù)根據(jù)需要可以進(jìn)行過(guò)濾灯节，如去掉在80%的樣本中表達(dá)值都為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>
得到的表達(dá)譜如下：
對(duì)名字進(jìn)行處理
#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進(jìn)行ID轉(zhuǎn)換
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>
最后的成熟體表達(dá)譜
數(shù)據(jù)是這個(gè)樣子的循头，mir變成了miR
愉快的保存
# 我喜歡RData格式绵估，省空間又方便下次使用加載。
save(miRNA_matrue_RPM2,file = paste(project_id,"_miRNA_matrue_RPM.RData",sep=""))
# 當(dāng)然也可以是你們習(xí)慣的格式
write.table(miRNA_matrue_RPM2,file = paste(project_id,"_miRNA_matrue_RPM.xls",sep=""),sep="\t",row.names = F,quote=F)</pre>
４丁Ｒ际俊！劃重點(diǎn)總結(jié)３ゾｕ锞取！
1.一個(gè)前體可以生成兩個(gè)成熟體或者其他的
如：hsa-mir-10a螟蒸，有兩個(gè)成熟體MIMAT0000253盒使，MIMAT0004555和precursor。
2.一個(gè)成熟體可以由多個(gè)前體生成
如：MIMAT0000062由hsa-let-7a-1七嫌，hsa-let-7a-2和hsa-let-7a-3三個(gè)前體生成少办。
3.ID轉(zhuǎn)換的時(shí)候，慎重選擇miRBase數(shù)據(jù)庫(kù)版本
如：VERSION ==22.0诵原，最后會(huì)發(fā)現(xiàn)只剩下了2197個(gè)ID英妓，如果選擇21版本，發(fā)現(xiàn)還是2213個(gè)ID绍赛。Because：GDC Data Portal 中的miRNA 數(shù)據(jù)是基于 GRCh38 和 mirbase21的蔓纠。也不能冒然使用最新版的數(shù)
