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
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
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!
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
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.
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");
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 fileTMP_FILE is a unique list of FASTA headers you do not want and one 'Sequence unavailable' lineOUT_FILE is the file without the 'Sequence unavailable' and their corresponding FASTA headersUse of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
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
Thanks, but it does not work at all.
Oh yes, I had thought that "Sequence unavailable" means that you have empty lines.
Anyways David and Siva have already given you answer, nevertheless, you can edit above script's last line
Instead of
(length $seq) > 1, write($seq) !~ "Sequence"It should work. Did you run it with your file? (
perl script.pl input.fasta)Would this not print everything? "Sequence unavailable" would also satisfy the condition.
Oh yes, I had thought that "Sequence unavailable" means that he has empty lines.
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.
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.
Oh! Yes. I remember. Its the post where you and Pierre were answering.
Dear seta, keep patience dear, its very important in science.
And if you take the code snippets, please try to understand them and don't just use them. Learning by doing, you know? ;)
Yeah, I don't mind though - less work for me :)