新一批的数据选的是v5-7区,之前silva v3-4区的训练库不合适了,所以要重新训练适合自己数据的库。
1.Obtaining and importing reference data sets
1.1 首先,第一步先下载好所需的文件:序列文件和taxonomy文件。
于是我得到了两个文件:
silva-138-99-seqs.qza和
silva-138-99-tax.qza
2.Extract reference reads
已有研究表明,当朴素贝叶斯分类器只对已测序的目标序列区域进行训练时,16S rRNA基因序列的分类准确率有所提高(Werner et al., 2012)。这可能不一定适用于其他标记基因(例如真菌ITS分类注释)。我们从 Moving Pictures tutorial 中得知,我们试图分类的序列reads是120个碱基的单端reads,用515F/806R引物对16S rRNA基因序列进行扩增。我们通过从参考数据库中提取匹配引物对的reads,然后将结果切片到120个碱基。
我们将silva-138-99-seqs.qza输入进去,并找到相应的引物(测序公司报告中有提供):
AACMGGATTAGATACCCKG
ACGTCATCCCCACCTTCC
qiime feature-classifier extract-reads \
--i-sequences silva-138-99-seqs.qza \
--p-f-primer AACMGGATTAGATACCCKG \
--p-r-primer ACGTCATCCCCACCTTCC \
--p-trunc-len 235 \
--p-min-length 100 \
--p-max-length 500 \
--o-reads ref-seqs.qza
经过了半个小时左右,得到了文件ref-seqs.qza。
其中的——p-trunc-len参数:应该只用于在查询序列被裁剪到相同长度或更短的情况下裁剪引用序列。成功连接的两端序列的长度通常是可变的。没有被截断到特定长度的单端读也可能在长度上是可变的。对于成对端读和未修剪的单端读的分类,我们建议训练一个分类器,在适当的引物位点上提取的序列,但没有修剪。
对于参数的选择,应参考demux.qzv文件中的交互质量图。
参考资料:
我的demux.qzv文件中是双端,左端为248,右端为235以后,质量下降很快,所以选择len 235
–p-trunc-len 235
接着,是min-length和max-length的确定,
–p-min-length 100
–p-max-length 400
如何来确定这个数字?
官网说的是:
上面的示例命令使用最小长度和最大长度参数来排除使用这些引物的模拟扩增物远远超出预期的长度分布。这样的放大可能是非目标命中,应该被排除。如果您将此命令用于您自己的使用,请确保选择适合标记基因的设置,而不是这里使用的设置。min-length参数应用于trim-left和trunc-len参数之后,max-length应用于_before_,因此一定要设置适当的设置,以防止有效序列被过滤掉。
请教过其他人后得知:
最小的应该是依据trim-left这个参数来确定,即碱基数剪去那些比较低质量的开头的序列,我的就是235-0,因为开头序列质量还不错。不过哪怕前面质量真的很好,我们一般也要留出几个bp的空间,所以200是ok的。
最长的话一般是双端序列加起来多一点,我的就是235+248=483,多留出一些空间,就选择500吧。
–p-min-length 200
–p-max-length 500 \
问:这是针对这批次的数据的demux来训练的,那要是下次换了新的数据,还得新换库?
答:不是的,因为同一家公司做出来的数据质量不会差很多,所以你的测序深度大概都在220到250之间,不会差很多,所以你范围留大一点,基本后面都能适用。所以,最后选择:
–p-min-length 100
–p-max-length 500 \
3.Train the classifier
现在,我们可以使用刚才创建的引用读取和分类,像下面这样训练Naive Bayes分类器。
qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads ref-seqs.qza \
--i-reference-taxonomy silva-138-99-tax.qza \
--o-classifier silvav57classifier.qza
经过大概一个半小时,得到了新的训练器silvav57classifier.qza,哇噻噻,真棒!
4.Test the classifier
最后,通过对自己的数据进行分类并将分类任务可视化来验证分类器的有效性。
qiime feature-classifier classify-sklearn \
--i-classifier silvav57classifier.qza \
--i-reads rep-seqs-dada2.qza \
--o-classification taxonomy.qza
qiime metadata tabulate \
--m-input-file taxonomy.qza \
--o-visualization taxonomy.qzv