TCGA miRNASeq 数据生存分析

第一步，在简易小工具上下载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'))

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