下载完RNA-Seq数据后，我们还需要下载一个拟南芥cDNA序列(记住是组，而不是全基因组，好了，你别说了，我记住了)
当然你完全不必真的去睡午觉，我们可以程序运行的时候准备一下tximport所需要的输入文件。
tximport可以纠正不同样本基因长度的潜在改变（比如说differential isoform usage）；能够用于导入 (Salmon,Sailfish,kallisto)程序的输出文件；能够避免丢弃那些比对到多个基因的同源序列，从而提高灵敏度
虽然tximport的参数看起来很多，但其实需要我们准备的就是两个files和tx2gene
吐槽: 原本的我以为sample从1到16会和194到206是一一对应，但是万万没想到，他的对应关系居然是下图这个样子的。看来凡是都不能想当然，一定要仔细看他们的对应关系。
由于后续要用DESeq2包做差异表达分析，所以需要用DESeqDataSetFromTximport这个函数，当然你还需要说明你的实验设计
最后找到了208个差异表达基因，105上调，103个下调；而文章是298个，148上调，156下调。
由于这个画图过程，后续会用大很多次，为了避免重复操作，我写成了一个函数，我你根据的流程自己写一个函数，不要用我的。
和文章的结论基本一致，早花同源异形基因在我的流程中只能检测到CAL, 而发育中期和后期的基因基本上表达量也特别低。不过文章说CAL3这个基因能被稳定被检测到，实际上表达量与这些基因没有明显区别（< 100)，我也不知道他如何得到这个结论的。
很尴尬，显著性上调的基因也都在我的列表中，但是显著性下调的，一个都没有，于是我把水平从0.01提高到0.05，重新检查下
好吧还是只找到一个，毕竟ELF4这个表达量实在是太低了，除非水平变成是0.1，不然是难以发现的啦。但是这个还是证明作者的猜想是对的。
为了验证这个结果是不是正确的，我使用Y叔之作clusterProfiler，进行GO富集分析。
结果在上调中并没有找到富集基因，但是在下调（day1 vs day0 和 day3 vs day0)中却发现了好多富集基因，虽然场面十分尴尬，但是至少说明下调的基因大多属于ribosome
或许我可能需要用文章说的AmiGO 2 versions 2.3.2 and 2.4.24进行GO富集分析，或许用DAVID，其他工具分析一遍，感觉又要许多主观能动性了。
在这次实战中，我学习了用tidyr, dplyr对数据进行预处理，重新看了一下ggplot2的图形语法，复习了之前差异表达分析的基本操作。
最新评论