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.
The fasta file looks like this:
terminase_large.fasta >YP_009300697.1 terminase large subunit [Arthrobacter phage Mudcat] MGLSNTATPLYYGQF... >YP_009208724.1 hypothetical protein ADP65_00072 [Achromobacter phage phiAxp-3] MSNVLLKQELDEWL... >YP_009148449.1 large terminase subunit [Delftia phage RG-2014] MSEPRKLVKKTLD... ...
I would like to end up with something like this:
mycobacterium_terminase.fasta >Mycobacterium phage JAMaL MVRKKPPPELE... >Mycobacterium phage Bruin MEVCGYTLDDI... >Mycobacterium phage Zaka MSLDNHLPELA... salmonella_terminase.fasta >Salmonella virus SETP7 MNVDITATEPQ... >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('>'): if '[' in line: name=line.split('[')[-1] name=name.split(' ', 1) outfile.write(name[:] + "\n") outfile.close()
And got to extract only the unique names with this:
lines_seen = set() outfile = open('species2.txt', "w") for line in open("species.txt", "r"): 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 = "firstname.lastname@example.org" for line in open("species2.txt", "r"): searchterm = "(terminase large subunit AND viruses[Organism]) AND" +line+ "AND refseq[Filter]" 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() out_handle = open('terminase_large_'+str(line[:-1])+'.fasta', 'w') out_handle.write(record.rstrip('\n'))
Can someone help me with it?