Qiime2+Origin绘制稀释曲线

稀释曲线(Rarefaction Curve)也称稀疏曲线，一般在微生物组研究中多用于评估测序量或样本量的饱和情况。利用dada2去噪获得的table文件，计算随机抽取n个（n小于测得reads序列总数）reads时出现的ASV数量，然后根据一组n值（一般为一组小于总序列数的等差数列）与其相对应的ASV数量作出曲线。

作图所需数据包含每个样本ASV的特征表，样本元数据。关注“环微分析”公众号，后台回复“稀释曲线”即可获取示例数据。

一、QIIME2实操

(1) 打开QIIME2，执行命令挂载共享文件夹qiimeshare；

sudo mount -t vboxsf qiimeshare share

​(2) 使用qiime diversity alpha-rarefaction可视化工具绘制稀释曲线，探索α多样性与采样深度的关系。执行命令，查看帮助文档；

qiime diversity alpha-rarefaction --help

二、帮助文档解析

使用说明：qiime diversity alpha-rarefaction

通过计算“min_depth”和“max_depth”之间的稀疏度，生成alpha稀疏曲线。计算中间深度的数目由"steps"参数控制，在每个稀疏深度处计算n个'iterations'，即迭代n次。如果提供了样本元数据，还可以根据元数据列中不同的值对样本进行分组。

输入文件：

--i-table 用于计算稀释曲线的特征表（必需文件）

--i-phylogeny 使用phylogeny align-to-tree-mafft-fasttree生成的有根树文件，系统发育度量（可选文件）

参数:

--p-max-depth 最大稀疏深度，要比最小深度大，从1开始（必需参数）；

--p-metrics 计算多样性指数，可以选择'pielou_e', 'shannon', 'faith_pd', 'ace', 'simpson', 'observed_features', 'chao1', 'goods_coverage'等，更多请看上图。如果没有定义该参数，默认计算observed_features和shannon，如果输入了树文件可以计算faith_pd（可选参数）；

--m-metadata-file 元数据文件（可选参数）；

--p-min-depth 最小稀疏深度，默认为1（可选参数）；

--p-steps 在最小深度和最大深度之间的稀疏深度数，从2开始，默认为10（可选参数）；

--p-iterations 计算每一步稀疏的迭代次数，次数越多，计算量也越大，默认为10（可选参数）；

输出文件：

 --o-visualization 名称（必须参数）

其它选项：

--output-dir 输出文件路径（可选参数）

三、实操小结

Alpha稀释曲线可视化工具在每个采样深度处计算Alpha多样性指数，范围介于最小采样深度--p-min-depth和最大采样深度--p-max-depth提供值之间。在每个采样深度生成10个抽样表，并对表中的所有样本进行alpha多样性指数计算（即迭代次数为10，在每个采样深度计算的稀疏表，通过--p-iterations来控制）。在每个采样深度为每个样本绘制平均多样性值，如果提供样本元数据--m-metadata-file参数，还可以基于元数据对样本进行分组。

四、运行分析

使用qiime diversity alpha-rarefaction在不同深度（在--p-min-depth和--p-max-depth之间）对ASV表进行子采样，并使用两个默认度量observed_features和shannon来计算alpha多样性（--p-metrics）。在每个采样深度，默认计算10个稀疏表以提供误差估计（--p-iterations）。

执行命令：

time qiime diversity alpha-rarefaction \ --i-table dada2-table-paired.qza \ --m-metadata-file sample-metadata.tsv \ --p-max-depth 18000 \ --o-visualization alpha-rarefaction-curves-1.qzv

稀疏分为两个步骤：首先，低于稀疏深度的样本被从特征表中滤掉；然后，对所有剩余样本进行无放回采样以达到指定的测序深度。

无放回抽取

无放回抽取是指每次抽出的一个个体不再放回总体中，下次再抽取时，总数比前一次少一个，每次抽取的概率发生变化，无放回抽取的各次抽取作为一个事件，他们不是相互独立的。例如，无放回抽取事件“无放回地逐个抽取n条序列”与“一次性任取n条序列”的概率相等。

稀疏深度

最大稀疏深度的值应根据dada2-table-paired.qzv文件中呈现的“每个样本的测序量”信息来确定。

一般来说，选择一个在中位数附近的值似乎很好用。如果得到的稀疏图中的线看起来没有变平，可以增加该值；如果由于大于最大采样深度而丢失了许多样本，则减少该值。通过不断调整，最终确定一个合适的值。因此，我们首先选择18000进行尝试。

输出结果：

alpha-rarefaction-curves-1.qzv | 稀疏曲线结果可视化

五、结果分析

① 首先选择Metric为observed_features，样本元数据列为place(不进行分组)。查看结果，如下图：

图表解读：

稀疏图(rarefaction plot)主要用于确定样品的丰度是否已被完全观察或测序。如果图中的线条在沿x轴的某个采样深度处看起来“平坦(level out)”（即斜率接近于零），这表明收集超过该采样深度的附加序列不太可能观测到新特征。如果绘图中的线条没有变平，这可能是因为尚未充分观察样本的丰富度（由于测序的序列太少），或者它可能是在数据中仍然存在许多测序错误（被误认为是新的多样性）。

本数据得到的稀疏曲线如上图，有一条曲线比较短，是因为18000超出了它的样本深度17375（sample8）。当序列深度达到10000以后，曲线趋于平行，说明在测序深度10000取样的数量合理，能保证绝大部分的特征均已被观察到，样本alpha多样性指标达到稳定。

② 样本元数据列调整为group(分组)进行查看。我们一共有9个样本，分成3组，每组3个样本，发现TS组的指数变平之后在16000之后突然上升。如下图：

图表解读：当通过元数据对样本进行分组时，此可视化结果底部的绘图结果非常重要。它说明了当特征表被细化到每个采样深度时，每个组中剩余的样本数量。如果给定的采样深度d大于样本s的总频率(针对样本s获得的序列数)，则不可能计算采样深度d下样本s的多样性。在顶部绘图将不可靠，因为它将计算基于相对少的样本。因此，当通过元数据对样本进行分组时，必须查看底部图表，以确定顶部图表中显示的数据是否可靠的。

我们设定的采样深度为18000，比最小深度17375（sample8）大，因此TS样本数在16000之后突然成2，观测到的特征数突然变大。调整参数，再次运行分析。

③ 增加步长，选择常见多样性指数，再次重新绘制曲线。执行命令：

time qiime diversity alpha-rarefaction \ --i-table dada2-table-paired.qza \ --m-metadata-file sample-metadata.tsv \ --p-metrics 'shannon' 'simpson' 'observed_features' 'chao1' 'goods_coverage' \ --p-max-depth 10000 \ --p-steps 25 \ --o-visualization alpha-rarefaction-curves-2.qzv

输出结果：

alpha-rarefaction-curves-2.qzv | 稀疏曲线结果可视化

解压缩qzv之后所含结果如下图所示：

六、疑问解答

① 有人可能想问，为什么不直接选好参数，设置好步长？

我们可以看到第一次以默认参数运行用了31s,而增加步长和选用多个多样性指数之后用了2min20s,所以我们以默认参数进行分析可以花费更少的时间尝试不同的采样深度，确定最终合适的参数，再进行深入的分析。

② 为什么要增加步长重新运行一次？

根据小编的经验，QIIME2可视化的图可能无法直接在文章中使用，因此需要根据稀释曲线的数据结果在其他软件重新绘制，然而默认10个间隔的抽样结果有时不能保证曲线平滑，因此可以增加间隔数。此外，也可以增加抽样次数。

七、使用origin绘制稀释曲线

根据observed_features.csv文件的数据内容在Origin中绘制曲线。

（1）使用excel打开数据，如下图所示：

（2）任意选中一格>按住Ctrl+A选中全部内容>复制内容>进行转置粘贴>整理成如下格式；

（3）计算出每个样本在各采样深度iter1-10的平均值，新增一个工作表格sheet1整理成如下格式，另存为xlsx类型文件；

（4）打开Origin，导入数据。数据>从文件导入>Excel>选择文件>调整导入sheet1；

（5）按住Ctrl+A选中所有数据>绘图>基础2D>样条图；

（6）调整美化图形，如下图所示：

八、总结——为什么要进行稀释曲线分析？

微生物组样本中的测序深度与群落中的原始生物量没有直接关系，但相对测序深度对观察到的群落具有很大影响。因此，对于大多数多样性度量，必须对数据进行标准化以解决样本之间不均匀的测序深度。目前最佳的方法是使用稀疏，通过二次取样进行标准化。

Alpha多样性通过抽平来计算（抽平可以对数据进行标准化），但是一次抽平有概率（小概率）在一定程度上评估错误的alpha多样性结果。所以现在通过多次抽平计算alpha多样性，并通过求取均值的方式来矫正alpha多样性，体现了在每一个采样深度迭代10次的重要性。

稀释曲线是对单个alpha多样性结果的补充，可以从不同梯度全面地分析和展示结果。基于不同深度或样本量水平上展示alpha多样性，更加有利于对微生物群体多样性的综合评估。

提示：如果你对几种常见的alpha多样性指数的计算公式很熟悉，你会发现alpha多样性指数均与观测到的OTU数目相关。

 这篇推文对你有帮助吗？喜欢这篇文章吗？喜欢就不要错过呀，关注本知乎号查看更多的环境微生物生信分析相关文章。亦可以用微信扫描下方二维码关注“环微分析”微信公众号，小编在里面载入了更加完善的学习资料供广大生信分析研究者爱好者参考学习，也希望读者们发现错误后予以指出，小编愿与诸君共同进步！！！

学习环境微生物分析，关注“环微分析”公众号，持续更新，开源免费，敬请关注！

转载自原创文章：

Qiime2+Origin绘制稀释曲线

最后，再次感谢你阅读本篇文章，真心希望对你有所帮助。感谢！