1
$\begingroup$

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.

Image of FASTA file

$\endgroup$
1
  • 2
    $\begingroup$ Please don't post images of text. Instead, copy/paste the text into your question and use the formatting tools to format it as code. Also, what operating system are you using? We need to know so we can know what tools and what versions of the tools will be available. $\endgroup$
    – terdon
    Nov 5, 2022 at 13:31

4 Answers 4

4
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ FWIW with GNU awk you could do awk 'BEGIN{RS=">"} {ORS=RT} !/Sequence unavailable/' file instead of piping to sed. $\endgroup$
    – Ed Morton
    Jan 11, 2023 at 23:40
3
$\begingroup$

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.

$\endgroup$
0
2
$\begingroup$

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])

$\endgroup$
4
  • $\begingroup$ Note that this will require loading the entire file into memory which might be an issue with large files. $\endgroup$
    – terdon
    Nov 5, 2022 at 13:00
  • $\begingroup$ I just processed a 22Gig file with a slightly more complicated command (removing loci). It worked fine. $\endgroup$
    – M__
    Nov 5, 2022 at 22:26
  • 1
    $\begingroup$ Thank you very much for your kind response. This is work fro me. $\endgroup$
    – Dr Sud
    Nov 6, 2022 at 9:42
  • $\begingroup$ Hi @DrSud great! I use it because its versatile (use it for different stuff). Please do remember to mark an answer "accepted" if you find it helpful and you can also now upvote all the answers given here. $\endgroup$
    – M__
    Nov 6, 2022 at 13:30
2
$\begingroup$

One way using :

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 '/^>/ { 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.

$\endgroup$
4
  • $\begingroup$ There's a bug in the sed code, sed: 1: "/^>/ { N }; /Sequence u ...": extra characters at the end of N command. OSX Darwin $\endgroup$
    – M__
    Nov 5, 2022 at 18:22
  • $\begingroup$ @M__ ah... my mistake. Should be good now I think. Thanks for testing it on OSX $\endgroup$
    – Steve
    Nov 5, 2022 at 21:54
  • $\begingroup$ Thanks but same error but as I'm not feeling motivated to test on Linux ... lets say it an OSX thing, most work on Linux. The awk works fine. $\endgroup$
    – M__
    Nov 5, 2022 at 22:29
  • $\begingroup$ @M__ I don't think I read your error message properly before. Sorry. Hopefully the above will work now.. $\endgroup$
    – Steve
    Nov 6, 2022 at 14:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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