? keeping only matching sequence of reads

[ashipde]

With seqkit grep we can filter reads based on whether their sequence matches a pattern. But is there a way to use seqkit on fastq input to also trim the reads so that they only contain the matching sequence?

E.g. for pattern "T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T":

Input:

@NB501936:87:HJCMMBGX5:1:22305:4332:2522 1:N:0:CGCTCATT+CTTCGCCT
CCAAGGGAAGGAGAGGTGAGATGATTCGCCTATTTGAGGTTAGCATCCGCGGTTGTGAGTATACACACGGTAGGT
+
AAAAAEEEEEEEEAEEEEEEEEEEEEEEEAAEEEAEEEEEEEEEEEEEEEEEAEE6EEAEEEEEEEEEEEEEE<E

Output:

@NB501936:87:HJCMMBGX5:1:22305:4332:2522 1:N:0:CGCTCATT+CTTCGCCT
TTCGCCTATTTGAGGTTAGCATCCGCGGTTGTGAGTATACACACGGT
+
EEEEEAAEEEAEEEEEEEEEEEEEEEEEAEE6EEAEEEEEEEEEEEE

[shenwei356]

$ cat t.fq | seqkit grep -s -r -P -p "T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T"
@NB501936:87:HJCMMBGX5:1:22305:4332:2522 1:N:0:CGCTCATT+CTTCGCCT
CCAAGGGAAGGAGAGGTGAGATGATTCGCCTATTTGAGGTTAGCATCCGCGGTTGTGAGTATACACACGGTAGGT
+
AAAAAEEEEEEEEAEEEEEEEEEEEEEEEAAEEEAEEEEEEEEEEEEEEEEEAEE6EEAEEEEEEEEEEEEEE<E

$ cat t.fq | seqkit grep -s -r -P -p "T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T" \
    | seqkit replace -s -p ".+(T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T).+" -r '$1'
[ERRO] editing FASTQ is not supported

$ cat t.fq | seqkit grep -s -r -P -p "T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T" \
    | seqkit fq2fa \
    | seqkit replace -s -p ".+(T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T).+" -r '$1'
>NB501936:87:HJCMMBGX5:1:22305:4332:2522 1:N:0:CGCTCATT+CTTCGCCT
TTCGCCTATTTGAGGTTAGCATCCGCGGTTGTGAGTATACACACGGT

[ashipde]

Thank you for the prompt response! Unfortunately I also need the phred scores.

[ashipde]

Awk alternative that serves my purpose:

cat t.fq | paste - - - - \
  | awk -F'\t' '{
      if(match($2, "T.{9}TGAGGTTAG.{9}TTGTGAGTA.{9}T")){
        $2 = substr($2, RSTART, RLENGTH);
        $4 = substr($4, RSTART, RLENGTH);
      }
      print $1 "\t" $2 "\t" $3 "\t" $4
    }' \
  | tr '\t' '\n'