基于简易小工具下载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)
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))
ind_gene <- which(rownames(z_rna) == "hsa-mir-377")
table(event_rna[ind_gene,]) #一组422个，一组123个
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'))
其实我觉得我的结果不是很好，上下曲线方基本重合。如果有大神能指出存在的问题，感激不尽！
