# Fungal ITS analysis tutorial **URL:** https://forum.qiime2.org/t/fungal-its-analysis-tutorial/7351 **Category:** Tutorials **Created:** [December 9, 2018, 2:56pm UTC](https://forum.qiime2.org/t/fungal-its-analysis-tutorial/7351 "2018-12-09T14:56:30Z") **Posts on this page:** 1 **Showing post:** 1
### Author: ![Nicholas\_Bokulich](https://forum.qiime2.org/user_avatar/forum.qiime2.org/nicholas_bokulich/32/19937_2.png) [@Nicholas\_Bokulich](https://forum.qiime2.org/u/Nicholas_Bokulich) #### Post date: [December 9, 2018, 2:56pm UTC](https://forum.qiime2.org/t/fungal-its-analysis-tutorial/7351/1 "2018-12-09T14:56:30Z")
# Fungal ITS analysis, mock communities, and more fun > _NOTE: This tutorial was written in a QIIME 2 2018.11 environment. It is not guaranteed to work with earlier or later versions of QIIME 2. This tutorial was compiled as a working exercise for a [QIIME 2 workshop](https://workshops.qiime2.org/microbiome-bioinformatics-qiime-2-2018-12-12/) in December 2018, and does not represent the only possible fungal ITS workflow with QIIME 2, or even a benchmarked protocol recommendation. See other tutorials, e.g., the [q2-itsxpress tutorial](https://forum.qiime2.org/t/q2-itsxpress-a-tutorial-on-a-qiime-2-plugin-to-trim-its-sequences/5780) for other fungal ITS analysis options in QIIME 2._ This tutorial details different modalities of analysis for fungal internal transcribed spacer (ITS) data with QIIME 2. Some of these concepts may generalize to other non-16S rRNA marker genes, such as: 1. issues with sequence length heterogeneity 2. fitting your own classifier using a custom database (must be in a QIIME 2-compatible format!) Additionally, this tutorial will give a nice overview of using QIIME 2 and mock communities (e.g., from [mockrobiota](http://mockrobiota.caporasolab.us)) to benchmark analysis pipelines. These concepts are useful for any type of marker gene, be it ITS, 16S rRNA, or anything else under the sun. :sunshine: Let's begin. ```auto mkdir fungi-tutorial cd fungi-tutorial ``` # Fit a classifier for the [UNITE database](https://unite.ut.ee/) ```auto wget https://files.plutof.ut.ee/doi/0A/0B/0A0B25526F599E87A1E8D7C612D23AF7205F0239978CBD9C491767A0C1D237CC.zip ``` ```auto unzip 0A0B25526F599E87A1E8D7C612D23AF7205F0239978CBD9C491767A0C1D237CC.zip ``` We need to scrub lowercase characters from the sequences to avoid errors downstream. This is because this version of the UNITE database contains some lowercase characters — this is not something that needs to be done to all databases prior to importing to QIIME 2. Some sequences also have pesky blank spaces at the ends of the lines, so we will remove blank spaces from these lines (credit to @angrybee for spotting this and adding the `sed` command below to fix this): ```auto awk '/^>/ {print($0)}; /^[^>]/ {print(toupper($0))}' developer/sh_refs_qiime_ver7_99_01.12.2017_dev.fasta | tr -d ' ' > developer/sh_refs_qiime_ver7_99_01.12.2017_dev_uppercase.fasta ``` ```auto qiime tools import \ --type FeatureData[Sequence] \ --input-path developer/sh_refs_qiime_ver7_99_01.12.2017_dev_uppercase.fasta \ --output-path unite-ver7-99-seqs-01.12.2017.qza ``` ```auto qiime tools import \ --type FeatureData[Taxonomy] \ --input-path developer/sh_taxonomy_qiime_ver7_99_01.12.2017_dev.txt \ --output-path unite-ver7-99-tax-01.12.2017.qza \ --input-format HeaderlessTSVTaxonomyFormat ``` Finally, let's just train the classifier. This will take some time! So we can just download one that I've made for you below 😉 ```auto qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads unite-ver7-99-seqs-01.12.2017.qza \ --i-reference-taxonomy unite-ver7-99-tax-01.12.2017.qza \ --o-classifier unite-ver7-99-classifier-01.12.2017.qza ``` As promised, you can just download this instead: ```auto wget https://s3-us-west-2.amazonaws.com/qiime2-data/community-contributions-data-resources/2018.11-unite-classifier/unite-ver7-99-classifier-01.12.2017.qza ``` # Download an ITS mock community from [mockrobiota](http://msystems.asm.org/content/1/5/e00062-16) We will use a fungal ITS mock community published by [Taylor et al., 2016](https://aem.asm.org/content/82/24/7217). First download the sequence data and metadata from mockrobiota: ```auto wget -O "mock-25-sample-metadata.tsv" https://raw.githubusercontent.com/caporaso-lab/mockrobiota/master/data/mock-25/sample-metadata.tsv wget https://s3-us-west-2.amazonaws.com/mockrobiota/latest/mock-25/mock-forward-read.fastq.gz wget https://s3-us-west-2.amazonaws.com/mockrobiota/latest/mock-25/mock-reverse-read.fastq.gz ``` Create a [fastq manifest file](https://docs.qiime2.org/2018.11/tutorials/importing/#fastq-manifest-formats) for importing these data into QIIME 2, `fastqmanifest.csv`: ```auto echo "sample-id,absolute-filepath,direction" > fastqmanifest.csv echo "Mock.1,$PWD/mock-forward-read.fastq.gz,forward" >> fastqmanifest.csv echo "Mock.1,$PWD/mock-reverse-read.fastq.gz,reverse" >> fastqmanifest.csv ``` ```auto qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path fastqmanifest.csv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33 ``` ```auto qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv ``` ## To trim or not to trim One issue with ITS (and other marker genes with vast length variability) is _readthrough_, which occurs when read lengths are longer than the amplicon itself! The polymerase will read through the amplicon, the primer, the barcode, and on into the adapter sequence. This is non-biological DNA that will cause major issues downstream, e.g., with sequence classification. So we want to trim primers from either end of the sequence to eliminate read-through issues. Enter cutadapt. Note that we trim the forward primer and the reverse complement of the reverse primer from the forward reads (the forward primers have already been trimmed in the raw reads, but we will demonstrate forward + reverse trimming here since attempting to trim the forward read will not hurt). We trim the reverse primer and reverse complement of the forward primer from the reverse reads. ```auto qiime cutadapt trim-paired \ --i-demultiplexed-sequences demux.qza \ --p-adapter-f AYTTAAGCATATCAATAAGCGGAGGCT \ --p-front-f AACTTTYRRCAAYGGATCWCT \ --p-adapter-r AGWGATCCRTTGYYRAAAGTT \ --p-front-r AGCCTCCGCTTATTGATATGCTTAART \ --o-trimmed-sequences demux-trimmed.qza ``` ```auto qiime demux summarize \ --i-data demux-trimmed.qza \ --o-visualization demux-trimmed.qzv ``` Compare the `demux-trimmed.qzv` summary to the `demux.qzv` summary that we made before trimming. It doesn't look like the sequences were trimmed, does it? According to [Taylor et al. 2016](https://aem.asm.org/content/82/24/7217) (the source of the primers and mock community used here), the predicted amplicon length range of the 5.8S-Fun + ITS4-Fun primer set is 267-511 bp (see Figure 1). So with the read lengths we are using here (250 nt) we should not see any read-through. ## Denoising Even though we have paired-end reads, the read quality in the mock-25 dataset is not good enough to support merging the expected amplicon length (See the "interactive quality plots" tab in `demux.qzv` to visualize how quality drops off in both forward and reverse reads.) We can use `qiime dada2 denoise-paired` to confirm that we lose the majority of our reads at the merging stage (I have already confirmed this). The best solution in this case is to proceed with `denoise-single` (using only the forward reads). ```auto qiime dada2 denoise-single \ --i-demultiplexed-seqs demux-trimmed.qza \ --p-trim-left 13 \ --p-trunc-len 160 \ --o-representative-sequences dada2-single-end-rep-seqs.qza \ --o-table dada2-single-end-table.qza \ --o-denoising-stats dada2-single-end-stats.qza ``` ```auto qiime metadata tabulate \ --m-input-file dada2-single-end-stats.qza \ --o-visualization dada2-single-end-stats.qzv ``` In your own data you should follow this process to determine whether you have enough sequence length to proceed with pairend-end sequences, or if you must proceed with single-end reads due to insufficient read overlap (as we have done above): 1. Determine the expected amplicon length for your primers 2. Check the amplicon length and assess the quality of your reads (use the `demux.qzv` or `demux-trimmed.qzv` summaries) 3. Are your (trimmed) reads long enough to overlap? dada2 expects a minimum overlap of 12 nt, so your reads are long enough if: length(forward read) + length(reverse read) ≥ 12 nt. 4. If your answer to #3 is "definitely not!", proceed with single-end reads as we have shown here. If your answer to #3 is "maybe?", give paired-end denoising a shot but pay very close attention to how many reads are being lost after merging (check out the stats file that is output by dada2). # Taxonomic classification Now we will use the taxonomy classifier that we trained above to predict the taxonomic assignment for each of our denoised sequence variants. ```auto qiime feature-classifier classify-sklearn \ --i-classifier unite-ver7-99-classifier-01.12.2017.qza \ --i-reads dada2-single-end-rep-seqs.qza \ --o-classification taxonomy-single-end.qza ``` ```auto qiime taxa barplot \ --i-table dada2-single-end-table.qza \ --i-taxonomy taxonomy-single-end.qza \ --m-metadata-file mock-25-sample-metadata.tsv \ --o-visualization taxa-bar-plots.qzv ``` [taxa-bar-plots.qzv](https://cdck-file-uploads-global.s3.dualstack.us-west-2.amazonaws.com/flex002/uploads/qiime21/original/2X/d/d979f52221189cb5c434a5291dc33ad69b414ad5.qzv) (328.4 KB) # Other analyses Now that you have a feature table and taxonomy classifications, you can typically proceed with analysis of fungal ITS (and other marker-gene data) in the same way that you would 16S rRNA gene data. With one major exception. Although this marker provides more accurate species identification than other fungal marker genes, alignments of ITS sequences are not useful for informing evolutionary distances among distantly related species (but can yield good alignments for closely related sequences). So you have two options for diversity and other potentially phylogenetic analyses: 1. do not use phylogenetic methods with your ITS data. Instead, use non-phylogenetic methods diversity analyses, such as Jaccard (qualitative) or Bray Curtis (quantitative) for estimating beta diversity. 2. Use [q2-ghost-tree](https://forum.qiime2.org/t/q2-ghost-tree-plugin-community-tutorial-for-creating-hybrid-gene-phylogenetic-trees/6139) to infer phylogeny via grafting ITS sequences. This will enable use of phylogenetic methods, including UniFrac. **Edit 2024: ghosttree is no longer maintained, so is no longer a solution for estimating phylogenies from ITS sequences** Many other tutorials (e.g., the [overview for beginners](https://docs.qiime2.org/2018.11/tutorials/overview/#fun-with-feature-tables)) outline some of the possible downstream analyses possible with _any_ feature table (in theory). So we will not revisit these options here. Let's talk about something else instead. # Evaluate accuracy _NOTE: If you are following this tutorial to analyze your own data (of real samples), stop here. The following section is specific for this tutorial and other mock communities. This is not a meaningful analysis of real samples._ We just analyzed a mock community, meaning that we know the _expected composition_ of that sample's taxonomic abundances and sequences. The expected composition would be precisely what our inputs were, converted to relative abundances. Taxonomic names are mapped to the nearest match in our reference sequence database of choice (the same one we used for taxonomic classification) — something that mockrobiota already does for us. So we can use mock communities like this to evaluate he accuracy of different bioinformatics methods or pipelines. Let's compare the methods we used to see how well they did! First we will grab the expected composition of this mock community: ```auto wget https://raw.githubusercontent.com/caporaso-lab/mockrobiota/master/data/mock-25/unite/7-1/99-otus/expected-taxonomy.tsv ``` That file was annotated with an older version of the UNITE database, so let's fix some of the annotations (to match the version we are using) before proceeding: ```auto sed 's/s __Tylospora_asterophora/s__ unidentified/' expected-taxonomy.tsv | sed 's/p __Zygomycota;c__ Mortierellomycotina_cls_Incertae_sedis/p __Mortierellomycota;c__ Mortierellomycetes/' | sed 's/c __Chytridiomycetes/c__ Spizellomycetes/' | sed 's/g __Coprinopsis/g__ Coprinopsis;s __unidentified/' | sed 's/g__ Tricholoma/g __Tricholoma;s__ Tricholoma_vaccinum/'> expected-taxonomy-mod.tsv ``` Convert to biom and import ```auto biom convert \ -i expected-taxonomy-mod.tsv \ -o expected-taxonomy.biom \ --table-type="OTU table" \ --to-json qiime tools import \ --type FeatureTable[RelativeFrequency] \ --input-path expected-taxonomy.biom \ --input-format BIOMV100Format \ --output-path expected-taxonomy.qza ``` Collapse our features on taxonomy and convert our feature table of observations to relative frequency: ```auto qiime taxa collapse \ --i-table dada2-single-end-table.qza \ --i-taxonomy taxonomy-single-end.qza \ --p-level 7 \ --o-collapsed-table dada2-single-end-table-collapsed.qza qiime feature-table relative-frequency \ --i-table dada2-single-end-table-collapsed.qza \ --o-relative-frequency-table dada2-single-end-table-relative.qza ``` Evaluate accuracy: ```auto qiime quality-control evaluate-composition \ --i-expected-features expected-taxonomy.qza \ --i-observed-features dada2-single-end-table-relative.qza \ --o-visualization dada2-single-end-evaluation.qzv ``` [dada2-single-end-evaluation.qzv](https://cdck-file-uploads-global.s3.dualstack.us-west-2.amazonaws.com/flex002/uploads/qiime21/original/2X/0/0703a427f6ea43daaeda9957de3c150fb060b117.qzv) (357.4 KB) What does it all mean? You can read more about TAR and TDR scores [here](https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-018-0470-z). Other metrics should be more familiar to you. What this is showing: We do pretty well! At least as far as TDR is concerned. TDR is high at species level, indicating a very low number of false negatives (see list at bottom). TAR increases with sequencing depth (as does # of observed species/expected), indicating an increasing number of false positives. What are we missing? Which species are being misclassified and why? BONUS: we can use this evaluation to compare multiple methods. How does dada2 single end with cutadapt trimming compare to no trimming or other methods? --- _[View the full topic](https://forum.qiime2.org/t/fungal-its-analysis-tutorial/7351)._