Removing adapter/primer sequence from reads while keeping flanking regions

Question

[Bryan]

hi folks,

I have some metagenomics data that was generated using an amplification protocol that integrated a known primer sequence into random sites throughout genomes in the sample. I'd like to selectively remove the primer sequence and keep all flanking regions of the reads, a la deleting the primer sequence, but not trimming the read, which is done separately.

An example:

Raw read

sequence1ADAPTERsequence2

Desired result:

sequence1sequence2

I haven't figured out a way to do this with cutadapt (though would love any advice I've missed) and I know BBDuk can mask the adapter sequence with Ns, but that isn't exactly what I'm looking for either. Any suggestions welcome!

thanks!

[bk11]

If you have input fastq files, you should be able to do using cutadapt-

cutadapt -a ADAPTER -o output.fastq input.fastq

Reference

[GenoMax]

This ADAPTER is in the middle of the read. OP wants to remove that and join the left and right pieces.

[GenoMax]

See if this helps: https://bioinf.shenwei.me/seqkit/usage/#amplicon

No standard scan/trim program is likely going to do this. Generally all sequence 3' of where the initial match is made to adapter is trimmed.

Answer 1

[shenwei356]

seqkit amplicon can only keep either 5' or 3' flanking sequence.

How about

seqkit replace -s -p ADAPTER

Assuming there's only one adapter sequence in each read.

$ echo -ne ">s\nsequence1ADAPTERsequence2\n"
>s
sequence1ADAPTERsequence2
$ echo -ne ">s\nsequence1ADAPTERsequence2\n" \
    | seqkit replace -s -p ADAPTER
>s
sequence1sequence2

[swbarnes2]

Will seqkit replace do the same kind of fuzzy matching that cutadapt will? At a glance, it doesn't look like it.

[shenwei356]

No fuzzy matching, this command only supports regular expression.

[GenoMax]

shenwei356 This should work if the adapter sequence is constant? If OP confirms then we can move this to an answer.

Answer 2

[swbarnes2]

It's kind of clunky, but you could use cutadapt to remove the adapter and everything after, save that to a fastq, use cutadapt to remove the adapter and everything before, save that to a second fastq, then join the fastqs together. (check for matching names, since the two files will probably not be in sync). Not sure what program will join fastqs like that, but you could write a script in something to do that.

Answer 3

[Bryan]

thanks for the speedy replies, all! I ended up just resorting to my roots with sed, which produces nearly identical results to shenwei356 seqkit replace answer. As neither supports degeneracy, I'm not sure that either stands out as more optimal, aside for some efficiency concerns, which I didn't evaluate.

$ echo -ne ">s\nsequence1ADAPTERsequence2\n" \
> | sed 's/ADAPTER//g'
>s
sequence1sequence2

I will note that this sed command performs global substitution, so it will replace multiple ADAPTER instances per read, if found, whereas it looks like shenwei356 answer will remove a single instance, in case that matters for folks referring to this in the future.

[GenoMax]

I moved your and Shenwei's comments to an answer. If you want to accept both to provide closure (green check mark) to this thread that would be great.

While sed is fine when you work with defined format like fastq there may be unexpected glitches with sed since it does not understand fastq format. shenwei356 's tool on the other hand does.

Is seq1ADAPTERseq2ADAPTERseq3 an edge case or do you actually expect to see that kind of construct. Because if not this is an erroneous read and should be thrown away.

[Bryan]

thanks for cleaning this up, GenoMax . It's possible to have a situation like seq1ADAPTERseq2ADAPTERseq3, but that would probably only occur in a pretty low biomass sample for this protocol and would need some extra QC anyway. I accepted shenwei356 answer since seqkit does understand fasta/q as you pointed out. Thanks all!