数据来源：Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Applied and Environmental Microbiology. 79(17):5112-20.
教程参考：MiSeq SOP[https://www.mothur.org/wiki/MiSeq_SOP]

科学问题

肠道微生物的正常变化对宿主健康的影响。在这个教程中作者截取了部分数据，试图回答的问题是断奶的前10天快速增加的体重是否对微生物组结构有影响以及与断奶140天到150天微生物组的比较。

数据集

数据集中包含41个文件，对于10个时间节点的雌鼠3和1个对照。F3D0表示femal 3 on day 0 (断奶的时间)。

首先进入MiSeq_SOP文件，然后在mothur中输入：

    set.dir(input=你的文件夹位置）
    make.file(inputdir=., type=fastq, prefix=stability)

1、Reducing sequencing and PCR errors

第一步是合并每个样品双端测序的文件和所有样品的数据。这一步的命令是用make.contigs。该命令提取reads和相应的得分，将反向read反向互补后一起拼接成contigs.
make.contigs(file=stability.files)
代码结果会显示每个样品的序列数目，同时每个样品文件分别会生成stability.trim.contigs.fasta和stability.contigs.groups文件。
可用summary.seqs查看结果：
summary.seqs(fasta=stability.trim.contigs.fasta)
双端测序的长度一般只有250bp，所以进一步用screen.seqs去除那些长度明显不对的序列：
screen.seqs(fasta=stability.trim.contigs.fasta, group=stability.contigs.groups, maxambig=0, maxlength=275)

2、Processing improved sequences

去除重复序列：
unique.seqs(fasta=stability.trim.contigs.good.fasta)
计算每个唯一序列出现在每个组中的次数：
count.seqs(name=stability.trim.contigs.good.names, group=stability.contigs.good.groups)
统计序列信息：
summary.seqs(count=stability.trim.contigs.good.count_table)
将序列比对到参考序列，这里用的是silva.bacteria.fasta，SILVA 有50000列，包含了18S rRNA和16S rRNA序列。作者认为这个数据库比greengenes要好。
pcr.seqs(fasta=silva.bacteria.fasta, start=11894, end=25319, keepdots=F, processors=8)
重新对文件命名：
rename.file(input=silva.bacteria.pcr.fasta, new=silva.v4.fasta)
比对到处理好的文件：
align.seqs(fasta=stability.trim.contigs.good.unique.fasta, reference=silva.v4.fasta)
查看比对信息：
summary.seqs(fasta=stability.trim.contigs.good.unique.align, count=stability.trim.contigs.good.count_table)
发现有些序列的起始和终止位置不正常，所以继续处理序列：
screen.seqs(fasta=stability.trim.contigs.good.unique.align, count=stability.trim.contigs.good.count_table, summary=stability.trim.contigs.good.unique.summary, start=1968, end=11550, maxhomop=8)
保留了1968到11550之间的序列，同时保证相同的多聚核苷酸不超过8。
去除序列首尾的gap characters：
filter.seqs(fasta=stability.trim.contigs.good.unique.good.align, vertical=T, trump=.)
重新去重：
unique.seqs(fasta=stability.trim.contigs.good.unique.good.filter.fasta, count=stability.trim.contigs.good.good.count_table)
接下来对序列进行预聚类，这样进一步降噪：
pre.cluster(fasta=stability.trim.contigs.good.unique.good.filter.unique.fasta, count=stability.tirm.contigs.good.unique.good.filter.count_table, diffs=2)
这里会对序列进行分组，然后排序，在两两比较，如果序列中有2个碱基不同，则会将其合并。
去除嵌合体(chimeras)：

chimera.vsearch(fasta=stability.trim.contigs.good.unique.precluster.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.count_table, dereplicate=t)
remove.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta, accnos=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.accnos)

去除质体中的序列，比如mitochondria等：

classify.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denono.vsearch.pick.count_table, reference=trainset9_032012.pds.fasta, taxonomy=trainset9_032012.pds.tax, cutoff=80)
remove.lineage(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.count_table, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.taxonomy, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota)

此步也可以在phyloseq中进行。
在进行后续的OTU和phylotypes注释时需要将序列文件中对照的序列去除：
remove.groups(count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.pick.count_table, fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.fasta, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.pick.taxonomy, groups=Mock)

OTUs

cluster.split(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.fasta, count=stability.trim.contigs.good.unique.good.filter.unque.precluster.denovo.vsearch.pick.pick.pick.count_table, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.pick.pick.taxonomy, splitmethod=classify, taxlevel=4, cutoff=0.03)

每个group中每个OTU的序列数目：
make.shared(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.opti_mcc.list, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.pick.pick.count_table, label=0.03)

Phylotypes

phylotype(taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.pick.pick.taxonomy)
上面会列出从属到界，如果只要属水平的文件：
make.shared(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.pick.pick.tx.list, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.pick.pick.count_table, label=1)
然后对这些OUTs分类到phylotypes：
classify.otu(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.pick.pick.tx.list, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.pick.pick.count_table, taxonomy=stability.trim.contigs.good.filter.unique.precluster.pick.pds.wang.pick.pick.taxonomy.label=1)

Phylogenetic

dist.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.fasta, output=lt)
clearcut(phylip=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.phylip.dist)

这里生成的树是用的全部的序列，而且节点的名称和OTUs里的不同，导致不能在phyloseq里使用。