Question: Removing identifier with "unavailable sequence" from FASTA file
1
gravatar for seta
2.7 years ago by
seta920
Sweden
seta920 wrote:

Dear all,

I have several FASTA nucleotide file, containing millions identifiers and their sequence. Among them, there are some identifier like this, ">AT1G01340|AT1G01340.2 Sequence unavailable". Would you please let me know how I can remove them? looking forward to hearing your helpful commands. thanks

blast rna-seq sequence genome • 1.7k views
ADD COMMENTlink modified 2.7 years ago by Manvendra Singh1.9k • written 2.7 years ago by seta920

wouldn't this work for you ? It would remove identifiers where sequence is not available

#!/usr/bin/perl
use strict;
use warnings;

$/="\n>";
while (<>) {
 s/>//g;
 my ($id, $seq) = split (/\n/, $_);
 print ">$_" if ((length $seq) > 1);
}

save it with .pl extension and run it on your file.

hth

ADD REPLYlink modified 2.7 years ago by Ram12k • written 2.7 years ago by Manvendra Singh1.9k

thanks, but it dose not work at all.
 

ADD REPLYlink written 2.7 years ago by seta920
1

Oh yes, I had thought that "Sequence unavailable" means that you have empty lines.

Anyways David and Siva have already given you answer, neverthless, you can edit above script's last line 

instead 

(length $seq) > 1

write

($seq) !~ "Sequence"

 

ADD REPLYlink written 2.7 years ago by Manvendra Singh1.9k

It should work. Did you run it with your file? (perl script.pl input.fasta)

ADD REPLYlink written 2.7 years ago by Ram12k
print ">$_" if ((length $seq) > 1);

Would this not print everything? "Sequence unavailable" would also satisfy the condition.

ADD REPLYlink written 2.7 years ago by Siva1.5k

Oh yes, I had thought that "Sequence unavailable" means that he has empty lines.

ADD REPLYlink written 2.7 years ago by Manvendra Singh1.9k

Yeah, it was not clear from the OP. I first thought that "Sequence unavailable" was part of the FASTA header until they posted an actual example.

ADD REPLYlink written 2.7 years ago by Siva1.5k

OP asked this as a follow up to another question I answered yesterday. It was midnight so I postponed answering the follow up question, and woke up to a new question addressing this specific issue - I guess OP got a bit impatient and opened a new post without giving it the full context.

ADD REPLYlink written 2.7 years ago by Ram12k
1

Oh! Yes. I remember. Its the post where you and Pierre were answering. 

Dear seta , keep patience dear, its very important in science.

ADD REPLYlink written 2.7 years ago by Manvendra Singh1.9k
3

And if you take the code snippets, please try to understand them and don't just use them. Learning by doing, you know? ;)

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by David Langenberger7.9k
1

Yeah, I don't mind though - less work for me :)

ADD REPLYlink written 2.7 years ago by Ram12k
3
gravatar for Manvendra Singh
2.7 years ago by
Manvendra Singh1.9k
Berlin, Germany
Manvendra Singh1.9k wrote:

I have modified the perl script which should work for your cause

#!/usr/bin/perl
use strict;
use warnings;

$/="\n>";

while (<>) {

 s/>//g;

  my ($id, $seq) = split (/\n/, $_);

  print ">$_" if ((length $seq) > 10 && $seq !~ "Sequence unavailable");

}

 

hth

ADD COMMENTlink written 2.7 years ago by Manvendra Singh1.9k

thanks a lot dear friend.

ADD REPLYlink written 2.7 years ago by seta920
2
gravatar for David Langenberger
2.7 years ago by
Deutschland
David Langenberger7.9k wrote:

If there is no sequence below these headers, you could try something like this:

grep -v "Sequence unavailable" file.fasta > newFile.fasta

It will delete all headers containing "Sequence unavailable". But before you do that, double-check that there is no sequence under these headers!

ADD COMMENTlink written 2.7 years ago by David Langenberger7.9k

oh! I provided above script where sequences is not available at all.

does he has the  "Sequence unavailable" as the text in fasta record or its just empty which he has named as sequence unavailable?

ADD REPLYlink written 2.7 years ago by Manvendra Singh1.9k

the "Sequence unavailable" exist as a text, my mean is exactly like,">AT 1G01340|AT1G01340.

Sequence unavailable". thanks for your help me in advance.

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by seta920

Thanks, using this command just words "sequence unavailable" is disappeared, but its identifier is remained. would please help me to remove both of them from a fasta file?

ADD REPLYlink written 2.7 years ago by seta920

Can you provide an example of your fasta file? Just... 5 entries with at least one being the "sequence unavailable"? :)

ADD REPLYlink written 2.7 years ago by David Langenberger7.9k

Does it look like 1, 2?

1.

 

>AT 1G01340|AT1G01340. 
Sequence unavailable

2. 

 

>AT 1G01340|AT1G01340. Sequence unavailable
ADD REPLYlink written 2.7 years ago by David Langenberger7.9k

Sure. that's it:

>AT1G07745|AT1G07745.2
GCCTCCTTGGCTATATATTCCACGCCCATCGTTTGGATGTGCAAGAAGAGAACCACATAA
ATGAATGATCAGATGGATTTCTTGTTTCGGTTTTACATGTAGAATTATACTGCAGCAATG
AAAATTAAAATTTGTTGAACAAATGAGTGGTGTGTATATAAGGATCACATGAATATAGAA
TCACATTTTAAAACTTTGCGCGCACACC
>AT1G08845|AT1G08845.2
Sequence unavailable
>AT1G11870|AT1G11870.1
TTTTTAATGGTTAGCTCTTTGCAATTCTTGCTGTGATTCCTCTTGTACTTGAGGTTTTGT
TCATTATACTTTTCACAGGCGCTGAGTTTTTTGTGTTCTTGTCTCAGAACATACCATTTA
TAGTTTATATCATCATAACCGATGTTCTTTGATTCAATAATAACGAAAGGGAAGTATACA
TT
>AT1G10670|AT1G10670.3
GCTTCTCTCAGTACCTTCTATGACCAAAACTTTGTCTGTGTTTTAGAGCCTTTATTTACT
GTGGTTAAGATTACTCAAGCTATAAGATACTTGCAATTTCTTGAAAACTTCTGTTGTTCG
ATTCTCTTTCCCCTAACGTTTTCTTCAGATTCAATAAATAATCGTTACTTTTTTTTCCCC
TCT

>AT1G08010|AT1G08010.2
Seqence unavailable


 

ADD REPLYlink written 2.7 years ago by seta920
4

Try this:

grep -v "$(grep -B 1 "Sequence unavailable" file.fa)" file.fa
ADD REPLYlink written 2.7 years ago by David Langenberger7.9k
1

Many thanks David, it's great. It wil be great if you please share me how your perfect command is applied for removing short sequences( say less than 10 bp) with their identifiers, ?

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by seta920

Well, that's not possible with grep. Therefore, you might want to look into perl, or awk.

ADD REPLYlink written 2.7 years ago by David Langenberger7.9k

I'm new in this field and need more help. It's my lucky if you could please share me the useful script to this end (removing short sequences( say less than 10 bp) with their identifiers from fasta file)

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by seta920
2

Where would be the fun, if you don't try it yourself? Sorry, but this forum is not for getting free ready-to-use scripts. At least not from me.... don't be too angry with me. :)

ADD REPLYlink written 2.7 years ago by David Langenberger7.9k

Please accept my apologize for many questions. I'm learning all day and night, but as I never pass any course for even some basic concepts, I try to ask just my problems from your experts. 

ADD REPLYlink written 2.7 years ago by seta920

It's fine that you are a learner, seta. It's just that most of us do not have the time to tweak and debug scripts for anyone (including ourselves at times). You have to try and implement changes yourself. Just learn a bit of grep, sed and awk and you'll be all set. You can always read up on Python if the situation demands a more complicated approach.

ADD REPLYlink written 2.7 years ago by Ram12k
2

Not a great idea, asking for a tailored script. Not many here have the time to tweak and debug scripts for others.

ADD REPLYlink written 2.7 years ago by Ram12k

If you modify one line in Manvendra Singh's Perl script, you will get what you want. Change 1 to 10 or whatever sequence length threshold you want. I also added a condition to ignore sequences that are "Sequence unavailable".

print ">$_" if ((length $seq) > 10 && $seq !~ "Sequence unavailable");
ADD REPLYlink written 2.7 years ago by Siva1.5k

may be I write full script as an answer 

ADD REPLYlink written 2.7 years ago by Manvendra Singh1.9k
1
gravatar for Siva
2.7 years ago by
Siva1.5k
United States
Siva1.5k wrote:

Here's one way to do with grep in two steps

grep -B 1 'Sequence unavailable' INPUT_FILE | sort | uniq > TMP_FILE

grep -v -f TMP_FILE INPUT_FILE > OUTPUT_FILE

-INPUT_FILE is your multi-FASTA sequence file

-TMP_FILE is a unique list of FASTA headers you do not want and one 'Sequence unavailable' line

-OUT_FILE is the file without the 'Sequence unavailable' and their corresponding FASTA headers

 

 

ADD COMMENTlink written 2.7 years ago by Siva1.5k

thanks siva. I tried the first one and it's work well. the command that David suggested is exactlt my mean

ADD REPLYlink written 2.7 years ago by seta920

I did not see David's comment before posting which I think is better as it combines these two steps in to one without needing a temp file.

ADD REPLYlink written 2.7 years ago by Siva1.5k

It is actually the same command. The command from Siva is just easier to read and understand. :) 

ADD REPLYlink written 2.7 years ago by David Langenberger7.9k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1183 users visited in the last hour