3

I have a FASTA file that has intentionally some sequences with wrong header (i.e absence of >) and some with good header. The file is well-formatted in the sense that the nucleotidic sequence is in one line.

Example :

2865958
AACTACTACAG
>hCoV-19/2832832
ACTCGGGGGG
28328332
ATTCCCCG
>hCoV-19/2789877
ACTCGGCCC

And I want to only keep the sequence with a correct header (i.e line that starts with >) like this :

>hCoV-19/2832832
ACTCGGGGGG
>hCoV-19/2789877
ACTCGGCCC

I've tried various method for it ( sed, grep, awk ) but no proper result :

awk '/^>/ { ok=index($0,"hCoV")!=0;} {if(ok) print;}' combined_v4.fa > combined_v5.fa

sed -n '/^>.*hCoV/,/^>/ {/^>.*hCoV/p ; /^>/! p}' combined_v4.fa > combined_v5.fa

grep -w ">" -A 1 combined_v4.fa > combined_v5.fa

Do you have an idea how to do that?

2
  • 2
    How about grep -A 1 '^>hCoV-'
    – Mr R
    Mar 23, 2021 at 12:35
  • why the -w (=--word-regexp) switch for grep? Mar 24, 2021 at 15:16

5 Answers 5

9

Tell grep too look for lines starting with >, and include the line following it:

grep -A1 --no-group-separator '^>' combined_v4.fa > combined_v5.fa

In case your version of grep does not support --no-group-separator, try this:

grep -A1 '^>' combined_v4.fa | grep -v '^--$' > combined_v5.fa
1
  • 4
    A shame --no-group-separator hasn't made it to man. It's documented in info grep invoking command-line context , though.
    – FelixJN
    Mar 23, 2021 at 13:15
3

Here's one way with sed:

sed -n '/^>/!d;N;p' file

Note that if the last line of input begins with ">", this will not print it.

With awk:

awk 'prev {print prev ORS $0; prev=""} /^>/ {prev=$0}' file

Wait until the line after ">" to print the two lines. Insead of a regex, index() could be used:

awk 'prev {print prev ORS $0; prev=""} index($0,">")==1 {prev=$0}' file

^> is a regular expression to match a string that begins with ">". index($0,">") returns the starting position of the string ">", which must to be equal to one (...==1) if the record starts with it.

To make sure it also has "hCoV":

sed -n '/^>/!d;/hCoV/!d;N;p' file
awk 'prev {print prev ORS $0; prev=""} /^>/ && /hCoV/ {prev=$0}' file
1
  • Or just sed '/^>/!d;N'. Then if the last line starts with >, whether it's printed will depend on the sed implementation and in the case of GNU sed, whether POSIXLY_CORRECT is in the environment or not. Mar 24, 2021 at 8:08
3

With awk, you could do:

awk '/^>/ {c=2} c-- > 0' file

to print 2 lines when > is seen.

1
  • 2
    you could use c&&c-- so it stops counting when c is 0 and so c can't wrap around if the input file is massive enough for the negative to become positive again. it's an edge condition for sure!
    – Ed Morton
    Mar 24, 2021 at 12:56
2

I propose this, using GNU sed:

# find the multiline pattern ^digits→ACGT
# and delete those lines
$ sed '/^[[:digit:]]\+$/,/^[ACGT]\+$/d' file 
>hCoV-19/2832832
ACTCGGGGGG
>hCoV-19/2789877
ACTCGGCCC
3
  • 2
    \+ in a BRE might be specific to GNU sed.
    – rowboat
    Mar 23, 2021 at 12:49
  • @rowboat thanls, I added it to the answer. Mar 23, 2021 at 12:54
  • 1
    The standard equivalent of GNU's \+ is \{1,\}. Or you could use + with -E (soon standard, and more widely supported than \+). Mar 24, 2021 at 8:10
1

With awk you can try:

awk 'f {print; f=0} /^>/ {print;f=1}' file
>hCoV-19/2832832
ACTCGGGGGG
>hCoV-19/2789877
ACTCGGCCC

  • This excellent post can help you:

"Printing with sed or awk a line following a matching pattern": https://stackoverflow.com/questions/17908555/printing-with-sed-or-awk-a-line-following-a-matching-pattern

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .