PacBio ITS workflow - some questions

[Rob_DNA]

Hello,

I've use QIIME2 quite a lot for Illumina paired end data, but now it is the first time I'll use it with PacBio full ITS data. The primers used are ITS9mun (GTACACACCGCCCGTCG) and ITS4ngsUni (CGCCTSCSCTTANTDATATGC). FYI: For this specific experiment, I need to trim the ITS2 region, even though I have full ITS reads.

I'm currently analyzing full ITS pacbio CCS reads and I currently use the following workflow:

1. First I orient the reads outside of QIIME2 so all start with the forward primer. (I do this because I use a ITSxpress step which trims the primers, and thus dada2 denoise-ccs cannot be used because --p-front cannot be specified - see later)

for f in *.fastq
do
    cutadapt \
        -j 80 \
        --revcomp \
        --action=none \
        -e 0.12 \
        -g '^GTACACACCGCCCGTCG' \
        --discard-untrimmed \
        -o "${f%.fastq}.oriented.fastq" \
        "$f"
done

2. Then of course data import into QIIME2.

qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path manifest.tsv \
--output-path raw.orient.qza \
--input-format SingleEndFastqManifestPhred33V2

3. Then I use itsxpress to extract the fungal ITS2 region:

 qiime itsxpress trim-single \
        --i-per-sample-sequences raw.orient.qza \
        --p-region ITS2 \
        --p-cluster-id 1.0 \
        --p-taxa F \
        --p-threads 80 \
        --o-trimmed itstrimmed.raw.orient.qza

4. Then I use dada2 denoise-single for denoising:

qiime dada2 denoise-single \
--i-demultiplexed-seqs itstrimmed.raw.orient.qza\
--p-trunc-len 0 \
--p-n-threads 80 \
--o-table table-DADA2.qza \
--o-representative-sequences seqs-DADA2.qza \
--o-denoising-stats stats-DADA2.qza

steps afterwards:

5. clustering on 97% ID using vsearch

6. taxonomic identification using the pre-trained UNITE database ( Releases · colinbrislawn/unite-train · GitHub ), but my questions are regarding these initial 4 steps.

I know there are several (old) posts about PacBio workflows in QIIME2, but PacBio does not seem to be fully supported yet, and hence I still have some questions:

1. Did I miss anything in my workflow?

2. In this post,Adam Rivers uses

reformat.sh in=./pacbio/pacbio_samp1.fastq.gz out=pacbio_samp1_reformat.fastq.gz mincalledquality=2 maxcalledquality=41 qin=33

to reformat the PacBio qual scores into Illumina scores, but is that actually necessary? itsxpressand dada2 seem to function fine without it.

3. Because itsxpress trims primers from the reads, I cannot use dada2 denoise-ccsbecause that requires a --p-frontparameter which is not present after itsxpress (see e.g., this post - I assume there is no fix possible as this is inherent to the mixed orientation of pacbio CCS reads?) . Is it wise to use dada2 denoise-single on these oriented and itsxpress processed reads? In the post from Adam Rivers, he uses denoise-single but that post is from 2018 so I'm not sure how representative that is currently. As I extracted the ITS2 region, the reads are relatively short (I guess 150-250 bp).

The dada2 stats look good (pilot pipeline with 3 samples):

[I'm using the qiime2 Amplicon 2024.10 distribution installed with conda].

[Adam_Rivers]

I think the Dada2-Qiime2 plug-in added support for denoise-ccs in version QIIME 2 2022.2. We have never supported that flag but it probably makes sense to do that. I think that --denoise-ccs models the error rate much more accurately for PacBio data. I can try to add this feature so you don't need to do these work arounds. Would you be willing to share sequence data that I can use for testing? I'm not sure I completely understand your specific use case. You are trying to extract just the ITS2 region even though you have long reads? I would think it would be more common to extract the full its or both ITS1 and ITS2 linked sequences, but I have not done this.

[Rob_DNA]

Hi Adam,

thanks for the quick reply!

@Adam_Rivers

I can try to add this feature so you don't need to do these work arounds.

What feature are you specifically referring to here? So basically I'd like to use denoise-ccs while trimming the ITS2 region from the reads. Another option would be to append the forward primer artifically to each oriented, ITS2 trimmed read and use that forward primer sequence as --p-frontin dada2 denoise-ccsbut i do not know if that would result in some issues within dada2?

@Adam_Rivers

I'm not sure I completely understand your specific use case. You are trying to extract just the ITS2 region even though you have long reads? I would think it would be more common to extract the full its or both ITS1 and ITS2 linked sequences, but I have not done this.

I now realize this might look a bit odd, but indeed for this specific experiment I only want the ITS2 region. I basically want to integrate this PacBio dataset with Illumina ITS2 data, hence I want to trim the ITS2 region for (some) comparability.

@Adam_Rivers

Would you be willing to share sequence data that I can use for testing?

yes, I'm using publicly available ITS sequencing data. I've added two fastq PacBio CCS files as test data. Is that enough? (file size wise I could not add more, but I could share more with an external link)

two_pacbio_fastq.7z (7.7 MB)

[Adam_Rivers]

I just added support for qiime itsxpress trim-single --trim-css option in Version 2.2.0. This release also improves documentation, Docker and Bioconda builds, and has some performance enhancements. ITSxpress can now orient ccs reads using Vsearch and attach dummy primer sequences to the ends that Dada2 can use to verify read orientation. Let me know if you encounter any issues.