Biostar Beta. Not for public use.
Extract specific fasta sequences from multiple multi-fasta files located in a directory using Perl.
0
Entering edit mode
19 months ago
MB • 20

I have a fasta sequence in a file called "seq.fasta", which I want to search into multiple multi-fasta files located in a directory "/home/fasta_sequences" and print all the matching sequences to a file "output.fasta".

seq.fasta file looks like:

>KNH123.1_gallus_gallus
NDEAHPWKPLRPGDIRGPCPGLNTLASHGYLPRNGVATPVQIINAVQEGLNFDNQAAVFA
TYAAHLVDGNLITDLLSIGRKTRLTGPDPPPPASVGGLNEHGTFEGDASMTRGDAFFGNN
HDFNETLFEQLVD


/fasta_sequences/ directory looks like:

Homo_sapiens.AG_v1.pep.all.fa
Aspergillus_aculeatus_atcc_16872.Aspac1.pep.all.fa
.
.


* Homo_sapiens.AG_v1.pep.all.fa* file looks like:

    >OJJ94681_Homo_sapiens
MKITPIPVGVVGENWANIHDPVERRRAQNRIAQRGYRRRLRKGEERNTAQASAATLPSPP
VLFPDPVAYSAYPSPEPTSSGSPTLTSTTEGRSGPAHGHTALPLALPLPDLLHEPTPGLT
HRAALHGHEWTPLHLAARQGHTRVVQLLLTSDVHHRGAVNSVLAVLLRAGGDPSLPDSAG
RTTRRGETAFHLAVQARQHAVVRVLLEHPASYYSVNAQDCWGRTPLHLACESNQQELVEM
LILAGAQLDIRDFEDQTPLHSACLGG
>OJJ94682_Homo_sapiens
AQERKARLYAVFGGQGNDEHYFEELRTLWQTYRPLVEELIIPASSLLQRESMGGGSSRFY
HNGLDALSWLEHPERTPASAYLIGAPISMPLIGLLQLAWYRVVGRVAGVTPAQLQTALAG
TTGHSQGILPATVIAVAQSWTQFDALALDALRVLLAIGACSQDQFAVAQLPPAIVANAEA
NVWLRAVKAPASGVQHRIPFSQRRPEIRNRFLPITAAFHSSYLAEVVPEVLARVPGLTIS
GEQLRIPVYCTRTGLDLRTRGAANLVPELVAMITEQPLVWPKATQFPGATHILAFGPEGS
AGIGSLTNRLKEGTGVRTILASGGHHSHSARAEQLGDRS
.
.


I know how to search for a string in multiple files

 #!/usr/bin/perl
{
local @ARGV = glob "/home/fasta_sequences/*.fa";
while (<>) {
print "$ARGV:$.:\$_" if /TTHACVGSYYLSEDQKKLMNIFFTYRP/;
} continue {
close ARGV if eof;
}
}


But I am unable to find a way to do the same with a complete fasta sequence.

Any help would be appreciated!

fasta perl • 421 views
1
Entering edit mode

Hello MB,

you haven't just have the problem to find the sequence in another file. You also have to find out in which sequence it is contained if you have multiple sequences in one file.

So my question is if you realy want/need to use perl or if something like standard unix tools or e.g. seqkit is also suitable?

fin swimmer

0
Entering edit mode

Thanks for your answer. I will try to use grep command from seqkit.

2
Entering edit mode
19 months ago
National Centre for Cell Science, Pune

You can use faSomeRecords.

faSomeRecords - Extract multiple fa records

usage:

faSomeRecords in.fa listFile out.fa

options:

-exclude - output sequences not in the list file.


in.fa:- input fasta file

0
Entering edit mode

Hello doctor.dee005,

OP wants to look up given sequence in other files. But faSomeRecords is looking for sequence names, isn't it?

fin swimmer

0
Entering edit mode

yes, faSomeRecords searches for headers only.

0
Entering edit mode

Thanks, that helped!

2
Entering edit mode
17 months ago
h.mon 25k
Brazil

Do you have to use Perl? CDHIT has an utility which does (almost) what you want:

CD-HIT-2D compares 2 protein datasets (db1, db2). It identifies the sequences in db2 that are similar to db1 at a certain threshold. The input are two protein datasets (db1, db2) in fasta format and the output are two files: a fasta file of proteins in db2 that are not similar to db1 and a text file that lists similar sequences between db1 & db2.