我的生信之路-2019（完整版）
学习这个过程，一个月可能就是一个样，所以生信之路每篇展示一年的生信学习记录。这是19年的记录，由于20年初和19年底之间是一个整体，所以一起算在内了，这次的记录长达15个月。不同于第一篇。
写在前面 （2020年5月）
没想到新冠病毒的影响力一直持续到了现在，无论是我们的何种工作，或多或少的受到了影响。生物信息是一个大的题目，我起名太大了，我所谓的生信，仅仅也就是多组学数据的高效运作和常见统计分析的高效处理。要学习的东西实在是太多了，一句话，两句话，说不清楚，一年的经历，足够写一长篇笔记，借用腾讯的平台，留给初学者一些文字。也是很好。
为什么我要继续补充我的生信之路呢2呢? 因为上次写的时间在19年11月，其实有两个月时间没有计算在内，再加上新冠病毒，让20年的一个季度和19年的学习变成了一个整体，我有必要进行一个汇总，也算是真正结束这一年多的学习。
写下前面 （2019年底）
现在已经是19年了，年底也快到了，各地的会议也在朋友圈此起彼伏。想想我做生物信息也三年了。故事总在发展着，去年我总结了我的生信之路，写了我是如何开始做生信的。今天我将为大家带来新的一年来的生信之路。近几个月来，我更新的速度比较慢，一方面是科研任务比较繁重，其实有许多东西学习需要很长时间才能和大家见面，所以我将邀请我的朋友五谷杂粮协助运营。努力为大家带来更好的知识分享和技术进步。
引子：
有一句话我想先讲：在这个年纪，最忙的永远是明天。微生信生物公众号成立也快两年了，到现在也在不断的见证着我学习的过程。但是被大家熟知的是我总是不定期的来给大家展示一些有意思的东西。不能很好的对大家负责，自由散漫，所以大家留下来的都是我的亲人呐。
开始先攻坚代谢组，无奈的是相关资料较少，学期起来十分吃力，尤其是在拿到物质丰度矩阵后，因为需要坐很多的后续分析，牵扯到pca排序，批量差异分析，机器学习等等，吃力的学习了一段时间后发现其实组学后续分析的方法都是类似的，所以我后来就准备从扩增子数据开始做。花费我很大功夫的确是也是扩增子测序分析。为什么呢？因为扩增子比较成熟，学习资源和文章都比较多，尤其是后续数据分析，我可以跟着别人的中文教程来模仿。事实也证明我的选择没有错，真正让我入门生信的也是扩增子测序的流程。其实转录组流程也挺成熟的，最起码在我当时入门生物信息学的16年已经如此了。所以，结合自己的课题，选择合适的组学做为自己生物信息之路的第一扇门。
宏基因组公众号给我我很大的帮助，我的入门也算是伴随着宏基因组公众号的开端。刘老师的无私奉献让我确实在最短的时间就入门了扩增子。到现在我也成为了宏基因组的一员，饮水思源，尤其是生物信息这个行业，所以，现在我跟着刘老师一起来帮助大家还在摸索的新生力量。开端在之前一篇公众号推文：我的生信之路已经将清楚了我在刚开始的模糊摸索之路。我在这里提一提的目的还是想告诉新人，沉下心来，总会入门的，我用了两年时间，可能你们只需要几个月就会喜欢上这门学科，走向正轨。
正文：
18年的冬天，南京还是没怎么下雪，天气也不算很冷，可能在这里经历过两个冬天后慢慢习惯上这里的冬天没有暖气的生活。伴随着朝阳，我也开始了一天的忙碌。
如果说之前的学习是为了学会如何做分析的话，现在我更多的是面对一堆数据如何选择合适的分析方法和流程整理出最合适的结果。这句话也从18年底一直指导我工作到了今天。
扩增子分析进一步学习
扩增子测序描述微生物群落已经快十年了，其实真正开始大规模在土壤学运用也在12年以后了，记得老院长曾在13年花了10万测了10个扩增子样品。现在想想就是10个宏基因组样本也用不了1万了。确实测序发展很快，一直很快。伴随着测序成本的降低和通量的提高，更多的有关扩增子测序的文章在16年以后铺天盖地的涌现。如果说一门新的技术的运用在刚开始是出于狂热期，也就是更多的去关注技术本身，而不是关注真正的科学问题，那么现在人们扩增子测序已经变得理性了，更多的是如何使用这门技术来解释自己的科学问题。慢慢的这门技术也走向了常规化。从这个状况来看，其实我开始学习扩增子测序正好处于研究的狂热后期，拽着这个最后的劲，我赶了一趟末班车。
我并不是专职的分析人员，我最需要的不仅仅是掌握分析技术，而是解释科学问题。其实再16年，我是不理解的，当时认为我会分析我就牛逼啦，会分析了，我还需要什么？这样麻痹了一年多，沉浸在自己勾画的小圈子里慢慢享受学会分析的乐趣。终于17年底这个美梦还是醒了。我需要写文章了，这是大多数人科研的目标，我们又何尝不是想要一篇篇高分文章呢！当时数据回来后我是很高兴的，蹭蹭蹭的我就将多种分析做完了，这些结果我拿着微信截图都摆到了ppt上，兴冲冲的就去汇报了。拿着一堆分析结果去汇报，我竟然啥都讲不出来，这是比较尴尬的，为什么我讲不出来呢？尴尬的来回翻了翻ppt，最后悻悻的下去了。师兄问：有什么结果吗？我记得当时懵了一下，说：群落有差异。哈哈哈，想起来往事确实挺有意思。但是我也陷入了深深的思考。我在之前更多的关注分析的过程了，更多的关注如何学会分析，如何run出来而不是出现error。而现在更多的关注结果如何，如何根据结果调整分析方法，或者更换流程。而这个理解也是我从学分析到学会分析的转变的开始。
当时我总结的扩增子流程：一方面参照刘老师的usearch流程，后续使用R语言做的基础多样性和差异分析，一方面参考超哥的指导弄了虚拟机，做一些高级分析。所以我当时整合了Qiime1的原始数据处理和后续刘老师公众号的多样性和差异分析流程。这两方面让我学会了使用生物信息两大重要的工具：R语言和linux系统。之后的工作我便大部分都布置在这两个环境了。Qiime1的使命似乎快要终结，但是usearch在作者的不断更新下持续换发的生命，usearch也从当初一个嵌合体算法变成如今的一整套处理扩增子序列的pipeline。那么qiime1没用吗？不是这样的，qiime1虽然已经不在维护了，但是当时它的理念是尽可能囊括扩增子数据分析的全部所需工具。所以，很多基于fastq文件和fasta文件的操作都很全面，甚至还有一些非标准格式的测序数据处理和转化工具。只是基于python2的环境和一些旧包的依赖让这个大的工具体系逐渐开始被替代。
新一代扩增子分析T体系
18年是Qiime2公测的一年，各大平台讨论过很多Qiime2。我也凑过去看了一些。也就是那个时候我开始有更新扩增子流程的打算。最主要的就是将基于聚类获得otu的方法更换为非聚类的方法。由于在usearch中的unoise3算法比较快，所以我首先测试并添加到我的流程当中。但是DADA2算法之前吵得很热。unoise3学习完成之后我就开始了DADA2的学习了。当时我只知道这个算法可以在qiime2上实现，所以就使用qiime2做了个测试，然而在虚拟机上的速度确实很慢，我也无法看到中间过程，这让我在等待的过程中很焦急。所以我就去了解到这是个R包，就学习了DADA2包，并使用R包做了一套扩增子原始数据处理流程。这个流程好的一点就是可以按照自己的需求进行细部的修改，尤其是错误模型的运行参数，在qiime2中修改有限。由于在win下的R并不能使用多线程，所以在linux才能显示出这套流程的真正速度。
借助这个算法，我接触到了phyloseq包，这是一篇16年发表的文献，整合了扩增子分析需要的五大文件，也就是我们经常用到的otu表格，样本分组文件，代表序列，进化树，物种注释，并以S4类对象存储。也就是说数据被打包了，这是一种新的数据保存策略，我认为最方便的就是可以同步修改数据，也就是可以同步过滤otu，合并otu，样本，甚至是进化树。当然存储对象就不是很好理解了，需要一定的入门学习。
当时我认定这个包很重要。所以就开始慢慢摸索使用这个包来做一些分析，将自己原来流程的输入全部使用phyloseq格式打包。慢慢学习逐渐替换掉自己别的方式的分析脚本。然而时间过得很快，又因为有大量的科研任务，所以学习的速度也不够快 ，还断断续续的，中间一段翻译了部分教程，也发过几次到微生信生物上，做的不够完善，就这样一点一点的学习时间到19年中了。六月份的天是美丽的，不热不冷，梅雨刚过，天气也不好不坏。头脑中积累的东西一直在推着我有这样的想法出现：我要做一套扩增子分析流程。虽然扩增子流程已经很成熟了，但是我们仍然无法理解庞大的微生物群落之间的联系，这不是单单扩增子可以解决的问题，宏基因组，多组学联合慢慢变成了趋势；将多组学整合到统一的框架，并且思考如何运用多组学数据解释群落问题早就变成了现实。我想以扩增子为入手点，逐渐展开。
想法是美妙的，但是大量的科研任务在身，让这个计划推迟了许久，终于在七月份，我启动了这个计划。开始了为期两周的魔鬼式写代码。从原始数据处理开始我一共构建了十几份代码，以phyloseq为框架，整合了DADA2算法为核心的数据预处理；以phyloseq为核心的原始数据的修正和数据标准化；以microbiome和phylodeq包为框架整合超过30种alpha多样性分析流程；以phyloseq包中超过40种算法为基础的八种排序分析；使用deseq2，edgr，非参数检验等多种分析方法构建的差异分析一体化流程；群落和环境因子关系探索的全套流程，包括显著性排序-基于多种目的筛选物种同环境因子相关等。同时，我有意构建高水平的可视化部分：基于ggplot展示物种分类关系的冲击图，和弦图等，基于R语言的高效graphlan注释系统结合graphlan工具展示完美进化树和物种进化树。基于R语言和gephi构建漂亮的二分网络。基于R语言的全套网络分析和可视化方案。就这样一期工作做完了。一期工作成果已经开源到github上了，一个叫easy_microbiome项目。希望大家多多参与，并提出自己的意见。后我又开始了二期工作：基于R语言多种机器学习在群落中的运用和评估；细菌真菌联合分析方案探索和网络构建：不同研究群落数据整合方案及其评估。而这些工作一直持续到了今天。前两期动工作都做的是基本流程，努力将目前已有的分析框架整合到同一个linux平台上，并形成可供分享的iso或者vdi镜像文件。
七月底在刘老师培训班的继续学习让我又有了新的思考，但是由于工作的需要，这些想法目前还在成熟中，也就放下了。慢慢的时间就到最近了。九月是金秋，南农的菊花展每年都很惊艳。而我也看了两年了。如此多的菊花确实是我在别的地方看不到的，每次去图书馆我都会看看再走。让心灵得到纯净。
一周完成转录组全套
最近由于课题需要，我需要做一些转录组的工作。所以我又断断续续的阅读一些转录组的知识，做为储备。九月中旬，本科舍友转录组数据下来了，出于帮忙，也出于练手，我答应帮忙构建一套转录组流程。但是马上假期就来了，同时本科舍友婚礼，让我好一把玩到了十月上旬结束。期间见到老院长，也在头疼一批转录组数据，我们聊了聊，也想帮老院长解决一些问题。可是研究任务也愈发加重，回来后我还要继续弄自己文章的分析和写作，一直忙活到现在。这些琐事压在心里似乎无法排解也无法解决，于是，这样的计划便形成啦，上周末两天，加上这几天每天晚上抽出两个小时，我急行军似的开始构建一套转录组分析流程。因为想到转录组已经很多人发教程了，我虽然没有做过，但是幸好之前了解过一些，就这样愣头愣脑的开始了。
但非模式物种的有参转录组还是参考比较少的。我是用hisat2比对基因组却没有索引文件，这个过程生信技能树和生信媛的推文给我了一些帮助，我得已完成索引够构建和比对工作。使用一些扩增子用不到的工具做了一些数据转化后就开始进行基于R语言的后续数据分析了。
说起来现在转录组一个双端测序一个文件不到10G，让我这个经常处理扩增子的感到了压力，我用的这台电脑32G内存8个线程，勉强完成唉，也不知道同行什么配置，sam文件每个基本都超过了15G，最终全部比对完花费了一晚上和一早上。占用了我500G硬盘空间，也怪我没有将sam文件及时转化为bam文件，毕竟可以减少至少三分之二的空间。终于在周末完成这件事情。后续基于count矩阵的排序分析，差异分析我就轻车熟路了。基于pca的排序，deseq2差异分析不到半个小时就完成了，因为有现成的代码。
但是最麻烦的还是来了-富集分析。我最终做了四个富集分析：GO，KEGG，GSEA，GWVA。其中go和kegg是常规富集，有很多参考资料，但是我的物种不是模式作物，所以旧麻烦一些，幸好网站上有基因对应的go注释，否则这一趟注释下来又得好一会，通过Y叔的教程手动构造背景数据库得以最终完成go富集。Y叔的几个包成了我的宝贝，最后依靠enricher函数完成GO富集分析 。
KEGG数据库中是可以找到我的物种的，但是我选择的基因组却和数据库中注释的基因id不同，在时间有限的情况下我是用了ko模式做KEGG富集分析。也就是说species选择了ko，将我自己的基因注释到了对应的ko号上进行操作。最终完成kegg分析，并将差异倍数成功映射到通路上，输出pdf和xml文件。在这里我有个问题：R语言什么工具可以输出kegg每个通路可交互的html文件，我一直没找到类似的工具，时间也有限，所以这个问题就没解决？希望大家伸出援手呀。
GWEA分析考虑到的GO分析的短板，让基因差异按照一定顺序排布，在Y叔推送帮助下通过仔细设置排序好genelist得以解决输入问题，剩下的就和go富集一样了；gwea富集我用了差异倍数，这也是不同于富集的地方，让我们可以识别上调和下调的通路，比传统使用全部上调或者下调基因基因做富集更加优越。
GWVA考虑到了丰度，分析内容和可视化也更加丰富，和picrust功能预测一样，都是得到每个通路通过计算得到的值来判断不同样本对应通路的变化，所以这里是可以得到通路矩阵的不同与以上三中富集方法，这里的矩阵可谓是打开了又一套分析流程，基于通路矩阵的排序，差异分析，注意这里不能使用deseq2或者edgr，使用linma，因为这不是表达矩阵啦。真正做起来才发现也是要构建背景基因集合的，这个背景基因不同于以上三种富集的背景，这是一个kist文件，每个对象以通路命名包含通路内的基因。找了google来搜索下载kegg基因及其通路的注释，最终构造成功。还好最后都完成了。
迅速给这些结果加上可视化和文本输出，做了个循环，便构建好了这整套的非模式作物，全部要自己构建背景注释基因来坐富集的RNA-seq流程。已经发布在微生信生物，但是还是显得粗糙。
游走于宏基因组
长期以来，我对宏基因组都是三分热度，过一段时间做一点，学一点，所以一直以来都没有什么大的进步，直到19年11月，science发表的甜菜枯萎病宏基因组数据，让我真正去啃这块硬骨头。为什么说宏基因组难啃呢？
硬件资源跟不上，需要最起码也对几十个线程，100G以上内存才可以小做一番。
内容十分多样，虽然都是测序数据，但是处理数据质控，拼接一样之外，和扩增子的后续流程差异太大。
软件没有扩增子那么成熟，流程衔接也不是很容易。
数据库众多，扩增子也就几个物种注释数据库，没什么要理解的地方，但是宏基因组有NB，NR，nog，KEGG，GO，等等与功能相关的大库，小库就有十几个，物种库也按照不同软件有许多。
这套数据库配制下来没有一个T是弄不好的。
宏基因组包括扩增子流程，包括转录组流程，单菌基因组流程，如果预测了代谢组后还要用一部分代谢组流程。
在19年11月，我终于将以上四个部分准备完毕：
老板升级了服务器
基础知识学习到位，各种软件了解到位，作为宏基因组成员，有刘老师大量的学习资料。
张磊老师帮助配制了全部数据库。
同时对张老师流程进行了试运行。
在11月之前，我已经完成了扩增子，转录组，代谢组流程的工作。
所以我开始了半年来的宏基因组的学习，其实这段时间额学习我是以11月发表的science作为模板学习的，这份教程长达1800行，已经作为独立的流程发表了。按照我近半年来的梳理，这份教程确实十分完善，比我见过的任何一份教程都完善，其次，这份教程还有众多的分析选项和思路在里面，值得我们不断的回味。大家可以试试学一下，目前我也就学习到了1000行，也是主要的流程。还有800行等着我，主要是后续功能和物种的后续处理，这部分我可能会相对熟悉。这里也向大家问好，这1800行最终在今年年底和大家见面，见一下目前我认为最完善的宏基因组分析的老一套典型流程。
这份教程用snakefile做的流程，大部分数据转换用的是python脚本，和小部分shell。不知道为什么，微生信生物总是这样，只要认定的事情，就一定要做到。我半年来开始了shell的学习，pyhton的学习，只为学习这一份代码。虽然这份代码已经好多步骤不能跑起来了，但是我都使用了新版本软件，或者替换程序作了修改，只为重现流程。这一件事情就花费了我半年的时间，相信出来之后大家就会看到微生信生物在对待宏基因组到底有多么认真。
就这样这些工作一直持续到了4月，期间我在宏基因组发布的MEGAN学习笔记，也是来自其中对的一部分，也是目前改版最大的一个部分。开学日期将近，我准备东西走向了学校，宏基因组的学习也告一段落了。这才真正完成了19年的生信学习记录。
写在后面 19年11月：
今年生物信息学习明显比之前轻松了，这是客观的，但是要学的东西却源源不断，扩增子让我了解了群落，了解了组学得后续分析内容和部分通用方法，那么宏基因组将带我们走向功能的研究。本来我想的是先构造一个适合自己的宏基因组流程。但迫于需求和研究需要，还是先做了转录组。今年年底之前我准备讲宏基因组也整上一遍。毕竟组学之间不是独立的，如何多组学整合是我们要不断理解和解决的问题。也是我的兴趣所在。
写在后面 20年5月：
越来越多的工作，让我意识到再也没有足够的整块时间，像这次年底学习一样了，我面对的是无数的实验，头疼的分析，还有最不擅长的文章的写作。总要走这一步的，精力划过去之后，在生信的时间就慢慢少了，随着我的成长，对于初学者的友好性也将逐渐消失，我最近已经在思考如何处理这些问题，和让微生信生物持续稳定发展。
最后还是那句话：学习永不止步，分享永不停歇！