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 at 12:35
  • why the -w (=--word-regexp) switch for grep? – törzsmókus Mar 24 at 15:16
8

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 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. – Stéphane Chazelas Mar 24 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 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 at 12:49
  • @rowboat thanls, I added it to the answer. – schrodigerscatcuriosity Mar 23 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 \+). – Stéphane Chazelas Mar 24 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.