上一篇我们聊了使用QIIME2 dada2的准备工作，发现大家对这一话题还比较感兴趣，问题比较多，感谢大家的关注~本来想一次性介绍完剩下的所有内容，但是发现qiime2的坑挺多，所以这次我还是先把数据导入和切除引物等内容的一些细节问题详细介绍一下（建议用电脑进行阅读）。
在读这篇文章之前，我假设大家已经看过我之前的微文，所以很多地方并没有过多解释。（没读过的朋友可以点击以下微文的链接进行阅读）
1 初学者如何深入解读16S rDNA扩增子测序数据，从而选择自己的分析步骤
2 技术贴 | 16S专题 |基于QIIME2 dada2插件的16S扩增子测序数据的分析流程详解（上）
Dada2方法要求的输入测序数据必须是已经经过区分样本（类似qiime1裂库）的测序数据，这个测序数据也必须是带质量信息的。对于双端测序数据，不能提前拼接，思维比较灵活的人，可能会考虑先拼接，再把拼接好的测序数据当做单端测序数据进行分析，但是这就违背了dada2的假设，dada2假设随着测序长度的增长，测序质量会稍微增加，后急剧下降，而拼接好的序列，测序质量应该是先降低后增加的。在数据量比较小的情况下，Dada2的输入数据最好提前将barcodes和引物这些不想要的adapters切除，虽然dada2也能粗略切除barcodes和引物，但是做不到cutadapt插件那样精细。
对于公司不同的返回数据，我们可能需要做不同的预处理，经过导入并区分样本成为QIIME2带语义类型的qza格式文件。在此之前，统一介绍一下metadata文件，后面就不解释了，这个文件类似以往提的map文件（公司一般会提供这个文件），除样品ID（categorical）、引物（categorical，非必要，如果是双端测序，正向和反向引物可以放一列，可用英文逗号隔开）、barcodes信息（categorical，barcodes一般在正向序列引物前面，如果正向、反向序列都有barcodes，qiime2中，我目前还没有发现能够处理的方法，barcodes放在反向序列引物前面还勉强能处理）之外，还可以放入样本的分组信息（categorical）或者环境因子（numeric）。图1是metadata例子，其中橘色箭头表示制表符（用记事本打开，制表符不会显示，类似空格），制表符在记事本里可用tab键输入，也可以在Excel里面预先编辑好，另存为“文本文件（制表符分隔）”。第一行是变量名称，“#”不可省略，第二行是变量类型，“#q2:types”不要变更，一二行顺序不可颠倒。
图1 metadata文件
1 导入数据和区分样本
1.1 EMP未区分样本、带质量信息的单端测序数据
按照EMP的标准，你的测序数据应该是phred33的质量体系，你的所有样本的序列是完全混合在一起的，barcodes另外放一个样本，你应该有三个文件:1）metadata文件，这里假设名称为“metadata.tsv”；2）序列主体fastq文件，这里假设文件名称为sequences.fastq.gz，没有barcodes（还记得之前教大家怎么判断barcodes和引物在不在序列里吗？所有样本测序数据都在一个fastq文件里，如果是fastq的“.gz”格式，不需要提前解压；3）barcodes fastq文件，这里假设这个文件名为“barcodes.fastq.gz”，这个文件和序列主体fastq文件的序列条数应该是一样的，只不过这个文件的序列只包括了barcodes，如图2，其中“NNNNNNNNNNNN”表示对应序列没有barcode，这样对应的序列在分析过程中会被剔除掉。
图2 barcodes文件
步骤和代码如下：
先新建一个文件夹，例如我们可以在共享文件夹（名为weishengtai，见上篇）里创建一个命名为emp-single-end-sequences的文件夹，将barcodes.fastq.gz和sequences.fastq.gz移动到这一文件夹下，导入的时候只需要指定文件夹，而不用指定文件。
导入序列数据代码如下（注意，windows下代码直接粘贴到Linux命令行，可能会出错，因为两个系统换行符不一样，最好一行一行手动输入，或者全部放到一行中，删除“\”）：
cd /media/sf_weishengtai/;
qiime tools import \
--type EMPSingleEndSequences \
--input-path emp-single-end-sequences \
--output-path emp-single-end-sequences.qza
--input-path 指定导入的文件所在目录（不用指定文件）；--type指定导入文件的语义类型，查看所有可导入的语义类型可使用命令qiime tools import --show-importable-types；--output-path指定输出文件的名字，文件默认放在了当前目录（/media/sf_weishengtai/）下，也可以自己用它指定其他路径。
区分序列所属样本代码如下：
qiime demux emp-single \
--i-seqs emp-single-end-sequences.qza \
--m-barcodes-file sample-metadata.tsv \
--m-barcodes-column BarcodeSequence \
--o-per-sample-sequences demux.qza
其中，sample-metadata.tsv应该放在在当前目录下（/media/sf_weishengtai/），或者用--m-barcodes-file指定路径；--m-barcodes-column指定metadata中，barcode所在列的列名（变量名）；另外，到这里大家应该看出规律了，一般qiime命令--i开头的指定输入文件，--i-seqs 指定的是上一步输出的文件；--o开头的指定输出文件名称，默认保存在当前目录下。这一步得到的输出文件demux.qza就可以作为dada2的输入文件了（也可以选择先切除引物）。
1.2 EMP未区分样本、带质量信息的双端测序数据
按照EMP的标准，你的测序数据应该是phred33的质量体系，你的所有样本的序列是完全混合在一起的，正向和反向序列分开放，barcodes最好是正向的，barcodes另外放一个样本，你有4个文件，比上面多了一个：1）metadata文件，假设名字为sample-metadata.tsv，至少要包含样品ID和barcodes列；2）正向序列fastq文件，没有barcodes，所有样本序列混合在一个文件，假设名字为forward.fastq.gz；3）反向序列文件fastq文件，没有barcodes，假设名字为reverse.fastq.gz；4）barcodes序列fastq文件，假设名字为barcodes.fastq.gz。
先新建一个文件夹，例如我们可以在共享文件夹（名为weishengtai，见上篇）里创建一个命名为emp-paired-end-sequences的文件夹，将barcodes.fastq.gz、forward.fastq.gz、reverse.fastq.gz都移动到emp-paired-end-sequences文件夹下，导入的时候只需要指定文件夹，而不用指定文件。
导入命令如下（先切换当前目录到emp-paired-end-sequences所在目录）：
cd /media/sf_weishengtai/; \
qiime tools import \
--type EMPPairedEndSequences \
--input-path emp-paired-end-sequences \
--output-path emp-paired-end-sequences.qza
其中，--input-path指定barcodes.fastq.gz、forward.fastq.gz、reverse.fastq.gz所在的文件夹名（如果这个文件夹不在当前目录下，需要指定绝对路径）。
区分序列所属样本代码如下：
qiime demux emp-paired \
--m-barcodes-file sample-metadata.tsv \
--m-barcodes-column BarcodeSequence \
--i-seqs emp-paired-end-sequences.qza \
--o-per-sample-sequences demux
其中，--i-seqs指定的是上一步命令的输出文件，如果barcodes是在反向引物上的，要另外添加参数--p-rev-comp-mapping-barcodes，如果正向、反向都有barcodes，抱歉，没有提供处理这种情况的参数。输出文件demux.qza（可用--o-per-sample-sequences参数指定名称）可以直接用于dada2，或者先切除引物。
1.3 区分样本的带质量信息的双端/单端测序数据
国内目前大多数公司给的测序数据形式都是这个，每个样本两个fastq文件，一个放正向序列，一个放反向序列（单端测序只有正向）。要是每个样本两个fasta文件（序列的质量信息已被删除）呢？不好意思，qiime2处理不了这种数据，当然你也可以写个脚本（如python），把序列的质量信息都填充成高质量值（既然已经经过质控，生成fasta，那么我们就可以假设这些序列的质量值都很高了），伪装成fastq文件，这样也能用下面的步骤导入，不过fasta最好还是用qiime1分析。要是双端测序数据一个样本只有一个文件（正向序列和反向序列没有分开放）呢？那你也可以写个脚本（如python）根据序列头的“1”和“2”来把正向序列和反向序列拆开，再用下面的步骤导入。这里，我们只提供大多数情况（每个样本两个fastq文件）的处理步骤：
首先你需要建立一个manifest文件，这个文件公司不会给你的，你需要自己手动书写，不过python使用比较熟练的，也可以用代码快速完成。这个文件的形式如图3所示：
图3 manifest文件
注意，manifest文件使用逗号分隔的，假设名称为manifest.csv。当你用Excel编辑好的时候，记得另存为“文本文件（逗号分隔）”，“#”开头的行可有可无，会被忽略，因为#后面都是注释信息。第一行是各变量的名字，不同于metadata，开头不能写“#”，否则就被忽略了。这个文件相当于一个csv文件，用Excel打开如图4所示，第一列就是样本ID；第二列就是序列文件的绝对路径，注意不是windows下的绝对路径，而是你虚拟电脑里的绝对路径，例如，你把sample1_R1.fastq.gz放在了共享文件夹下（weishengtai），那么sample1_R1.fastq.gz的绝对路径应该是“/media/sf_weishengtai/sample1_R1.fastq.gz”；第三列是序列的方向，不清楚的解压fastq后打开，看序列头是“1”还是“2”，参考之前的微文，单端测序数据都是“forward”。
图4 manifest文件用Excel打开
对于双端测序数据，写好这个manifest文件以后，你就可以用下面的命令导入了：
qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path manifest.csv \
--output-path paired-end-demux.qza \
--source-format PairedEndFastqManifestPhred33
对于单端测序数据，用以下命令：
qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path manifest.csv \
--output-path single-end-demux.qza \
--source-format SingleEndFastqManifestPhred33
其中--input-path指定manifest.csv的路径，这里假设manifest.csv在当前目录（还记得之前教大家怎么改当前目录吗）下，这是唯一的输入文件，记得不要移动序列的位置就行；--source-format指定文件格式，注意质量体系（还记得之前教大家怎么推测质量体系吗？大多数情况是phred33）。这步输出的文件（single-end-demux.qza/ paired-end-demux.qza）直接可以用于dada2（或者先切一下引物）。
1.4 不区分样本的带barcodes、带质量信息的测序数据
其实目前国内很多公司给的应该是区分样本（不同样本的序列放在不同文件里）的带barcodes和引物的fastq测序数据，这种情况你可以用1.3方法导入，再用dada2 --ptrim-trim-left参数大致切一下barcodes和引物；或者你可以选择稍微精确一点但是费事一点的方法，先把所有样本的序列合并到一个文件（如果双端，正向、反向要分开），再用下面的方法完成数据导入和区分样本。
这个类型的数据，应该是不区分样本的，所有样本测序数据放在一个文件，没有切除barcodes和引物（如果是双端测序数据，则正向、反向分开，只有正向序列fastq里有barcodes）。
对于单端数据，先导入：
qiime tools import \
--type MultiplexedSingleEndBarcodeInSequence \
--input-path forward.fastq.gz \
--output-path multiplexed-seqs.qza
注意类型和1.1不一样，是MultiplexedSingleEndBarcodeInSequence，--input-path指定你单端数据的路径，这里假设在当前目录（所以没必要输入绝对路径）。接下来用cutadapt插件区分样本：
qiime cutadapt demux-single \
--i-seqs multiplexed-seqs.qza \
--m-barcodes-file metadata.tsv \
--m-barcodes-column Barcode \
--o-per-sample-sequences demultiplexed-seqs.qza \
--o-untrimmed-sequences untrimmed.qza
--i-seqs指定的是上一步的输出，--m-barcodes-column指定的是metadata.tsv中，barcodes在所在列的列名。
对于双端数据，与上面类似，不过双端数据导入时，--type指定的类型应该是MultiplexedPairedEndBarcodeInSequence，--input-path指定的是forward.fastq.gz和reverse.fastq.gz所在的目录（而不是文件，而且区别于1.2，没有barcodes.fastq）。
很遗憾的一点是，qiime cutadapt demux-paired目前只支持barcodes存在正向序列的情况，不支持barcodes在反向序列，也不支持正向、反向序列都有barcodes的情况，期待它的更新~如果你恰好幸运，barcodes只存在正向序列，可以用qiime cutadapt demux-paired --help命令查看一下怎么使用参数。下面是一个例子
qiime cutadapt demux-paired \
--i-seqs multiplexed-paired-seqs.qza \
--m-forward-barcodes-file metadata.tsv \
--m-forward-barcodes-column Barcode \
--o-per-sample-sequences demultiplexed-seqs.qza \
--o-untrimmed-sequences untrimmed.qza
这里介绍的qiime cutadapt demux方法输出的结果文件（--o-per-sample-sequences指定的文件）也是可以直接用于dada2的，或者先切一下引物。
2 切除引物
上面的1.1~1.4所介绍方法导入的并区分样本的数据，可以先切除引物或者其他不想要的片段，再用作dada2的输入文件。切除引物用的是qiime2 cutadapt插件。可用--help参数查看帮助。
切除单端测序引物的一个例子：
qiime cutadapt trim-single \
--i-demultiplexed-sequences DemuxSeq.qza \
--p-front CCTACGGGNGGCWGCAG \
--o-trimmed-sequences trimmed-seqs.qza
注意这里的引物是直接书写的（在一个试验里，所有样本的引物应该是一样的），而不是用metadata指定的，所以metadata里其实可以不用放引物。
切除双端测序数据引物的例子：
qiime cutadapt trim-paired \
--i-demultiplexed-sequences demux.qza \
--p-front-f CCTACGGGNGGCWGCAG \
--p-front-r ADAPTER2SEQUENCE \
--o-trimmed-sequences trimmed-seqs.qza
其中，--p-front-f指定正向引物，--p-front-r指定反向引物，注意不要用--p-adapter-r或者--p-adapter-f参数，它们匹配的是3’末端的序列。
本文由Bayegy原创，由董小橙、江舜尧编辑。
联系客服
