Removing particular line along with Fasta header in fasta file

Question

I have a fasta file of 3'UTR sequences of the genes. Some of the gene do not have a 3'UTR sequence and reported with sequence unavailable but there is FASTA header.

I need to delete or remove all the "sequence unavailable" lines along with their header, i.e. all blank fasta entries.

Please suggest a command with awk or sed command.

Answer 1

You can use awk here, by setting the input (RS) and output (ORS) record separator, this is what defines a "line" in awk, to >. Then, tell it to print records that do not match Sequence unavailable. For instance, using this input file:

$ cat file
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>baz
Sequence unavailable
>bar1
AGAGATAGAGAGATATA
>asld
Sequence unavailable
>asd
ACTTC

You can do:

$ awk  'BEGIN{RS=ORS=">"}!/Sequence unavailable/' file 
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC
>

That will add a trailing >, so you need one more step to remove it:

$ awk  'BEGIN{RS=ORS=">"}!/Sequence unavailable/' file | sed 's/>$//'
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC

---

Alternatively, using the FastaToTbl and TblToFasta scripts I have posted here, you can do:

$ FastaToTbl file | grep -v 'Sequence unavailable' | TblToFasta 
>foo 
AGAGTAGATAGAAGATA
>bar 
AGAGATAGAGAGATATA
>bar1 
AGAGATAGAGAGATATA
>asd 
ACTTC

Both approaches have the advantage of not needing to read the entire file to memory so can be used on files of any size.

Answer 2

Try this, if you have GNU grep (the default on Linux systems):

grep --no-group-separator -i -B 1 '^[ACGTN]\+$' input.fasta > filtered.fasta

Explanation: find all lines that are composed of only A/C/G/T/N (case insensitive search), and show that line plus the preceding line. Don't put any breaks between matches that are not consecutive.

This should work on both multi-line and single line fasta files.

Answer 3

If there are no line breaks within the sequence data (the format I use) then the following will work.

perl -0777 -nle 'while (/^(>.+\n^[AGCT-]+)\n/gm){print $1,"\n"}' myfile.fa

---

Test

>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>baz
Sequence unavailable

Output

>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA

---

However If there are line breaks within the sequence data, please use the following ...

perl -0777 -nle 'while (/^(>.+\n[AGCT-]{10}.+\n[^>]+)/gm){print $1}' myfile.fa

Then pipe to a file,

perl -0777 -nle 'while (/^(>.+\n[AGCT-]{10}.+\n[^>]+)/gm){print $1}' myfile.fa > myfile_minus_blanks.fa

This works due to a multi-line m

Test ...

>foo
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
AGAGATAGAGAGATATA
>baz
Sequence unavailable

Out ...

>foo
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
AGAGATAGAGAGATATA

Note If you've mixed bases, degenerate nucleotides,

AGCT- should be replaced with AGCTRYWSMKNBDHV (omitting - [not needed])

Answer 4

One way using sed:

sed -e '/^>/ { N; }; /Sequence unavailable/ { d; }' file.fa

For each header line, append the next line of input into the pattern space. Then, if the pattern space contains "Sequence unavailable" delete the pattern space and start a new cycle.

---

Or using awk:

awk '/^>/ { f=$0; next } /^Sequence unavailable/ { next } f { print f; f="" }1' file.fa

For each header line, use a variable to store the current record. We can use next here to force awk to immediately stop processing the current record and move on to the next record, which is just the next line, since, by default, awk's record separator is a newline character. Then, stop processing the current record and move on to the next record if that line begins with "Sequence unavailable". Finally, print the header line if it exists and unset the variable that we used to store it if it does. Once we've printed the header, we can then just print the current record. We can do this by appending 1 to our command string. The 1 on the end just returns true, and by doing so, it performs the default action, which is to print the current record.