前言

很多时候，我们分析完差异表达基因后，发现会得到一大堆差异基因，常见的做法有降低Pvalue的阈值，挑选fold-change最大的基因，做通路富集然后挑重要通路上的基因。本文主要讲述在差异基因表达分析后，用随机森林的方法进一步找出差异基因里面，对表型影响最大的基因，作为我们的下游研究。

本文主要参考下列文章

对方向扔了一个赤果果的生信分析SCI思路
机器学习之随机森林（R）randomFordom算法案例
机器学习算法之随机森林的R语言实现-表达芯片示例

本文思路

TCGA下载乳腺癌数据，把病人分成生存时间小于三年和生存时间大于五年两组，两组进行差异表达分析，然后用随机森林找出最重要的几个基因。

分析步骤

数据下载

library(TCGAbiolinks)
library(dplyr)
library(DT)
library(SummarizedExperiment)

#数据下载
#下面填入要下载的癌症种类
# if (!requireNamespace("BiocManager", quietly=TRUE))
#   install.packages("BiocManager")
# BiocManager::install("TCGAbiolinks")
# 
# 
# if (!requireNamespace("BiocManager", quietly=TRUE))
#   install.packages("BiocManager")
# BiocManager::install("TCGAbiolinksGUI", dependencies = TRUE)
# 
request_cancer=c("BRCA")
for (i in request_cancer) {
  cancer_type=paste("TCGA",i,sep="-")
  print(cancer_type)
  #下载临床数据
  clinical <- GDCquery_clinic(project = cancer_type, type = "clinical")
  write.csv(clinical,file = paste(cancer_type,"clinical.csv",sep = "-"))
  
  #下载rna-seq的counts数据
  query <- GDCquery(project = cancer_type, 
                    data.category = "Transcriptome Profiling", 
                    data.type = "Gene Expression Quantification", 
                    workflow.type = "HTSeq - Counts")
  
  GDCdownload(query, method = "api")
  expdat <- GDCprepare(query = query)
  count_matrix=assay(expdat)
  write.csv(count_matrix,file = paste(cancer_type,"Counts.csv",sep = "-"))
  
}

用DESeq2包进行差异分析分析

#载入包，并启用8核并行运算
library(DESeq2)
library(BiocParallel)
register(MulticoreParam(8))
#clinical=read.csv("TCGA-BRCA-clinical.csv",header = T,row.names = 1)
#count_matrix=read.table("TCGA-BRCA-Counts.csv",row.names = 1,sep=",")


#首先对病人进行分组
poor_OS_patients=clinical[clinical$days_to_death <= 1095 & clinical$vital_status == "dead" & is.na(clinical$days_to_death) == F ,]
poor_OS_patients[,1]=gsub("-",".",poor_OS_patients[,1])
good_OS_patients=clinical[clinical$days_to_death >= 1825 | clinical$vital_status == "alive",]
good_OS_patients[,1]=gsub("-",".",good_OS_patients[,1])

#处理表达矩阵
count_matrix_t=as.data.frame(t(count_matrix))
#去掉normal的数据
count_matrix_t1=cbind(row.names(count_matrix_t),count_matrix_t)
count_matrix_t1=count_matrix_t1[substr(count_matrix_t1[,1],14,15) < 10,]
#分割成两个矩阵
count_matrix_t2=cbind(count_matrix_t1[,1],count_matrix_t1)
colnames(count_matrix_t2)[1]=colnames(poor_OS_patients)[1]
count_matrix_t2[,1]=substr(count_matrix_t2$submitter_id,1,12)
poor_matrix=merge(count_matrix_t2,poor_OS_patients,by=colnames(poor_OS_patients)[1])
good_matrix=merge(count_matrix_t2,good_OS_patients,by=colnames(poor_OS_patients)[1])

#构建差异分析矩阵

database=cbind(t(poor_matrix),t(good_matrix))[2:56717,]
patients_barcode=database[1,]
colnames(database)=patients_barcode
database=database[-1,]
mode(database)="numeric"
condition <- factor(c(rep("poor_OS",length(poor_matrix[,1])),rep("good_OS",length(good_matrix[,1]))))
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
nrow(dds)
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds)
dds <- DESeq(dds,parallel = T)
res <- results(dds)
summary(res) 
#提取FDR<0.01且|logFC|>1的差异表达基因
table(res$padj<0.01 & abs(res$log2FoldChange)>1)
res <- res[order(res$padj),]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
signresdata<-na.omit(resdata[resdata$padj<0.01 & abs(resdata$log2FoldChange) > 1,])

对以上得到的差异基因进行随机森林分类

#随机森林
# install.packages("randomForest")
# install.packages("ROCR")
# install.packages("Hmisc")
# source("http://bioconductor.org/biocLite.R")
# biocLite("genefilter")

library(randomForest)
library(ROCR)
library(genefilter)
library(Hmisc)
#对差异基因表格进行处理
signres_norm_matrix=signresdata[,-2:-7]
row.names(signres_norm_matrix)=signres_norm_matrix[,1]
signres_norm_matrix=signres_norm_matrix[,-1]
signres_norm_matrix_forest=as.data.frame(t(signres_norm_matrix))

#制作性状表格
trait=as.data.frame(row.names(signres_norm_matrix_forest))
trait=cbind(trait,substr(trait$`row.names(signres_norm_matrix_forest)`,1,12))
colnames(trait)[2]="barcode"
good_trait=cbind(good_OS_patients[,1],rep("good",length(good_OS_patients[,1])))
poor_trait=cbind(poor_OS_patients[,1],rep("poor",length(poor_OS_patients[,1])))
merge_trait=as.data.frame(rbind(good_trait,poor_trait))
colnames(merge_trait)[1]="barcode"
trait=merge(trait,merge_trait,by="barcode")
trait=trait[,-1]
colnames(trait)[1]="barcode"
# barcode=as.character(trait$barcode)
# status=trait$V2
# trait=data.frame(barcode,status)
forest_trait=cbind(row.names(signres_norm_matrix_forest),signres_norm_matrix_forest)
colnames(forest_trait)[1]="barcode"
forest_trait=merge(trait,forest_trait,by="barcode")
colnames(forest_trait)[2]="trait"
row.names(forest_trait)=forest_trait[,1]
forest_trait=forest_trait[,-1]
#将数据集分为2/3训练集和1/3测试集,并查看数据集基本属性。
ind=sample(2,nrow(forest_trait),replace=TRUE, prob=c(0.67,0.33))
train=forest_trait[ind==1,]
test=forest_trait[ind==2,]
#选取mtry节点值,选择误差最小时对应的i值
#每次都set.seed一下，保证数据的可重复性
errset=1:(length(names(train))-1)
for (i in 1:length(names(train))-1){
  set.seed(350)
  mtry_fit=randomForest(trait ~ .,train,mtry=i)
  err=mean(mtry_fit$err.rate)
  print(err)
  errset[i]=err
}
best_mtry=which(errset==min(errset),arr.ind=TRUE)
#选择ntree，选择图像趋于稳定时的树值
set.seed(350)
model <- randomForest(trait~., train, mtry=best_mtry[1], ntree = 1000)
plot(model)
#大约400棵树比较稳定了

#训练集进行随机森林
set.seed(350)
rf <- randomForest(trait~., train, mtry=best_mtry[1], ntree = 400, importance = TRUE)

#查看变量的重要性
importance <- importance(x=rf)

#绘制变量的重要性图
varImpPlot(rf)

#最后验证并预测

pred1<-predict(rf,data=train)
Freq1<-table(pred1,train$trait)
sum(diag(Freq1))/sum(Freq1)
plot(margin(rf,test$trait))
#全体数据用于训练集
set.seed(350)
rf_output=randomForest(trait ~ ., forest_trait, importance = TRUE, mtry=best_mtry[1], ntree = 400, proximity=TRUE)

#查看变量的重要性
importance <- importance(x=rf_output)

#绘制变量的重要性图
varImpPlot(rf_output)

image.png

最后得到的结果，左边是基于精确度的重要性排序，右边是基于基尼系数的重要性排序，最右上的重要性最高。附：重要性的表格

按精确度排序

image.png

按基尼系数排序

image.png

后续待做的工作

以全部数据作为训练集弄出来重要的基因以后，还需要寻找合适的数据，作为测试集，去验证我们的分析是否具有普适性。数据我现在懒得找。

总结

本文的矩阵处理过程略显捉急，自己都觉得有点蠢，不过能够达到目的就行了，反正我也不要求代码简洁，也不要求运行效率多高，出来正确的结果就行。
另外，我只是依样画葫芦，对真正的随进森林算法的理解非常浅薄，还要日后恶补统计学知识，准备抽空看一下《统计学完全教程》，据说很不错。