要保证当前文件夹下面有了742KO1.count等4个文件，就是用htseq对比对的bam文件进行处理后的输出文件。
library(DESeq) #加载数据 K1=read.table(“742KO1.count”,row.names=1) K2=read.table(“743KO2.count”,row.names=1) W1=read.table(“740WT1.count”,row.names=1) W2=read.table(“741WT2.count”,row.names=1) #列名 data=cbind(K1,K2,W1,W2) #如果是htseq的结果，则删除data最后四行 n=nrow(data) data=data #如果是bedtools的结果，取出统计个数列和行名 kk1=cbind(K1$V5) rownames(kk1)=rownames(K1) K1=kk1 #差异分析 colnames(data)=c(“K1″,”K2″,”W1″,”W2″) type=rep(c(“K”,”W”),c(2,2)) de=newCountDataSet(data,type) de=estimateSizeFactors(de) de=estimateDispersions(de) res=nbinomTest(de,”K”,”W”) #res就是我们的表达量检验结果 到这里，理论上差异基因的分析已经结束啦！后面只是关于R的bioconductor包的一些简单结合使用而已 library(org.Mm.eg.db) tmp=select(org.Mm.eg.db, keys=res$id, columns=c(“ENTREZID”,”SYMBOL”), keytype=”ENSEMBL”) #合并res和tmp res=merge(tmp,res,by.x=”ENSEMBL”,by.y=”id”,all=TRUE) #go tmp=select(org.Mm.eg.db, keys=res$ENSEMBL, columns=”GO”, keytype=”ENSEMBL”) ensembl_go=unlist(tapply(tmp[,2],as.factor(tmp[,1]),function(x) paste(x,collapse =”|”),simplify =F)) #为res加入go注释， res$go=ensembl_go[res$ENSEMBL]#为res加入一列go #写入all——data all_res=res write.csv(res,file=”all_data.csv”,row.names =F) uniq=na.omit(res)#删除无效基因 sort_uniq=uniq[order(uniq$padj),]#按照矫正p值排序 #写入排序后的all_data write.csv(res,file=”all_data.csv”,row.names =F) #标记上下调基因 sort_uniq$up_down=ifelse(sort_uniq$baseMeanA>sort_uniq$baseMeanB,”up”,”down”) #交换上下调基因列位置 final_res=sort_uniq[,c(12,1:11)] #写出最后数据 write.csv(final_res,file=”final_annotation_gene_bedtools_sort_uniq.csv”,row.names =F) #然后挑选出padj值小于0.05的数据来做富集 tmp=select(org.Mm.eg.db, keys=sort_uniq[sort_uniq$padj<0.05,1], columns=”ENTREZID”, keytype=”ENSEMBL”) diff_ENTREZID=tmp$ENTREZID require(DOSE) require(clusterProfiler) diff_ENTREZID=na.omit(diff_ENTREZID) ego <- enrichGO(gene=diff_ENTREZID,organism=”mouse”,ont=”CC”,pvalueCutoff=0.05,readable=TRUE) ekk <- enrichKEGG(gene=diff_ENTREZID,organism=”mouse”,pvalueCutoff=0.05,readable=TRUE) write.csv(summary(ekk),”KEGG-enrich.csv”,row.names =F) write.csv(summary(ego),”GO-enrich.csv”,row.names =F)
原文来自：http://www.bio-info-trainee.com/867.html