看白癜风哪家医院好 http://wapyyk.39.net/bj/zhuanke/89ac7.html
一份高度定制化个性化的富集分析教程，将解决你做富集分析的全部问题。

写在前面

上周末两天，加上这几天每天晚上抽出两个小时，我急行军似的开始构建一套转录组分析流程。因为想到转录组已经很多人发教程了，我虽然没有做过，但是幸好之前了解过一些，就这样开始了。但非模式物种的有参转录组还是参考比较少的。我是用hisat2比对基因组却没有索引文件，这个过程生信技能树和生信媛的推文给我了一些帮助，我得已完成索引够构建和比对工作。使用一些扩增子用不到的工具做了一些数据转化后就开始进行基于R语言的后续数据分析了。说起来现在转录组一个双端测序一个文件不到10G，让我这个经常处理扩增子的感到了压力，我用的这台电脑32G内存8个线程，勉强完成唉，也不知道同行什么配置，sam文件每个基本都超过了15G，最终全部比对完花费了一晚上和一早上。占用了我G硬盘空间，也怪我没有将sam文件及时转化为bam文件，毕竟可以减少至少三分之二的空间。终于在周末完成这件事情。后续基于count矩阵的排序分析，差异分析我就轻车熟路了。基于pca的排序，deseq2差异分析不到半个小时就完成了，因为有现成的代码。但是最麻烦的还是来了-富集分析。我最终做了四个富集分析：GO，KEGG，GSEA，GSVA。其中go和kegg是常规富集，有很多参考资料，但是我的物种不是模式作物，所以旧麻烦一些，幸好网站上有基因对应的go注释，否则这一趟注释下来又得好一会，通过Y叔的教程手动构造背景数据库得以最终完成go富集。Y叔的几个包成了我的宝贝，最后依靠enricher函数完成GO富集分析。KEGG数据库中是可以找到我的物种的，但是我选择的基因组却和数据库中注释的基因id不同，在时间有限的情况下我是用了ko模式做KEGG富集分析。也就是说species选择了ko，将我自己的基因注释到了对应的ko号上进行操作。最终完成kegg分析，并将差异倍数成功映射到通路上，输出pdf和xml文件。在这里我有个问题：R语言什么工具可以输出kegg每个通路可交互的html文件，我一直没找到类似的工具，时间也有限，所以这个问题就没解决？希望大家伸出援手呀。GSEA分析就相对高档一点了，在Y叔推送帮助下通过仔细设置排序好genelist得以解决输入问题，剩下的就和go富集一样了；gsea富集我用了差异倍数，这也是不同于富集的地方，让我们可以识别上调和下调的通路，比传统使用全部上调或者下调基因基因做富集更加优越。GSVA据说很好唉，我想试试，为什么呢？看到这个富集的简介，我感觉和扩增子的picrust功能预测差不多嘛！都是得到每个通路通过计算得到的值来判断不同样本对应通路的变化，所以这里是可以得到通路矩阵的不同与以上三中富集方法，这里的矩阵可谓是打开了又一套分析流程，基于通路矩阵的排序，差异分析，注意这里不能使用deseq2或者edgr，使用linma，因为这不是表达矩阵啦。真正做起来才发现也是要构建背景基因集合的，这个背景基因不同于以上三种富集的背景，这是一个kist文件，每个对象以通路命名包含通路内的基因。找了google来搜索下载kegg基因及其通路的注释，最终构造成功。还好最后都完成了。迅速给这些结果加上可视化和文本输出，做了个循环，便构建好了这整套的非模式作物，全部要自己构建背景注释基因来坐富集的RNA-seq流程。在此晒出来，有经验的朋友给我指正指正。

#设置工作目录，所有文件放到此目录下path=getwd()path=paste(path,"/result",sep="")dir.create(path)mytheme-theme_bw()+#theme_classic()+#scale_color_manual(values=mi,guide=guide_legend(title=NULL))+#scale_fill_manual(values=mi,guide=guide_legend(title=NULL))+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),plot.title=element_text(vjust=-8.5,hjust=0.1),axis.title.y=element_text(size=20,face="bold",colour="black"),axis.title.x=element_text(size=24,face="bold",colour="black"),axis.text=element_text(size=20,face="bold"),axis.text.x=element_text(colour="black",size=24),axis.text.y=element_text(colour="black",size=24),legend.text=element_text(size=20,face="bold"))+#theme(legend.position=c(0.1,0.2))+theme(strip.text.x=element_text(size=15,angle=0),strip.text.y=element_text(size=12,face="bold"),strip.background=element_rect(colour="blue",fill="#CCCCFF"))

library("DESeq2")#-----------------------------------------导入数据---------------------------------------------------countdata-read.delim("CountMatrix.txt",row.names=1)head(countdata,10)#构造分组文件SampleType-sapply(strsplit(colnames(countdata),"_"),`[`,1)SampleTypecoldata=data.frame(ID=colnames(countdata),SampleType)head(coldata)Desep_group-as.character(levels(coldataSampleType))Desep_groupaaa=