miRNA

TCGA-miRNA的生存分析

2019-12-16  本文已影响0人  Minka__

引言:之前写了mRNA的生存分析,对于miRNA来说基本一致,只不过因为从TCGA上下载的miRNA的reads文件和mRNA有一点不一样,需要额外处理一下。
TCGA-mRNA的生存分析

TCGA上直接下载的miRNA的表达量数据如图所示:


miRNA.png

此处采用RPM值,处理结果如下


miRNA_exp.png
library(SummarizedExperiment)
library(TCGAbiolinks)
library(survival)
library(survminer)
library(maftools)
TCGAbiolinks:::getProjectSummary("TCGA-LIHC") 
clinical <- GDCquery_clinic(project = "TCGA-LIHC", type = "clinical")
query <- GDCquery(project = "TCGA-LIHC", 
                  experimental.strategy = "miRNA-Seq",
                  data.category = "Transcriptome Profiling", 
                  data.type = "miRNA Expression Quantification"
                  )
GDCdownload(query)
LIHC_miseq <- GDCprepare(query)
rownames(LIHC_miseq)<-LIHC_miseq$miRNA_ID
rpm_miR<-colnames(LIHC_miseq)[grep('reads_per_million',colnames(LIHC_miseq))]
miR<-LIHC_miseq[,rpm_miR]
colnames(miR) <- gsub("reads_per_million_miRNA_mapped_","", colnames(miR))
miR_matrix<-as.matrix(miR)
#-------------这里第一部分结束----------下载完miRNA表达数据以及clinical数据
#下面内容和mRNA的生存分析内容一样
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(miR_matrix),typesample = c("TP"))
gene_exp <- miR_matrix[c("hsa-mir-4745"),samplesTP]
names(gene_exp) <-  sapply(strsplit(names(gene_exp),'-'),function(x) paste(x[1:3],collapse="-"))
clinical$GENE <- gene_exp[clinical$submitter_id]
#-----------------第二部分结束------整合完最终需要的文件-------------
df<-subset(clinical,select=c(submitter_id,vital_status,days_to_death,days_to_last_follow_up,GENE))
df$os<-ifelse(df$vital_status=='Alive',df$days_to_last_follow_up,df$days_to_death)
#alive的样本采用days_to_last_follow_up
df <- df[!is.na(df$GENE),]#去掉表达量为0的样本
df$exp=''
df[df$GENE >= mean(df$GENE),]$exp <- "High"
df[df$GENE <  mean(df$GENE),]$exp <- "Low"
df[df$vital_status=='Dead',]$vital_status <- 2
df[df$vital_status=='Alive',]$vital_status <- 1
df$vital_status <- as.numeric(df$vital_status)
#--------------第三部分结束----对miRNA的高低表达量进行确定------
fit <- survfit(Surv(os, vital_status)~exp, data=df) # 根据表达建模
ggsurvplot(fit,pval=TRUE)


Rplot01.png
a<-ggsurvplot(fit,legend.title = "Expression",palette = c('red','black'), 
           pval = TRUE,
           risk.table = TRUE,
           tables.height = 0.2,
           tables.theme = theme_cleantable(),
           xlab='Time(days)'
)
ggsave('/home/zhang/xxxx/',print(a)) #输出图片
1.png
上一篇下一篇

猜你喜欢

热点阅读