Demultiplexing PacBio sequences

[uloeber]

Hi all,
I have the following problem:
I have a dataset containing circurlarized and linear DNA. It was amplified with two pairs of primers using inverse PCR. I'd like to dissect the data by the following scheme:
Primer1a----------Primer1b
Primer2a----------Primer2b
Primer1a---------- or ----------Primer1a
Primer1b---------- or ----------Primer1b
Primer2a---------- or ----------Primer2a
Primer2b---------- or ----------Primer2b

I tried to use cutadapt for this issue, but it seems that for the paired primers I only get the sequences with e.g. 1a at 5', 1b at 3' and have to specify it the other way around.
cutadapt -g ^P1 -a P2\$ --trimmed-only -e 0.05 --no-trim >out.fasta in.fasta
I need an approach where I can specify a certain error rate as well as the possibility to search for example within the first and last 40 bp for my 20 bp primer.
Do you have any suggestions???
Thanks in advance!

[Brian Bushnell]

You can use BBDuk to look for this kind of thing:

bbduk.sh in=reads.fq out=unmatched.fq outm=matched.fq ref=primer1a.fa k=20 restrictleft=40 edist=3 -da

That will allow an edit distance of 3 and only look for matching kmers in the first 40 bp of the reads. You can subsequently trim the reads in another pass:

bbduk.sh in=matched.fq out=trimmed.fq ref=primer1a.fa k=20 restrictleft=40 edist=3 ktrim=l -da

For the adapters on the right end, you'd use "ktrim=r" and "restrictright=40".

[uloeber]

Hi Brian,
Thanks for your reply! I forgot to mention, I want to keep the reads untrimmed. And would your software care about the "both primers present" issue? Otherwise I have again to check every possibility separately for every sample. That's what I try to avoid. But again, thank you very much! I appreciate your help. Bbnorm is one of my favorites.
Best wishes,
Ulrike

[Brian Bushnell]

Hi Ulrike,

The first command will not do any trimming, only filtering. And no - unfortunately, BBDuk will not handle both ends at once, so you would need 2 passes. However, if you have a lot of different primers, you can put them all in a file and demultiplex with Seal like this:

seal.sh in=reads.fq pattern=out_%.fq ref=primers.fa restrictleft=40 k=20 edist=3 -da

That would produce one output file per primer sequence. Might save some time.