I want to search sequences with length greater than 100. Will this program work if I change the sign >
ie. if($seq->length >= 100) {
thank you.raghul
Hi to all, I have a fasta file with lots of sequences with the description not like NCBI.But the file has keyword mentioning the length of the sequence. I want to extract sequences only of, say upto 100 bp in length.How can I do this with perl?
the format is as follows:
>gene1 group5 length=84
AaaaaatGTCTGGATGAGTTCATCCTGTAAAAaTTGCTGTCTGATACAAATACTTTGCTT
AGTCCAGTTAAATCTTCACTACTTTTGTGCACTGAAAGGCTAGCTTTtCTTCCAAAGCGG
TTTTCAATAATTCCTCTGACGCCTCCTTTTTtAGAGTATTTATTGTGTCTTCTATTTCCT
thanx raghul
Yes the length of the sequence is is not as mentioned. I did that for brevity(sorry!). I have many sequences of varying length. I am trying to write a program counting the strings. But I also felt it will be easy if I use the keyword length in FASTA description line to extract sequences within a range of values. thanx once again.raghul
NCBI fasta description is as follows,It has gi no.followed by reseq ID & organism etc which is NCBI way of describing the sequence but Mine is not so when u observe.
>gi|159476307|ref|XM_001696201.1| Chlamydomonas reinhardtii strain CC-503 cw92 mt+
CACAGTACCTTTCTGGTCAGCTGCACTGCATTGCTTTGTGACTAGTGAAGCTTCGACAGCTCACTGCGGA
CATTCCAAAATTGCTGTAACTCGACATTGATTTAACTACAGTATGCTGTTATATCCATAGCGCAAGAGAG
CTTGCGGCTTGCCTCCCCTCCATGCTCTTGTAGTCTGAGCCTATCCAGCTGCCTCGTCGCCGTTTGCAAA
GTTTTATTACTGAGACACAAGTAGCAGGGGCCGAGCAGGCAGCTGCCTGCGAGGCCGGTGAACCACGCGG
Most programming languages have libraries to parse Fasta (and many other sequence formats). Here's a quick example using the Bio::SeqIO library from BioPerl. It will create a new file of sequences 100 bp or less in length.
#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
my $seqin = Bio::SeqIO->new(-file => "myfile.fa", -format => "fasta");
my $seqout = Bio::SeqIO->new(-file => ">myfile_100.fa", -format => "fasta");
while(my $seq = $seqin->next_seq) {
if($seq->length <= 100) {
$seqout->write_seq($seq);
}
}
Similar solutions exist in BioPython, BioRuby, BioJava. There are pure "non-library" solutions, but I always recommend learning how to use a Bio* library. You can also search this site for "fasta parsing" to see more advice.
Since neilfws mentioned Biopython, here is a solution using it, although this does not satisfy the OP requirements:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
from Bio import SeqIO
fasta_file = sys.argv[1] # Input fasta file
min_length = int(sys.argv[2]) # Minimum length of sequence
result_file = sys.argv[3] # Output fasta file
fasta_sequences = SeqIO.parse(open(fasta_file),'fasta')
end = False
with open(result_file, "w") as f:
for seq in fasta_sequences:
if len(seq.seq.tostring()) > min_length:
SeqIO.write([seq], f, "fasta")
Should be used as follows (after saving to script.py and making the script executable):
script.py input.fasta 100 output.fasta
Cheers
Thanks for the answer. I use perl & currently learning to use Bioperl. I dont know python. So i am yet to try it. Thank u very much.raghul
perl -ne '/length=(\d+)/; print if $1>=100' file.fa
Seems to work. Neil's BioPerl method is far more reliable though ;)
Newer BLAST version(s) have blastdbcmd which is the successor to fastacmd. You can use range option in it. It selects the range of a sequence to extract.
blastdbcmd -range "< String >"
Range of sequence to extract (Format: start-stop)
In your case 1-100
Hi Thanx for the suggestion. Mine is local file with lot of sequences.blastdbcmd only works for blast formatted databases. My sequences do not follow NCBI fasta descriptions(u can find the differnce as explained in the question), so I have problems with formatting & then using it. If you have a solution for this,please let me know.
thanks once again raghul
">gene1 group5 length=84 " with above descriptions you can still format it in blast and use it give it a try
Ruby one-liner for the first option (with length in the title):
ruby -ne 'BEGIN{$/="\n>"}; length = $_.match(/length=(\d+)/)[1].to_i; print $_ if length <= 100' < in.fasta > out.fa
I have suggested one option in this thread
I have tested it by my own. It worls like gem
Log in to answer this question.
This example is somewhat confusing, since the length clearly is not 84.
This example is somewhat confusing, since the length clearly is not 84. Also, what is meant by "the description not like NCBI" ? Do you mean the format is not like the NCBI description of FASTA format ?
OK, so the sequence is valid fasta format; that's what I needed to know. The length description will work (as in Cass' answer, below), provided that it is correct, otherwise better to use the calculated length of the sequence.
thanx it worked for a newbie like me