I am using the Miseq SOP to analyze 18S sequences and using Silva 132 for alignment and classification. I’m looking at the v9 region using Earth Microbiome primers 1391f and EukB, expected amplicon length around 260+/-50. After creating a reference alignment to analyze my mock community, it looks like the sequence itself, including the primers, is about 174 bp (close to the published length of 168bp for S. cerevisiae, Stoek et al 2010, cited on Earth Microbiome website in reference to these primers).
I am trying to use an S. cerevisiae sequence from that reference alignment and follow the mothur “customize you reference alignment” protocol for 16S data, as listed below. However, the summary table I’m getting is all zeroes. What am I missing here?
Sequence in Scerevisiae_v9.fasta:
GTACACACCGCCCGTCGCTAGTACCGATTGAATGGCTTAGTGAGGCCTCAGGATCTGCTTAGAGAAGGGGGCAACTCCATCTCAGAGCGGAGAATTTGGACAAACTTGGTCATTTAGAGGAACTAAAAGTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCA
Forward GTACACACCGCCCGTC
Reverse Complement of Reverse Primer GTAGGTGAACCTGCAGAAGGATCA
my guess is that your reference alignment for some reason does not match your sequences and they are therefore filtered out. The way I go with creating such a custom reference alignment would be like this:
Take one typical FASTA sequence from your “YourSampleName.trim.contigs.good.fasta” file and paste it into the SILVA aligner (https://www.arb-silva.de/aligner/). This way you are sure that your custom alignment will fit to your amplicons.
The result of that alignment you should “view in wasabi” (mark the finished job and the button appears), which gives you the start & end positions of your sequence within the vast 50K nucleotide SILVA alignment. For your Scerevisiae_v9.fasta sequence this would for example be start: 41790 / end: 43275.
Use the command: mothur>pcr.seqs(fasta=silva.nr_v132.align, start=41790, end=43275, taxonomy=silva.nr_v132.tax, keepdots=F) and you will have (hopefully) your custom alignment and its corresponding tax file.