16SrRNA 基因通常作为分子标记进行微生物群落结构的研究，但是它有一些明显的限制，比如16S rRNA基因在物种中会有多个拷贝，而且，由于16S rRNA基因的进化速率较慢，在物种间保守，会存在多个物种的基因完全相同的情况，而且由于基因水平转移的发生，即使亲缘关系较远的物种，也可能出现基因序列完全相同的情况，更进一步讲，我们分析时通常只采用16S rRNA基因的某些区域，这导致物种间扩增出来的片段完全相同的概率大大增加；
而一些蛋白编码基因，特别是一些参与重要的信号通路的基因，比如参与氮循环的的基因，这些基因出现水平转移的概率小，也可以作为分子标记来研究微生物的群落结构 。FunGene 是一个免费的数据库，收录了许多功能基因的序列，而且提供了一些工具对功能基因进行分析。
FunGene的序列来源于GeneBank 数据库，而GeneBank 数据库是有冗余的，所以FunGene 也会有冗余现象，所以在下载完序列之后，需要去冗余。在去冗余的过程中需要注意，uncultured 代表无法单独培养的一些菌， 比如 uncultured bacterium 代表无法单独培养的细菌，这个是一类物种的简称，不是指同一个物种。
一个官方的fungene分析流程来自他的官网 RDP Tutorials分析流程有在线版FunGene Pipeline和本地版（托管在GitHub上）rdpstaff。本地版结合RDPTools可以做更多事情。
fungene
fungenepipeline主要有以下五个步骤：质控、蛋白质校正、比对、聚类和多样性分析。
具体地，在其官网上分为了12步
关键的一步就是framebot，即氨基酸的校正与过滤。为什么功能基因要用氨基酸序列分析呢？
- 功能基因多样性研究的目的是功能生态学方向的研究
- 功能基因扩增子序列相较于rDNA更复杂（GC%高）
- 部分功能基因在DNA水平具有较低的基因分辨率和相似性
- 取出非功能基因序列、含终止子的序列以及氨基酸的校正（移码突变、读框移位等）
功能基因多样性分析一般计算出多样性指数jaccard_and_sorensen、shannon_chao1等就可以做样本间和组间的差异分析了。这个pipeline并没有对氨基酸做物种的注释。
taxadiva
下面我们来看两个以nifh为例的已发表的pipeline。一个是 2017发表在Applied and Environmental Microbiology上的，并没有基于蛋白序列来分析而是直接拿扩增的nifh序列来物种注释。文章提供的pipeline托管在github上taxadiva，
我们看他的分析流程上并无新意，但是他的注释过程却值得我们注意：
三种注释情况的ROC曲线来判断注释的结果。
注释的数据库nifH Sequence Database是比较早的数据库（2014，contains 32,954 aligned nifH sequences with associated metadata）
为了使用UNIFRAC进行β多样性分析，pipeline重建了来自每个OTU的代表性序列的系统发育。 首先，使用QIIME命令pick_rep_set.py挑选一组代表性序列，并将这些序列与parallel_align_seqs_pynast.py对齐nifH参考比对。 使用QIIME命令filter_alignment.py消除所有间隙的对齐列，并使用QIIME命令make_phylogeny.py制作系统发育树。 在QIIME中计算α和β多样性指标。 为了创建丰度条形图，TaxADivA脚本用于为OTU分配分类并输出OTU表，然后将其转换为hdf5格式并用作QIIME脚本summarize_taxa_through_plots.py的输入。
当然有了以上结果，后面的统计分析也都可以用R，python、spss来是实现了。
NifMAP
Evaluation of Primers Targeting the Diazotroph Functional Gene and
Development of NifMAP – A Bioinformatics Pipeline for Analyzing
nifH Amplicon Data
NifMAP（“NifH MiSeqIllumina Amplicon Analysis Pipeline“），用Hidden-Markov模型将同源基因过滤掉（氨基酸序列的校正）。
（1）使用QIIME的join_paired_ends.py将成对的原始MiSeq读数组装成重叠群【contigs】（Caporaso等，2012）。
（2）利用HMMER中的HMMQuestH命令，将合并后的重叠群与基于HMM的HMM（HMMN-NUCUC1160SnIFH.HMM）进行比对过滤。不管模型匹配分数如何，通过这一步骤的所有reads都被接受用于下一步骤。
（3）使用UPARSE（Edgar，2013）对序列进行嵌合过滤和聚类。
首先，使用-derep_fullength命令对contigs进行去重，并删除单独的唯一序列。 然后使用-cluster_otus命令（using 3% radius）确定OTU质心。
通过使用-usearch_global命令（以0.97％的同一性，以下称为OTU97）将过滤的重叠群（在去重复之前）映射到OTU质心来确定OTU的丰度。
（4）使用FrameBot（Wang等人，2013）针对nifH蛋白参照集进行氨基酸翻译和（潜在的）移码校正。
（5）使用HMMER中的hmmscan命令，针对nifH基因（bchX，chlL，bchL和parA）的同源基因针对HMM nifH_ChlL_bchX.hmm（见上文）进行过滤。 只保留了与bchX和chlL-bchL模型相比，对nifH模型得分最高的OTU代表序列。
（6）OTU分类和系统发育构建：剩余的OTU代表序列使用BLASTP（Camacho等，2009）RefSeq数据库进行分类学注释（Pruitt等，2005）。
其实到这里功能多样性差异检验等的分析就可以做了，这流程的一个比较新的点是他用校正的nifh基因构建进化树。
OTU代表也使用RAxML中的进化树算法（EPA）实现放置在基础树上（Stamatakis，2014）。 基础树生成如下：
（1）提取含有nifH_2014April04.arb数据库中nifH基因的氨基酸序列信息的所有条目。 其中包括来自基因组的1971个条目和来自非基因组测序来源的39,258个条目。
（2）过滤掉所有短于133AA的序列，并将剩余的序列去除重复。
（3）使用CD-HIT将剩余序列聚集在90％同一性的聚类阈值（Fu等人，2012）。 此外，使用CD-HIT将用于构建bchX和chlL-bchL HMM的序列聚集在80％同一性的阈值，并与聚类的nifH序列合并。 然后使用MAFFT L-INS-i将组合的数据集与用于构建HMM Zehr_2014_1812genomes_nifH_AA.hmm的基因组起源的1812个氨基酸序列的比对集合进行比对，以保持与nifH_2014April04.arb数据库相容的比对数。 最后，使用RAxML重建基于CAT模型的自举最大似然树（ bootstrapped maximum likelihood tree ）。
为了将OTU代表放置在基础树上，使用MAFFT对齐用于构建基础树的对齐，然后使用RAxML使用EPA将序列添加到基础树。 pipeline的步骤1和3与Herbold等人描述的程序相同。 （2015年），而第2,4,5和6步是这项工作的新内容。 用于再现步骤2,4,5和6的HMM，基础树和shell脚本可在https://github.com/roey-angel/NifMAP上公开获得。
NifH的无根RAxML-EPA树推导的氨基酸序列（具有至少133个位置）。 nifH和非nifH（同源基因）簇以不同的颜色表示。 来自“nifH_2014April04.arb”数据库（Heller等，2014）的代表性序列用黑色表示。 本研究中来自三个引物对的OTU代表的系统发育位置描述于末端节点：用引物对Ueda19F-R6获得的序列以红色显示，IGK3-DVV衍生的序列以蓝色显示，以及F2-R6衍生的序列 黄色 最外面的条代表每个OTU的相对丰度[log（x％+ 1）]。
文章用靶向nifH基因的引物对（Ueda19F-R6，IGK3-DVV和F2-R6）对不同生态的样本做了比较，也是对此流程的验证。
总结
- nifH功能基因同16S一样属于微生物的扩增子研究，但是扩增的区域和研究的目的不同。
- 功能基因研究目前的流程不像16S那么成熟，数据库也各有构建的方法。
- 功能基因多样性研究流程大都参考扩增子分析流程与工具，主要变化体现在氨基酸的校正以及数据库的应用。
参考：
rdpstaff
fungene
FunGene Pipeline
NifMAP
taxadiva
RDP Tutorials
Functional Gene Unsupervised Workflow
FunGene 功能基因数据库
Using RDPTools Output with Phyloseq
refseq介绍
nifh
nirs
amoA_AOA
amoA_AOB
图片备注：
图2 | nifH引物对基于读数比例在环境样品中扩增nifH和非nifH（同源）基因的性能
（A）和OTU（B）。
环境样本包括：
（a）山毛榉森林土壤;
（b）草甸土壤;
（c）根际和（d）Arrhenatherum elatius的根样本;
（e）根际和（f）Oryza sativa的根样本;
（g）沿海，亚北极生物土壤结皮（BSC）;
（h）温和的平衡计分卡;
（i）高山BSC;
（j）半干旱平衡计分卡;
（k）干旱BSC;
来自（l）Great Belt和（m）Roskilde Fjord的河口样本。
关于样品的更多细节可以在“材料和方法”部分中找到。将分类为同源基因（例如bchX，chlL，bchL和parA）的读数或OTU概括为“非nifH”。
图6 | 基于使用不同引物对的定量PCR，环境样品中每ng DNA标准误差的平均nifH拷贝数。 环境样本包括：（a）山毛榉森林土壤; （b）草地
泥; （c）根际和（d）Arrhenatherum elatius的根样本;（e）中
根际和（f）Oryza sativa的根样本; （g）沿海，亚北极
生物土壤结皮（BSC）; （h）温和的平衡计分卡; （i）高山BSC; （j）半干旱
BSC; （k）干旱BSC; 来自（l）Great Belt和（m）Roskilde的河口样品
峡湾。
校正拷贝数以排除使用从具有特异性引物对的样品的扩增子测序获得的信息共同扩增的非nifH基因。
图3 | 基于BLASTP搜索不同引物对和环境样品的nifH序列的分类学分类。 热图的颜色
表示每个样本的分类学类别的平均相对丰度。 环境样本包括：（a）山毛榉森林土壤; （b）草甸土壤; （C）
根际和（d）Arrhenatherum elatius的根样本; （e）根际和（f）Oryza sativa的根样本; （g）沿海，亚北极生物土壤结皮（BSC）; （H）
温带BSC; （i）高山BSC; （j）半干旱平衡计分卡; （k）干旱BSC; 来自（l）Great Belt和（m）Roskilde Fjord的河口样本。
