for a file full of reads?
My goal is to split a sequence at a specific site into two separate sequences. Searching for the site should be a bit fuzzy due to sequencing-pipeline (basecalling on MinION) error.
Example:
Assume a sequence as below. X, Y, Q and Z are sequence nucleotides not necessary for understanding the problem but are useful for demonstration purposes.
XXXXXXXXXXXXXXYYYYYACTCATAQQQQQQQQQZZZZZZZZZZZZZ
|-----|
I would like to find site ACTCATA (with fuzzy matching) and split the sequence into
XXXXXXXXXXXXXXYYYYY
and
QQQQQQQQQZZZZZZZZZZZZZ
with optionally discarding the matched sequence.
Bonus points if this is done on fastq files where data on quality of reads is also split into new strings.
This could probably be accomplished the pedestrian way in biopython but was wondering if I missed a tool that does what I describe above.
1 answer
$ cat read.fq
@seq
XXXXXXXXXXXXXXYYYYYACTCATAQQQQQQQQQZZZZZZZZZZZZZ
+
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
$ cat read.fq | seqkit locate -p ACTCATA
seqID patternName pattern strand start end matched
seq ACTCATA ACTCATA + 20 26 ACTCATA
$ cat read.fq | seqkit subseq -r 1:19
@seq
XXXXXXXXXXXXXXYYYYY
+
GGGGGGGGGGGGGGGGGGG
$ cat read.fq | seqkit subseq -r 27:-1
@seq
QQQQQQQQQZZZZZZZZZZZZZ
+
GGGGGGGGGGGGGGGGGGGGGG
Writing a script is the better way for more than one reads :P
Doesn't this expect an exact match?
regular expression (default) and motif containing degenerate bases like N(-d) are supported:
http://bioinf.shenwei.me/seqkit/usage/#locate
Log in to answer this question.
Are you looking for adapter sequences? If so: https://github.com/rrwick/Porechop
Thank you @WouterDeCoster. I may end up using this in another part of the pipeline.