Subset FASTA file by species name I have a problem: I've managed to download a massive fasta file of 1500 sequences, but now I want to split them into separate fasta files based on the genus.
EDIT
The fasta file looks like this:
terminase_large.fasta
&gt;YP_009300697.1 terminase large subunit [Arthrobacter phage Mudcat]
MGLSNTATPLYYGQF...

&gt;YP_009208724.1 hypothetical protein ADP65_00072 [Achromobacter phage phiAxp-3]
MSNVLLKQELDEWL...

&gt;YP_009148449.1 large terminase subunit [Delftia phage RG-2014]
MSEPRKLVKKTLD...
...

I would like to end up with something like this:
mycobacterium_terminase.fasta
&gt;Mycobacterium phage JAMaL
MVRKKPPPELE...

&gt;Mycobacterium phage Bruin
MEVCGYTLDDI...

&gt;Mycobacterium phage Zaka
MSLDNHLPELA...

salmonella_terminase.fasta
&gt;Salmonella virus SETP7
MNVDITATEPQ...

&gt;Salmonella virus SE2
MEGGSDSLIAM...
and so on...

So that the label on the alignment would have the name of the phage and not the protein's ID.
I've managed to extract the genus names of my organism with this:
outfile = open('species.txt', 'w')
with open('terminase_large.fasta') as fd:
    for line in fd:
        if line.startswith('&gt;'):
            if '[' in line:
                name=line.split('[')[-1]
                name=name.split(' ', 1)[0]
                outfile.write(name[:] + &quot;\n&quot;)
outfile.close()

And got to extract only the unique names with this:
lines_seen = set()
outfile = open('species2.txt', &quot;w&quot;)
for line in open(&quot;species.txt&quot;, &quot;r&quot;):
    if line not in lines_seen:  # not a duplicate
        outfile.write(line)
        lines_seen.add(line)
outfile.close()

(Can I merge those two scripts together?)
Now, my genus names look like this:
Arthrobacter
Achromobacter
Delftia
....

I tried automating my script to get the Entrez data, but it gives me the 'Supplied id parameter is empty' message
My code looks like this:
from Bio import Entrez
Entrez.email = &quot;example@example.org&quot;

for line in open(&quot;species2.txt&quot;, &quot;r&quot;):
    searchterm = &quot;(terminase large subunit AND viruses[Organism]) AND&quot; +line+ &quot;AND refseq[Filter]&quot;
    searchResultHandle = Entrez.esearch(db=&quot;protein&quot;, term=searchterm, retmax=1000)
    searchResult = Entrez.read(searchResultHandle)
    ids = searchResult[&quot;IdList&quot;]

    handle = Entrez.efetch(db=&quot;protein&quot;, id=ids, rettype=&quot;fasta&quot;, retmode=&quot;text&quot;)
    record = handle.read()

    out_handle = open('terminase_large_'+str(line[:-1])+'.fasta', 'w')
    out_handle.write(record.rstrip('\n'))

Can someone help me with it?
Splitting into multiple files and changing the IDs can be easily done:

perl -pe 'if(/&gt;/){/\[(.*?)\]\s*$/; $_="&gt; $1\n"}' file.fa | 
    awk '(/^&gt;/){name=$2} {print &gt;&gt; name".fa"}'


That assumes all your FASTA headers have [foo bar baz] as the last element of a line. It will create a file called foo.fa (the bacterium's name) with all sequences saved there. 

Explanation


perl -pe : run the script given by -e on each line of the input file, and print the resulting line. 
if(/&gt;/) : if this line starts with a &gt;. \
/\[(.*?)\]\s*$/ : match an opening bracket (\[), then capture (that's what the parentheses do, they capture a pattern so we can refer to it as $1) everything until the first ] (.*?\])
$_="&gt; $1\n" : the $_ special variable in Perl is (in this case) the current line. So, $_=foo means "make the current line read foo. Since the -p prints each input line, changing the value of $_ means the changed value will be printed. So here, we are printing &gt;, whatever was in the square brackets ($1) and a newline character. 


The output of the perl command alone on your example input file is:

$ perl -pe 'if(/&gt;/){/\[(.*)\]\s*$/; $_="&gt; $1\n"}' file.fa
&gt; Arthrobacter phage Mudcat
MGLSNTATPLYYGQF...

&gt; Achromobacter phage phiAxp-3
MSNVLLKQELDEWL...

&gt; Delftia phage RG-2014
MSEPRKLVKKTLD...


So, we now pass it through awk which does:


(/^&gt;/){a=$2} : if this line starts with an &gt;, save the second field (the bacterial species) as the variable name. 
{print &gt;&gt; name".fa"} : print each line into a file whose name is the current value of the variable name with a .fa. extension. 




If you prefer python scripts to the one-liner approach, you can do the same thing with:

#!/usr/bin/env python
import re
outFile = None
with open("file.fa", "r") as inFile, open(species + ".fa", 'a') as outFile:
    for line in inFile:
        line = line.rstrip()
        if line.startswith('&gt;'):
            regex = re.compile('.*\[((.+?)\s+.*?)\].*')
            matches = regex.search(line)
            species = matches[2]
            outFile.write('&gt;%s\n' % matches[1])
        else:
            outFile.write("%s\n" % line)




As for your script, you've got the right idea, but have a small bug. You forgot to remove the \n from your input file, so it looks for Arthrobacter\n instead of Arthrobacter. The golden rule of debugging is "print all the things". If you add print("Searchterm: ",searchterm)  you will see:

Searchterm:  (terminase large subunit AND viruses[Organism]) ANDArthrobacter  
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDAchromobacter  
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDDelftia  
AND refseq[Filter]


So, you need to remove the newline characters and add a space like so (I also made it a bit more "pythonic" and conforming to the Python syntax guidelines):

#!/usr/bin/env python
from Bio import Entrez
Entrez.email = "example@example.org"

with open("species2.txt", "r") as in_handle:
    for line in in_handle:
        line = line.rstrip()
        searchterm = ("(terminase large subunit AND viruses[Organism]) " +
                      "AND %s AND refseq[Filter]" % line)
        print("Searchterm: ", searchterm)
        searchResultHandle = Entrez.esearch(db="protein",
                                            term=searchterm, retmax=1000)
        searchResult = Entrez.read(searchResultHandle)
        ids = searchResult["IdList"]

        handle = Entrez.efetch(db="protein", id=ids,
                               rettype="fasta", retmode="text")
        record = handle.read()

        with open('terminase_large_' + str(line[:-1]) + '.fasta', 'w') as out_handle:
            out_handle.write(record.rstrip('\n'))