Extract multiple fasta sequences from a fasta file based on sequenes

[entomology]

I've searched in the forum and google as well, but most of cases are extract sequences from fasta file based on id list. but now I only have several sequences without ids and i want to extract from a fasta file based on these sequences, does anyone have clues for it? Actually, that means I need the sequences identifier (or id) from the fasta file.

Thank you very much for your help.

[Brian Bushnell]

You can extract sequences that share kmers with your sequences with BBDuk:

bbduk.sh in=a.fa ref=b.fa out=c.fa mkf=1 mm=f k=31

This will print to C all the sequences in A that share 100% of their 31-mers with sequences in B.

You can also do something more precise with Dedupe, as it allows arbitrary set operations; so, let me know if the above method is insufficient.

[cmbetts]

There's almost certainly lots of ways to do this depending on what tools you're comfortable with.

In R/Bioconductor, you could read in the fasta using the bioconductor ShortRead package, and then use vcountPattern to identify the hits to your query sequences and write those as a new fasta file.

It's been a long time since I used it, but BioPython also has some nice iterators for going through fasta files, and would be better suited for a bigger fasta file. You would essentially write a little loop that iterates through the fasta and queries each record for your desired sequence. If it finds the sequence, write the record to a new file, otherwise move on to the next record.

I'm sure some folks might have some grep based methods for the commandline as well

[entomology]

3xs for the info.

Originally Posted by Brian Bushnell

You can extract sequences that share kmers with your sequences with BBDuk:

bbduk.sh in=a.fa ref=b.fa out=c.fa mkf=1 mm=f k=31

This will print to C all the sequences in A that share 100% of their 31-mers with sequences in B.

You can also do something more precise with Dedupe, as it allows arbitrary set operations; so, let me know if the above method is insufficient.

[entomology]

Actually, I do prefer command based under shell environment. grep seems a great way to do it, thanks for tips.

Originally Posted by cmbetts

There's almost certainly lots of ways to do this depending on what tools you're comfortable with.

In R/Bioconductor, you could read in the fasta using the bioconductor ShortRead package, and then use vcountPattern to identify the hits to your query sequences and write those as a new fasta file.

It's been a long time since I used it, but BioPython also has some nice iterators for going through fasta files, and would be better suited for a bigger fasta file. You would essentially write a little loop that iterates through the fasta and queries each record for your desired sequence. If it finds the sequence, write the record to a new file, otherwise move on to the next record.

I'm sure some folks might have some grep based methods for the commandline as well

[maubp]

Here are my Python scripts to do this, with Galaxy wrappers:

This filters the FASTA file (loads all the IDs, then goes through the FASTA file once):
https://github.com/peterjc/pico_gala...q_filter_by_id

This indexes the FASTA file, then goes through the IDs one by one:
https://github.com/peterjc/pico_gala...q_select_by_id

The difference is which order do you want? The order in the FASTA file (faster), or the order in the ID file (slower).

[entomology]

Seems these script also deal with a list of ids, then using these ids to fetch sequences in a fasta file. I wanna use sequences directly and get a subset of fasta from a big fasta file base on provided sequences. Thank you for your information as well.

Originally Posted by maubp

Here are my Python scripts to do this, with Galaxy wrappers:

This filters the FASTA file (loads all the IDs, then goes through the FASTA file once):
https://github.com/peterjc/pico_gala...q_filter_by_id

This indexes the FASTA file, then goes through the IDs one by one:
https://github.com/peterjc/pico_gala...q_select_by_id

The difference is which order do you want? The order in the FASTA file (faster), or the order in the ID file (slower).

[maubp]

What scripting/programming language(s) are you learning?

[entomology]

Actually, I do more wet lab most of the time. Sometimes, I'll also do simple small rna analysis with ready-to-use perl script and simple shell script. It's enough most of the time with my work, but sometime it's not easy for me.

Originally Posted by maubp

What scripting/programming language(s) are you learning?

[maubp]

Sadly my shell scripting skills are minimal, and my Perl worse, so I can't really help directly.

[entomology]

No worry, I'll try to use grep to deal with the problem .

Originally Posted by maubp

Sadly my shell scripting skills are minimal, and my Perl worse, so I can't really help directly.

[GenoMax]

Brian's solution should work. Did you try it?

While I like grep and its variants it may not always work for something as intricate as deciphering nucleotide patterns, specially if your sequences wrap around on multiple lines.

[entomology]

Yes, I've tried bbduk.sh.

bbduk.sh in=a.fa ref=b.fa out=c.fa mkf=1 mm=f k=31

But my situation is that b.fa is not fasta file, it contain one sequence per line. I just want the sequence in b from a.fa, than make a new fasta file (c.fa).

since my b.fa is not a fasta file, so bbduk.sh give some error:

Exception in thread "Thread-9" java.lang.RuntimeException: Error parsing read from text.

Originally Posted by GenoMax

Brian's solution should work. Did you try it?

While I like grep and its variants it may not always work for something as intricate as deciphering nucleotide patterns, specially if your sequences wrap around on multiple lines.

[GenoMax]

If your sequences are one on each line then use the following command to convert them to a fasta format file (change file names as needed)

$ awk -F "\n" 'BEGIN{counts=1}{print ">"counts"\n"""$0""; counts++}' your_file > new_file_as_fasta

Then use the file with BBDuk.

[entomology]

Thank you for the code. It can change my sequences to fasta file. And I try bbduk.fas again, but the result is not as expected. An example will be more easier to understand. there are two fasta

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


ref.fas(after I make fasta using your awk script)
>1
aaaaaa
>2
tttttt
>3
gggggg
>4
cccccc

I use "bbmap/bbduk.sh in=original.fas ref=ref.fas out=out.fas mkf=1 mm=f k=21"

out.fas is like this
>5-234
cggcgtcacgttggttgttga

actually, I want a fasta like this

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

Just like fetch the id from the original.fas

Originally Posted by GenoMax

If your sequences are one on each line then use the following command to convert them to a fasta format file (change file names as needed)

$ awk -F "\n" 'BEGIN{counts=1}{print ">"counts"\n"""$0""; counts++}' your_file > new_file_as_fasta

Then use the file with BBDuk.

[gsgs]

I usually write a small basic program for such problems.

post/send the file, I send the result ?

[Brian Bushnell]

Oh, I did not realize your sequences were so tiny; I assumed they were much longer. BBDuk is probably not appropriate for this situation, as it makes the implicit assumption that long kmers are relatively unique, which is not the case with 6-mers.

[GenoMax]

@entomology: Try the following

Use the original file of sequences (i.e. not the fasta format but just sequence, one on each line).

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

[SES]

Here is a simple script that uses an iterator to fetch records by sequence. This would be likely faster and less error-prone than grep:

#!/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;
open my $ini, '<', $infilei or die $!;
while (my ($id, $seq) = fasta_it(\*$ini)) {
    $hash{$seq} = $id;
}
close $ini;

open my $inj, '<', $infilej or die $!;
while (my ($id, $seq) = fasta_it(\*$inj)) {
    if (exists $hash{$seq}) {
	print join "\n", ">".$hash{$seq}, "$seq\n";
    }
}
close $inj;

sub fasta_it {
    my ($fh) = @_;
    
    local $/ = "\n>";
    return unless my $entry = $fh->getline;
    chomp $entry;

    my ($id, $seq) = split /\n/, $entry, 2;
    defined $id && $id =~ s/>//g;
    return ($id, $seq);
}

Here is the gist for easier download: https://gist.github.com/sestaton/889cba88b5279a58d997

The output:

perl fetch_by_seq.pl i.fas j.fas 
>1123-11234
aaaaaa
>232-23424
tttttt
>416-2
gggggg
>13424241234-23423
cccccc

Depending on the size of the original file you may want to think about using an SQLite database but this should work fine for most uses.

[entomology]

I've upload the two file, the expected output is like this:

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

Thanks!

Originally Posted by gsgs

I usually write a small basic program for such problems.

post/send the file, I send the result ?