Removing reads based on a pattern in the sequence name

Question

[fufuyou]

Hi, I have some reads like as:

>seq_5150639_x4
GGAACGAGATCGTCTGCAGTTGGC
>seq_5150619_x40
AACCGCCTGTAGAAATGCATGATT

X4 or X40 indicate how many reads same with the read. I want to remove lower than 10. How can I remove them? Thanks, Fuyou

[venu]

Do you mean you want to remove reads whose ids containing x10, x9 ... ?

[Pierre Lindenbaum]

the anwser is awk, again.

please, validate your previous questions: sequence head change! ; Add a sequence name! ; ...

Answer 1

[Alex Reynolds]

Assuming the FASTA input is single-line:

1. Parse each FASTA record header with awk and get the last part of the record header, split by the underscore (_)
2. Strip the first character from the last part and cast the rest to an integer
3. If the integer is greater than or equal to the threshold (10), set a flag value to a true-like value and print the header
4. Print the FASTA record sequence if the flag is set to a true-like value

For example:

$ awk '{ \
    if ($0 ~ /^>/) { \
      n = split($0, a, "_");  \
      r = int(substr(a[n], 2)); \
      f = 0; \
      if (r >= 10) { \
        f = 1; \
        print $0; \
      } \
    } \
    else if (f) { \
      print $0; \
    } \
  }' input.fa > output.fa

[fufuyou]

Thanks, It works well. Fuyou