rykahsay
commented
4 days ago Here is the script ---> /software/glygen/filter-refseq.py, also showing parts of it below. As you can see below, the downloaded files (excluding ["9606", "10090", "10116", "9823"]) cannot have tax_id in them.
def main():
species_obj = {}
in_file = "/data/projects/glygen/generated/misc/species_info.csv"
load_species_info(species_obj, in_file)
skip_list = ["9606", "10090", "10116", "9823"]
taxid2downloadfile = {
"559292":"/data/projects/glygen/downloads/ncbi/refseq/current/refseq_protein_all_fungi.gpff",
"7227":"/data/projects/glygen/downloads/ncbi/refseq/current/refseq_protein_all_insects.gpff",
"9031":"/data/projects/glygen/downloads/ncbi/refseq/current/refseq_protein_all_vertebrates.gpff",
"3702":"/data/projects/glygen/downloads/ncbi/refseq/current/refseq_protein_all_plants.gpff",
"44689":"/data/projects/glygen/downloads/ncbi/refseq/current/refseq_protein_all_amoebozoa.gpff"
}
for k in species_obj:
if species_obj[k]["is_reference"] != "yes" or k.isdigit() == False:
continue
tax_id = str(species_obj[k]["tax_id"])
long_name = species_obj[k]["long_name"]
if tax_id in skip_list:
continue
in_file = taxid2downloadfile[tax_id]
if os.path.isfile(in_file) == False:
continue
out_file = "/data/projects/glygen/downloads/ncbi/refseq/current/refseq_protein_all_%s.gpff" % (tax_id)
FW = open(out_file, "w")
record_count = 0
for record in SeqIO.parse(in_file, "genbank"):
organism = record.annotations['organism']
if organism.lower().find(long_name.lower()) != -1:
SeqIO.write(record, FW, "gb")
record_count += 1
FW.close()
cmd = "chmod 775 %s" % (out_file)
x = subprocess.getoutput(cmd)
print ("Created file: %s" % (out_file))
kmartinez834
Create script to separate records from refseq_proteinall*.gpff files. The script should accept an organism argument so it returns only records from the desired organism.
FYI @katewarner