PacBio 16S full-length amplicone analysis pipeline using QIIME2

Software requirement

  • QIIME2
  • parallel
  • seqtk

QIIME2 analysis pipeline

0. activate conda environment

conda activate qiime2-2019.10

1. flip sequences to the same direction

mkdir raw_data_rc/
mkdir raw_data_cat/
parallel -j 8 'seqtk seq -r {} > raw_data_rc/{/.}_rc.fastq' ::: rawdata/*.fastq
parallel -j 8 --link 'cat {1} {2} > raw_data_cat/{1/.}_cat.fastq' ::: rawdata/*.fastq ::: raw_data_rc/*_rc.fastq

2. trim primers

mkdir trimmed_reads/

parallel -j 8 'cutadapt -g AGRGTTYGATYMTGGCTCAG...AAGTCGTAACAAGGTARCY \
                --discard-untrimmed --no-indels -j 1 -m 1200 -M 1800 \
                -o trimmed_reads/{/.}_trim.fastq {}' \
               ::: raw_data_cat/*_cat.fastq

3. import sequences as a qiime2 artifact

mkdir reads_qza

qiime tools import \
    --type SampleData[SequencesWithQuality] \
    --input-path PacBioCCSmanifest \
    --output-path reads_qza/trimmed_reads.qza \
    --input-format SingleEndFastqManifestPhred33

4. DADA2 denoising

qiime dada2 denoise-single \
   --i-demultiplexed-seqs reads_qza/trimmed_reads.qza \
   --p-trunc-len 0 \
   --p-n-threads 8 \
   --p-n-reads-learn 100000 \
   --output-dir dada2_output

5. convert artifact to readable format

qiime tools export --input-path dada2_output/denoising_stats.qza --output-path dada2_output
qiime tools export --input-path dada2_output/table.qza --output-path dada2_output_exported
qiime tools export --input-path dada2_output/representative_sequences.qza --output-path dada2_output_exported
biom convert -i dada2_output_exported/feature-table.biom -o dada2_output_exported/feature-table.txt --to-tsv

6. assign taxonomy

time qiime feature-classifier classify-consensus-vsearch \
    --i-query dada2_output/representative_sequences.qza \
    --i-reference-reads /home/share/Q2_Silva99_DB/Silva99_16SRepSeq.qza \
    --i-reference-taxonomy /home/share/Q2_Silva99_DB/Silva99_16STaxonomy.qza \
    --p-maxaccepts 1 --p-strand "plus" \
    --p-threads 10 \
    --verbose \
    --output-dir taxa
qiime tools export --input-path taxa/classification.qza --output-path taxa

7. annotate feature table

sed -i -e '1 s/Feature/#Feature/' -e '1 s/Taxon/taxonomy/' taxa/taxonomy.tsv
biom add-metadata -i dada2_output_exported/feature-table.biom -o dada2_output_exported/feature-table_w_tax.biom --observation-metadata-fp taxa/taxonomy.tsv --sc-separated taxonomy
biom convert -i dada2_output_exported/feature-table_w_tax.biom -o dada2_output_exported/feature-table_w_tax.txt --to-tsv --header-key taxonomy

Diversity analysis & visualization

MARco R package

Please cite: Liu, P.-Y. poyuliu/MARco: MARco: Microbiome Analysis RcodeDB (Version v1.0). Zenodo. Available online: http://doi.org/10.5281/zenodo.4589898 (accessed on Day Month Year).

Barplot
Alpha diversity
Beta diversity