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
3 answers
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
•
written
by
Manvendra Singh
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!
Oh! I provided above script where sequences is not available at all.
Does he have the "Sequence unavailable" as the text in fasta record or its just empty which he has named as sequence unavailable?
•
written
by
Manvendra Singh
Can you provide an example of your fasta file? Just... 5 entries with at least one being the "sequence unavailable"? :)
•
written
by
David Langenberger
Does it look like 1, 2?
1.
>AT 1G01340|AT1G01340.
Sequence unavailable
2.
>AT 1G01340|AT1G01340. Sequence unavailable
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
Try this:
grep -v "$(grep -B 1 "Sequence unavailable" file.fa)" file.fa
Well, that's not possible with grep. Therefore, you might want to look into perl, or awk.
•
written
by
David Langenberger
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. :)
•
written
by
David Langenberger
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");
Maybe I'll write the full script as an answer
•
written
by
Manvendra Singh
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_FILEis your multi-FASTA sequence fileTMP_FILEis a unique list of FASTA headers you do not want and one 'Sequence unavailable' lineOUT_FILEis the file without the 'Sequence unavailable' and their corresponding FASTA headers
It is actually the same command. The command from Siva is just easier to read and understand. :)
•
written
by
David Langenberger
Log in to answer this question.
-
The VEP output for the variants of whole genome sequencing
written by seta
-
Getting Ensembl gene ID via biomart
written by seta
-
Getting allele frequency and zygosity information from several vcf file
written by seta
-
Why the ftp download of Ensemble is not available?
written by seta
-
Please help with comparing several gene ontology (GO) list
written by seta
-
How to deal with high percentage of non-specific matching during read mapping on Trinity assembly?
written by seta
-
What is the fastest way to extract a sequence ID from huge multiple FASTA file based on given seque…
written by Andrzej Zielezinski
-
How to remove some fasta sequences by header information from a large fasta file, any command and s…
written by seta
-
Please share the best tool for de novo motif discovery in large dataset
written by seta
-
Remove unwanted characters from FASTA file
written by seta
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.
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 :)
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"