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
    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!

fasta perl • 421 views
ADD COMMENTlink
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

ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
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

listFile:- file containg headers

ADD COMMENTlink
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

ADD REPLYlink
0
Entering edit mode

yes, faSomeRecords searches for headers only.

ADD REPLYlink
0
Entering edit mode

Thanks, that helped!

ADD REPLYlink
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.

ADD COMMENTlink

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.3.1