[entomology]

Yes, that's does work. And it's more flexible and robust with different situations. Thank you for the great work. I've got a lot of fantastic script in bbmap!

Yep, it keeps the sequences in a.fasta that match sequences in b.fasta, and retains the names. You could alternatively run "filterbysequence.sh in=b.fasta ref=a.fasta out=c.fasta" to keep the names from b.fasta.

[Brian Bushnell]

Yep, it keeps the sequences in a.fasta that match sequences in b.fasta, and retains the names. You could alternatively run "filterbysequence.sh in=b.fasta ref=a.fasta out=c.fasta" to keep the names from b.fasta.

[entomology]

Thank you for your great work. I've tried the script with below two fasta files.

a.fasta

>1
aaaaaa
>2
tttttt
>3
gggggg
>4
cccccc



b.fasta

>1123-11234
aaaaaa
>wer
atgcca
>ad
ctaacg
>232-23424
tttttt
>323-342
cacaaa
>416-2
gggggg
>13424241234-23423
cccccc
>5-234
cggcgtcacgttggttgttga


running the script
filterbysequence.sh in=a.fasta ref=b.fas out=c.fasta ow=true

I've got c.fasta exactly same with a.fasta, is it supposed to replace identifier of the sequences in a.fasta from b.fasta?

So, I was inspired by this thread to add something into BBTools that could accomplish this. Thus, there's yet another method, "filterbysequence.sh". Usage:

filterbysequence.sh in=a.fasta ref=b.fasta out=c.fasta

c.fasta will then contain all sequences shared by a.fasta and b.fasta. It supports case-matching or case-insensitive operation, and reverse-complement-aware or forward-only operation. And it can either do an exclusion or inclusion filter. Also, it can optionally reduce very large sequences down to their 128-bit hash-codes for low-memory operation (so, for example, you could easily filter sequences against nt or RefSeq microbial quickly in a small amount of memory to see if they are already present before adding yet another copy of E.coli, which is something NCBI absolutely needs to do). And it's very, very fast.

[Brian Bushnell]

So, I was inspired by this thread to add something into BBTools that could accomplish this. Thus, there's yet another method, "filterbysequence.sh". Usage:

filterbysequence.sh in=a.fasta ref=b.fasta out=c.fasta

c.fasta will then contain all sequences shared by a.fasta and b.fasta. It supports case-matching or case-insensitive operation, and reverse-complement-aware or forward-only operation. And it can either do an exclusion or inclusion filter. Also, it can optionally reduce very large sequences down to their 128-bit hash-codes for low-memory operation (so, for example, you could easily filter sequences against nt or RefSeq microbial quickly in a small amount of memory to see if they are already present before adding yet another copy of E.coli, which is something NCBI absolutely needs to do). And it's very, very fast.

[entomology]

Another good method, thank you.

only those sequences with a "-" in the name ?

filtentr original 1 - > k1

[filters sequences (including name-lines) with - in the first entry]

names k1 5 5 > result
[only the raw sequence data, omitting the names]

using selfwritten simple programs filtentr.c and names.c

[entomology]

Yes, this time it works though the output is not as expected .

more out.fas
>1123-11234
--
gggggg
>13424241234-23423
>1123-11234
aaaaaa
ctaacg
>232-23424
>232-23424
tttttt
tttttt
>323-342
>416-2
gggggg
cacaaa
>416-2
>13424241234-23423
cccccc

Are you running bash shell? If you are not then try explicitly going into bash like this

$ /bin/bash
$ while read i ; do grep -B 1 $i original.fas ; done < sequence_file > out.fas

[entomology]

Admire you with the enthusiasm on programming, which makes an excellent programmer. Different running environment makes tons of unexpected problems that's why a robust code is really not easy to produce .

This is irrelevant now that the question is answered, but I can reproduce the "can't find getline..." error mentioned above that resulted from my first script.

The simple solution is to put "use IO::File;" at the top of the script. Apparently the behavior of Perl's I/O has changed (I'm using 5.22). Also, there are often OS-specific differences with scripting languages because the "perl" (or whatever) that the vendor distributes is not the same as what you compile yourself. So, it's possible there could be other tweaks required based on the system/versions, but it would be kind of a waste to go down that path since the problem is solved. Okay, now it's time to move on back to work...

[gsgs]

only those sequences with a "-" in the name ?

filtentr original 1 - > k1

[filters sequences (including name-lines) with - in the first entry]

names k1 5 5 > result
[only the raw sequence data, omitting the names]

using selfwritten simple programs filtentr.c and names.c

[SES]

This is irrelevant now that the question is answered, but I can reproduce the "can't find getline..." error mentioned above that resulted from my first script.

The simple solution is to put "use IO::File;" at the top of the script. Apparently the behavior of Perl's I/O has changed (I'm using 5.22). Also, there are often OS-specific differences with scripting languages because the "perl" (or whatever) that the vendor distributes is not the same as what you compile yourself. So, it's possible there could be other tweaks required based on the system/versions, but it would be kind of a waste to go down that path since the problem is solved. Okay, now it's time to move on back to work...

[GenoMax]

Forgive my poor programming skill, still I got some error message as below

-bash: syntax error near unexpected token `do'

Are you running bash shell? If you are not then try explicitly going into bash like this

$ /bin/bash
$ while read i ; do grep -B 1 $i original.fas ; done < sequence_file > out.fas

[SES]

Thank you very much for your help. you are right, i get the order wrong. now it works really good with the updated version without warning.

the command test gives the result as below:

perl -v

This is perl, v5.10.1 (*) built for x86_64-linux-thread-multi

perl -MIO::Handle -e 'print IO::Handle->VERSION'
1.28

thanks again for your patience and useful help!!!

Great, thanks for the info. I'm creating a CentOS 6 image now to test your errors above.

[entomology]

Thank you very much for your help. you are right, i get the order wrong. now it works really good with the updated version without warning.

the command test gives the result as below:

perl -v

This is perl, v5.10.1 (*) built for x86_64-linux-thread-multi

perl -MIO::Handle -e 'print IO::Handle->VERSION'
1.28

thanks again for your patience and useful help!!!

You must be giving the files in the wrong order because it works for me. It needs to be "perl script.pl original.fas other.fas." Though, the approach I posted is odd since there is something wrong with readline on your system and it will generate warnings. Here is a different approach that should work the same:

#!/usr/bin/env perl

use strict;
use warnings;
use File::Basename;

my $usage   = "perl ".basename($0)." seqsi.fas seqsj.fas > seqs_out.fas";
my $infilei = shift or die $usage;
my $infilej = shift or die $usage;

my $hash = index_seq($infilei);
open my $inj, '<', $infilej or die $!;
{
    local $/ = '>';

    while (my $line = <$inj>) {
	chomp $line;
	my ($seqid, @seqparts) = split /\n/, $line;
	my $seq = join '', @seqparts;
	next unless defined $seqid && defined $seq;
	if (exists $hash->{$seq}) {
	    print join "\n", ">".$hash->{$seq}, "$seq\n";
	}
    }
}
close $inj;

sub index_seq {
    my ($file) = @_;
    open my $in, '<', $file or die $!;

    my %hash;
    {
	local $/ = '>';

	while (my $line = <$in>) {
	    chomp $line;
	    my ($seqid, @seqparts) = split /\n/, $line;
	    my $seq = join '', @seqparts;
	    next unless defined $seqid && defined $seq;
	    $hash{$seq} = $seqid;
	}
    }
    close $in;

    return \%hash;
}

By the way, can you tell me the output of "perl -v" and "perl -MIO::Handle -e 'print IO::Handle->VERSION'". Thanks. I am also using CentOS, so I'm curious about that message.

[SES]

I'm using centos

2.6.32-573.7.1.el6.centos.plus.x86_64

when i change as you recommended, this time I got result. but still with some error and the result seems not as expected.

./fetch_fasta.test seq.test seq.test.original.fas
Value of <HANDLE> construct can be "0"; test with defined() at ./fetch_fasta.test line 31.
>1
aaaaaa
>2
tttttt
>3
gggggg

I expected the result like this:

>1123-11234
aaaaaa
>232-23424
tttttt
>416-2
gggggg
>13424241234-23423
cccccc

You must be giving the files in the wrong order because it works for me. It needs to be "perl script.pl original.fas other.fas." Though, the approach I posted is odd since there is something wrong with readline on your system and it will generate warnings. Here is a different approach that should work the same:

#!/usr/bin/env perl

use strict;
use warnings;
use File::Basename;

my $usage   = "perl ".basename($0)." seqsi.fas seqsj.fas > seqs_out.fas";
my $infilei = shift or die $usage;
my $infilej = shift or die $usage;

my $hash = index_seq($infilei);
open my $inj, '<', $infilej or die $!;
{
    local $/ = '>';

    while (my $line = <$inj>) {
	chomp $line;
	my ($seqid, @seqparts) = split /\n/, $line;
	my $seq = join '', @seqparts;
	next unless defined $seqid && defined $seq;
	if (exists $hash->{$seq}) {
	    print join "\n", ">".$hash->{$seq}, "$seq\n";
	}
    }
}
close $inj;

sub index_seq {
    my ($file) = @_;
    open my $in, '<', $file or die $!;

    my %hash;
    {
	local $/ = '>';

	while (my $line = <$in>) {
	    chomp $line;
	    my ($seqid, @seqparts) = split /\n/, $line;
	    my $seq = join '', @seqparts;
	    next unless defined $seqid && defined $seq;
	    $hash{$seq} = $seqid;
	}
    }
    close $in;

    return \%hash;
}

By the way, can you tell me the output of "perl -v" and "perl -MIO::Handle -e 'print IO::Handle->VERSION'". Thanks. I am also using CentOS, so I'm curious about that message.

[entomology]

I'm using centos

2.6.32-573.7.1.el6.centos.plus.x86_64

when i change as you recommended, this time I got result. but still with some error and the result seems not as expected.

./fetch_fasta.test seq.test seq.test.original.fas
Value of <HANDLE> construct can be "0"; test with defined() at ./fetch_fasta.test line 31.
>1
aaaaaa
>2
tttttt
>3
gggggg

I expected the result like this:

>1123-11234
aaaaaa
>232-23424
tttttt
>416-2
gggggg
>13424241234-23423
cccccc

Interesting, never seen that. What OS are you using? To fix the Perl issue, replace this line:

return unless my $entry = $fh->getline;

with

return unless my $entry = <$fh>;

should do the trick.

[SES]

I've got below message when I running the script, do i miss some module?

Can't locate object method "getline" via package "IO::Handle" at ./fetch_fasta.test line 30

Interesting, never seen that. What OS are you using? To fix the Perl issue, replace this line:

return unless my $entry = $fh->getline;

with

return unless my $entry = <$fh>;

should do the trick.