In silico extraction of V regions from greengenes2 - AngusBall/The-Wailing-Somnambulant GitHub Wiki
Lets extract the V region from different primers from aligned sequences
Qiime2 DOES NOT WORK
This isn't too hard ngl\
conda activate rachis-qiime2-2026.7
cd Documents/story/Maastricht/Taiga_in_G_Major/data/On\ taxonomy\ calling/Raw\ from\ greengenes2/
Then we need to convert our aligned backbone into different V sections.
qiime feature-classifier extract-reads \
--i-sequences aligned_16s_backbone.qza \
--p-f-primer GTGGGTAA \
--p-r-primer GGCTAAT \
--p-min-length 100 \
--p-max-length 400 \
--p-n-jobs 7 \
--o-reads extracted_v4_backbone.qza
NOTES: do NOT use the reverse complement primer sequence for the reverse primer
Okay So What regions do we actually need?
We have the V123 region from the primers 27F/519R
We have the V3(V4) region from primers like 341/519, 347/806, 341/785, 341/806
We have the V4 region from primers like 530/803, 515/806 and modifications/derivatives
We have the V45 region from primers like 519/926, 515/928
And finally we have the V678 region from primers 926/1392
In total, I'm deciding, based on vibes, and quantity of samples per primer to create 6 models (unless this blows up my computer).
Model 1: V123 region with primers 27F/519R. These are the only primers in this section so easy choice.
Model 2: V3 region with primer 341/519. This is the shortest V3 region.
Model 3: V34 region with primers 341/806. This is because most samples in the V34 region come from this primer set.
Model 4: V4 region with primers 515/806 from the EMP cause duh.
Model 4.5: I will check what read lengths I have for the V4 region since there was some single read samples in this area. I will make some concatenated sequences if this is necessary (see cDEPP paper, they did this)
Model 5: V45 region with primers 515/926, most samples in this region come from this primer set
Model 6: V678 region: using primers 926/1392. Its the only option \
The code is thus:
We need to import the aligned backbone, this is found by unzipping the backbone.sepp file\
qiime tools import \
--input-path GG2_aligned_fullbackbone.fasta \
--output-path GG2_aligned_fullbackbone.qza \
--type 'FeatureData[AlignedSequence]'
Then we can process it\
Model 1: V123 with 27F/519R
qiime feature-classifier extract-reads \
--i-sequences GG2_aligned_fullbackbone.qza \
--p-f-primer AGAGTTTGATCMTGGCTCAG\
--p-r-primer GWATTACCGCGGCKGCTG\
--p-n-jobs 7 \
--o-reads extracted_v123_backbone.qza
Model 2: V3 with 341/519
qiime feature-classifier extract-reads \
--i-sequences GG2_aligned_fullbackbone.qza \
--p-f-primer CCTAYGGGDBGCWSCAG\
--p-r-primer GWATTACCGCGGCKGCTG\
--p-n-jobs 7 \
--o-reads extracted_v3_backbone.qza
Model 3: V34 region with primers 341/806.
qiime feature-classifier extract-reads \
--i-sequences GG2_aligned_fullbackbone.qza \
--p-f-primer CCTAYGGGDBGCWSCAG\
--p-r-primer GGACTACNVGGGTHTCTAAT\
--p-n-jobs 7 \
--o-reads extracted_v34_backbone.qza
Model 4: V4 region with primers 515/806
qiime feature-classifier extract-reads \
--i-sequences GG2_aligned_fullbackbone.qza \
--p-f-primer GTGYCAGCMGCCGCGGTAA\
--p-r-primer GGACTACNVGGGTHTCTAAT\
--p-n-jobs 7 \
--o-reads extracted_v4_backbone.qza
Model 5: V45 region with primers 515/926,
qiime feature-classifier extract-reads \
--i-sequences GG2_aligned_fullbackbone.qza \
--p-f-primer GTGYCAGCMGCCGCGGTAA\
--p-r-primer CCGYCAATTYMTTTRAGTTT\
--p-n-jobs 7 \
--o-reads extracted_v45_backbone.qza
Model 6: V678 region: using primers 926/1392
qiime feature-classifier extract-reads \
--i-sequences GG2_aligned_fullbackbone.qza \
--p-f-primer CCTATCCCCTGTGTGCCTTGGCAGTCTCAGAAACTYAAAKGAATTGRCGG\
--p-r-primer CGTATCGCCTCCCTCGCGCCATCAGACGGGCGGTGTGTRC\
--p-n-jobs 7 \
--o-reads extracted_v678_backbone.qza
SOOOOOO about all that, apparently feature-classifier extract-reads does not work for FeatureData[AlignedSequence]. RIP. Mother is maybe an option. lets try this one more time.
MOTHER
Downloading is really easy. I'm using the windows version.
Generic example. This is my first idea as to run. input the aligned backbone, the oligo primers and keep all the sequences.
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v123.oligos, nomatch=keep)
Where GG2_aligned_fullbackbone.fasta is the alligned greengenes full 16S region backbone and v3.oligos are the primers you use in the format here. While I'm adding links here the info on pcr.seqs
nomatch= keep keeps sequences that don't match to the primer set. In theory these sequences would therefore not be amplified by the PCR but I don't want to remove them on a poor justification such as that.
HOWEVER, most sequences did not end up being a match. When I ran: (default on nomatch is to remove)
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v123.oligos)
Only 60 909 sequences passed! Approx. 1/6 of the sequences. This is likely because I required exact primer matching to the er well primers. This is no worries however. Because from these 60 909 sequences its very easy to determine the extent of the V3 region from these primers as compared to the total alignment. Quite literally start counting dots from the left to get your start cut value and then count till the end of the sequence to get your end value.
er by count I mean use notepad++
in the first example from the v123 region, every sequence started 20 nucleotides from the left and ended at 555 nucleotides from the left
We can then use direct measurements to get a final alignment with all the sequences
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, start=20, end=555, nomatch=keep)
Except this doesn't work for some reason!!! So many sequences still get busted up. This is due to leading and trailing dots or something?? You can see others with this issue here. The recommendation is to use chop.seqs. sounds good enough to me.
chop.seqs(fasta=GG2_aligned_fullbackbone.fasta, keep=front, countgaps=true, numbases=555)
chop.seqs(fasta=GG2_aligned_fullbackbone.chop.fasta, keep=back, countgaps=true, numbases=535)
then rename GG2_aligned_fullbackbone.chop.chop.fasta to something meaningfull.
PS. I spot checked this a couple times using the original PCR primer sequences. PS2. why in the second one you use numbases=535?? We are keeping bases from the back up to 20 bases from the front so its 555-20=535+1.
Sweet, now as they say: find what you liek and then do it on purpose
V3 with 341/519
30k sequences passed primer, 387-555
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v3.oligos)
chop.seqs(fasta=GG2_aligned_fullbackbone.fasta, keep=front, countgaps=true, numbases=555)
chop.seqs(fasta=GG2_aligned_fullbackbone.chop.fasta, keep=back, countgaps=true, numbases=168)
V34 region with primers 341/806. 99% sequences passed, 387-823
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v34.oligos)
chop.seqs(fasta=GG2_aligned_fullbackbone.fasta, keep=front, countgaps=true, numbases=823)
chop.seqs(fasta=GG2_aligned_fullbackbone.chop.fasta, keep=back, countgaps=true, numbases=436)
V4 region with primers 515/806 90% of seqs, 570-823
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v4.oligos)
chop.seqs(fasta=GG2_aligned_fullbackbone.fasta, keep=front, countgaps=true, numbases=823)
chop.seqs(fasta=GG2_aligned_fullbackbone.chop.fasta, keep=back, countgaps=true, numbases=254)
V45 region with primers 515/926, 570-959
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v45.oligos)
chop.seqs(fasta=GG2_aligned_fullbackbone.fasta, keep=front, countgaps=true, numbases=959)
chop.seqs(fasta=GG2_aligned_fullbackbone.chop.fasta, keep=back, countgaps=true, numbases=390)
V678 region: using primers 926/1392 AS is this primer had no matches, gonna allow pdiffs and rdiffs to be one
pcr.seqs(fasta=GG2_aligned_fullbackbone.fasta, oligos=v678.oligos, pdiffs =1, rdiffs=1)
chop.seqs(fasta=GG2_aligned_fullbackbone.fasta, keep=front, countgaps=true, numbases=XXX)
chop.seqs(fasta=GG2_aligned_fullbackbone.chop.fasta, keep=back, countgaps=true, numbases=XXX)