Extract specific fasta sequences from multiple multi-fasta files located in a directory using Perl.

Question

[MB]

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
    QQQHDQVLPPPPLSLPPPLQSYSLHPPTPPTLPSDVSPYTSPADPLMFTGEPWDPPLSLD
    VLFPDPVAYSAYPSPEPTSSGSPTLTSTTEGRSGPAHGHTALPLALPLPDLLHEPTPGLT
    HRAALHGHEWTPLHLAARQGHTRVVQLLLTSDVHHRGAVNSVLAVLLRAGGDPSLPDSAG
    RTTRRGETAFHLAVQARQHAVVRVLLEHPASYYSVNAQDCWGRTPLHLACESNQQELVEM
    LILAGAQLDIRDFEDQTPLHSACLGG
    >OJJ94682_Homo_sapiens
    AQERKARLYAVFGGQGNDEHYFEELRTLWQTYRPLVEELIIPASSLLQRESMGGGSSRFY
    HNGLDALSWLEHPERTPASAYLIGAPISMPLIGLLQLAWYRVVGRVAGVTPAQLQTALAG
    TTGHSQGILPATVIAVAQSWTQFDALALDALRVLLAIGACSQDQFAVAQLPPAIVANAEA
    HGEGFPGPMLSVADFPVAQVRSYVADVNAHLPEDARIEIALINGPHNIVLAGPPLSLHGF
    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!

[finswimmer]

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

[MB]

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

Answer 1

[Dattatray Mongad]

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

listFile:- file containg headers

[finswimmer]

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

[Dattatray Mongad]

yes, faSomeRecords searches for headers only.

[MB]

Thanks, that helped!

Answer 2

[h.mon]

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.