TCGA miRNASeq 数据生存分析

基于简易小工具下载TCGA表达数据做生存分析。软件工具:R。参考教程:https://www.biostars.org/p/153013/。有建议或问题希望指出,互相促进互相学习:)


第一步,在简易小工具上下载TCGA miRNASeq 和临床数据


1、由于教程中是在FireBrowse ( http://gdac.broadinstitute.org/ )中下载的KIRC(肾透明细胞癌)数据,所以我在这里也下载了KIRC数据。

下载miRNASeq时,将文本下载到一个指定文件夹后,点击“合并文件”,我们会在指定文件夹中找到合并后的文本文件Merge_matrix.txt。

2、合并完后的数据列为geneID,行为Barcode(TCGA对样本的分类),但是,后面我做数据时发现中间有一行数据是无用的,显示read_count(检查用的COAD RNASeq数据也出现了一样的情况),以前的简易小工具教程中没有提到过,可能是在合并文件时有点小瑕疵。

     下面是我找此行数据的代码,然后我直接在原数据中删除了那一行。

a<-as.numeric(as.matrix(rna)

which(is.na(a))

3、接下来下载临床数据,同RNASeq 数据一样下载到一指定文件夹,点击“Clinical”,会得到文件Clinical_matrix.txt。


第二步利用R处理文件数据,做RNASeq数据的生存分析


1、首先需要安装相应的R包:

需要的第一个包肯定是“survival”,它里面的Surv函数能直接对数据进行生存分析。

library(survival)

我们还需要 “limma”包,需要用到里面的voon函数,对数据进行标准化。

source("http://bioconductor.org/biocLite.R")

biocLite("limma")

library(limma)

2、处理RNASeq文件数据:

1)导入文件(记得改变工作目录setwd(目录)):

rna <- read.table('miRNA/Merge_matrix.txt',nrows=60489,header=T,row.names=1,sep='\t') 

2)数据中有许多“0”数据,剔除超过50%的表达值为0的样本:

rna <- read.table('RNA/Merge_matrix.txt',nrows=60490, header=T,row.names=1,sep='\t')
rem <- function(x){
   x <- as.matrix(x)
   x <- t(apply(x,1,as.numeric))
   r <-as.numeric(apply(x,1,function(i) sum(i == 0)))
   remove <- which(r >dim(x)[2]*0.5)
   return(remove)
}
remove<- rem(rna)
rna<- rna[-remove,]

3)区别正常和肿瘤样本:根据TCGA样本分类的原则,第4个参数指样本类型,“Tumor types range from 01 - 09, normal types from 10 - 19 and control samples from 20 - 29”例如TCGA-CM-4746-01,第4个参数是01,所以是肿瘤样本。可以看出第4个参数第1位即总第14位如果是“0”,则为肿瘤样本,“1”则为正常样本。

table(substr(colnames(rna),14,14))  #得到545个肿瘤样本,71个正常样本
n_index <- which(substr(colnames(rna),14,14) == '1')
t_index <- which(substr(colnames(rna),14,14) == '0')

4)对数据进行标准化:

voom(): Transform count data to log2-counts per million, estimate the mean-variance relationship and use this to compute appropriate observational-level weights. The data are then ready for linear modelling.

vm <- function(x){
  cond <- factor(ifelse(seq(1,dim(x)[2],1) %in% t_index, 1, 0))
  d <- model.matrix(~1+cond)
  x <- t(apply(x,1,as.numeric))
  ex <- voom(x,d,plot=F)
  return(ex$E)
}
rna_vm <- vm(rna)
colnames(rna_vm) <- gsub('\\.','-',substr(colnames(rna),1,12))
hist(rna_vm)  #检查数据是否正常分布

5)处理数据: 

将数据转化为z-score(标准分数),算出每个病人每个基因表达值距离平均数多少个标准差,如果z>+1.96(大约p=0.05或标准差为2)认为出现不同表达

公式:z = [(value gene X in tumor Y)-(mean gene X in normal)]/(standard deviation X in normal)

scal <- function(x,y){
  mean_n <- rowMeans(y)  # 正常样本表达量均值
  sd_n <- apply(y,1,sd)  # 正常样本表达量标准差
  res <- matrix(nrow=nrow(x),ncol=ncol(x))
  colnames(res) <- colnames(x)
  rownames(res) <- rownames(x)
  for(i in 1:dim(x)[1]){
    for(j in 1:dim(x)[2]){
      res[i,j] <-(x[i,j]-mean_n[i])/sd_n[i]
    }
 }
  return(res)
}
z_rna <- scal(rna_vm[,t_index],rna_vm[,n_index])
rownames(z_rna) <- sapply(rownames(z_rna), function(x) unlist(strsplit(x,'\\|'))[[1]]) #用基因名设置行名
rm(rna_vm)
event_rna <- t(apply(z_rna, 1, function(x) ifelse(abs(x) >1.96,1,0)))


3、处理临床数据:

1)导入数据:

clinical <-t(read.table('Clinical/Clinical_matrix.txt',nrows=538,row.names=1,header=T,sep='\t'))

2)匹配z-rna 和clinical geneID:

ID=clinical[0,]
sum(colnames(ID) %in% colnames(z_rna)) #得到516个病人数据

Surv用法:Surv(time,event),需要生存时间和状态结果

3)判断状态的向量:

status= clinical[12,]
table(status)  #Alive  360   Dead 177
event <- ifelse(status == 'alive', 0,1)

4)时间数据在表中第14列:

time= as.numeric(clinical[13,])  
5)因为需要在两组数据中相同数量的患者来配对,所以对两组数据进行匹配:
ind_tum <- which(unique(colnames(z_rna)) %in% colnames(ID))
ind_clin <- which(colnames(ID) %in% colnames(z_rna))

6)生存分析需要选任一基因进行分组:

ind_gene <- which(rownames(z_rna) == "hsa-mir-377")
table(event_rna[ind_gene,])  #一组422个,一组123个

4、生存分析:

s <- survfit(Surv(time[ind_clin],event[ind_clin])~event_rna[ind_gene,ind_tum])
s1 <-tryCatch(survdiff(Surv(time[ind_clin],event[ind_clin])~event_rna[ind_gene,ind_tum]),error = function(e) return(NA))
pv <- round(1-pchisq(s1$chisq,length(s1$n)-1),3)[[1]]

生存分析曲线:

plot(survfit(Surv(as.numeric(as.character(time))[ind_clin],event[ind_clin])~event_rna[ind_gene,ind_tum]),col=c(1:3),frame=F,lwd=2,main=paste('KIRK',rownames(z_rna)[ind_gene],sep='\n'))
x1 <- ifese ( is.na(as.numeric(summary(s)$table[,'median'][1])),'NA',as.numeric(summary(s)$table[,'median'][1]))
x2 <- as.numeric(summary(s)$table[,'median'][2])
if(x1 != 'NA' && x2 != 'NA'){ lines(c(0,x1),c(0.5,0.5),col='blue') lines(c(x1,x1),c(0,0.5),col='black') lines(c(x2,x2),c(0,0.5),col='red')}
legend(1800,0.995,legend=paste('p.value =',pv[[1]],sep=''),bty='n',cex=1.4)
legend(max(time[ind_clin],na.rm = T)*0.7,0.94,) legend=c(paste('NotAltered=',x1),paste('Altered=',x2)),bty='n',cex=1.3,lwd=3,col=c('black','red'))

其实我觉得我的结果不是很好,上下曲线方基本重合。如果有大神能指出存在的问题,感激不尽!




  • 发表于 2017-05-22 05:46
  • 阅读 ( 16558 )
  • 分类:方案研究

18 条评论

请先 登录 后评论
不写代码的码农
金晓妍

3 篇文章

作家榜 »

  1. 祝让飞 118 文章
  2. 柚子 91 文章
  3. 刘永鑫 64 文章
  4. admin 57 文章
  5. 生信分析流 55 文章
  6. SXR 44 文章
  7. 张海伦 31 文章
  8. 爽儿 25 文章