數(shù)據(jù)來源：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]
科學問題
腸道微生物的正常變化對宿主健康的影響媒殉。在這個教程中作者截取了部分數(shù)據(jù)毙死，試圖回答的問題是斷奶的前10天快速增加的體重是否對微生物組結(jié)構(gòu)有影響以及與斷奶140天到150天微生物組的比較陕习。
數(shù)據(jù)集
數(shù)據(jù)集中包含41個文件脉让，對于10個時間節(jié)點的雌鼠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
第一步是合并每個樣品雙端測序的文件和所有樣品的數(shù)據(jù)乒省。這一步的命令是用make.contigs巧颈。該命令提取reads和相應(yīng)的得分，將反向read反向互補后一起拼接成contigs.
make.contigs(file=stability.files)
代碼結(jié)果會顯示每個樣品的序列數(shù)目袖扛，同時每個樣品文件分別會生成stability.trim.contigs.fasta和stability.contigs.groups文件砸泛。
可用summary.seqs查看結(jié)果：
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)
計算每個唯一序列出現(xiàn)在每個組中的次數(shù)：
count.seqs(name=stability.trim.contigs.good.names, group=stability.contigs.good.groups)
統(tǒng)計序列信息：
summary.seqs(count=stability.trim.contigs.good.count_table)
將序列比對到參考序列晾嘶，這里用的是silva.bacteria.fasta妓雾，SILVA 有50000列娶吞，包含了18S rRNA和16S rRNA序列。作者認為這個數(shù)據(jù)庫比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)
發(fā)現(xiàn)有些序列的起始和終止位置不正常妒蛇，所以繼續(xù)處理序列：
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)
去除質(zhì)體中的序列烈钞，比如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中進行。
在進行后續(xù)的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的序列數(shù)目：
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)
這里生成的樹是用的全部的序列毯欣，而且節(jié)點的名稱和OTUs里的不同，導致不能在phyloseq里使用臭脓。
