In [1]:
from Bio import Entrez
In [ ]:
note changes in final py script
In [2]:
Entrez.email = 'hughcross@yahoo.com'
In [31]:
paraFile = open('/Users/hughcross/Analysis/DEVELOPMENT/obi_tools/rerun_wolves/obi2qiimeParameters.txt')
In [32]:
final_ranks = []
rank_abbrev = {}
for line in paraFile:
    if line.startswith('#'):
        continue
    else:
        line = line.strip('\n')
        parts = line.split('=')
        if parts[0] == 'email':
            email = parts[1]
        elif parts[0] == 'file':
            filename = parts[1]
        else:
            rank_abbrev[parts[0]]=parts[1]
            final_ranks.append(parts[1])

print(final_ranks)
print(rank_abbrev)
paraFile.close()
        

['class', 'order', 'suborder', 'family', 'genus', 'species']
{'f': 'family', 's': 'species', 'o': 'order', 'so': 'suborder', 'g': 'genus', 'c': 'class'}
In [33]:
email
Out[33]:
'hugh.cross@otago.ac.nz'
In [3]:
# final ranks 
ranks = ['class','order','suborder','family','genus','species']
abbrev = {'class':'c', 'order':'o','suborder':'so','family':'f','genus':'g','species':'s'}
In [96]:
ranks = ['phylum','class','order','family','genus','species']
abbrev = {'phylum':'p','class':'c','order':'o','family':'f','genus':'g','species':'s'}
In [4]:
obiFile = open('/Users/hughcross/Analysis/DEVELOPMENT/obi_tools/rerun_wolves/wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.sort.tab')
In [5]:
# parse file and get 23 (taxid) and 0 (id) and 2 (conf)
obidict = {}
#taxidDict = {}
taxidList = []
for line in obiFile:
    if line.startswith('id'):
        continue
    else:
        line = line.strip('\n')
        parts = line.split('\t')
        otu = parts[0]
        otuID = otu.replace('#','')
        taxid = parts[23]
        taxidList.append(taxid)
        conf = parts[2]
        obidict.setdefault(otuID, {})['taxid']=taxid
        obidict.setdefault(otuID, {})['confidence']=conf
print(len(obidict))
print(len(taxidList))
obiFile.close()

26
26
In [6]:
taxset = set(taxidList)
finTaxList = list(taxset)
print(len(taxset))
print(len(finTaxList))

12
12
In [ ]:
# use entrez to download taxonomy
In [7]:
# determine ranges for taxonomy searches
if len(finTaxList) < 10000:
    multiple_size = len(finTaxList)
    batch_size = multiple_size
    remainder = []
else:
    multiple_size = int(len(finTaxList)/10000)*10000
    batch_size = 10000
    remainder = finTaxList[multiple_size:]   # len(finTaxList) - multiple_size

print(multiple_size)
print(batch_size)
print(len(remainder))

12
12
0
In [13]:
taxa_dict = {}
for i in range(0,multiple_size,batch_size): # need to make this adjustable
    j = i+batch_size
    current = finTaxList[i:j]
    search1 = Entrez.read(Entrez.epost('Taxonomy', id=','.join(current)))
    webenv= search1['WebEnv']
    queryk = search1['QueryKey']
    fetch1 = Entrez.efetch(db='Taxonomy', retmode='xml', webenv=webenv, query_key=queryk)
    recs = Entrez.parse(fetch1)
    for record in recs:
        Species = record['ScientificName']
        Sp_Id = record['TaxId']
        low_rank = record['Rank']
        lineage_list = record['LineageEx']
        for line in lineage_list:
            name = line['ScientificName']
            rank = line['Rank']
            if rank in ranks:
                taxid = line['TaxId']
                taxa_dict.setdefault(Sp_Id, {})[rank]=name
                #taxid_dict[taxid]=name
        #taxa_dict.setdefault(Sp_Id, {})['species']=Species
        taxa_dict.setdefault(Sp_Id, {})[low_rank]=Species
    #taxid_dict[Sp_Id]=Species
    print(len(taxa_dict))
    
    # now get the rest
if len(remainder) > 0:
    search1 = Entrez.read(Entrez.epost('Taxonomy', id=','.join(remainder)))
    webenv= search1['WebEnv']
    queryk = search1['QueryKey']
    fetch1 = Entrez.efetch(db='Taxonomy', retmode='xml', webenv=webenv, query_key=queryk)
    recs = Entrez.parse(fetch1)
    for record in recs:
        Species = record['ScientificName']
        Sp_Id = record['TaxId']
        low_rank = record['Rank']
        print(low_rank)
        lineage_list = record['LineageEx']
        for line in lineage_list:
            name = line['ScientificName']
            rank = line['Rank']
            if rank in ranks:
                taxid = line['TaxId']
                taxa_dict.setdefault(Sp_Id, {})[rank]=name
                #taxid_dict[taxid]=name
        #taxa_dict.setdefault(Sp_Id, {})['species']=Species
        taxa_dict.setdefault(Sp_Id, {})[low_rank]=Species
    #taxid_dict[Sp_Id]=Species
    print(len(taxa_dict))
        
print('final count: ',len(taxa_dict))

12
final count:  12
In [15]:
taxa_dict[obidict['HELIUM_000100422_612GNAAXX:7:13:6954:130390/1_CONS_SUB_SUB']['taxid']]
Out[15]:
{'class': 'Mammalia',
 'family': 'Sciuridae',
 'order': 'Rodentia',
 'suborder': 'Sciuromorpha',
 'tribe': 'Marmotini'}
In [24]:
obidict['HELIUM_000100422_612GNAAXX:7:13:6954:130390/1_CONS_SUB_SUB']['confidence']
Out[24]:
'0.949494949495'
In [16]:
# final ranks 
final_ranks = ['class','order','suborder','family','genus','species']
rank_abbrev = {'class':'c', 'order':'o','suborder':'so','family':'f','genus':'g','species':'s'}
In [17]:
qii = {}
for k,v in obidict.items():
    
    lineage = taxa_dict[v['taxid']]
    taxlist = []
    for rank in final_ranks:
        abbrev = rank_abbrev[rank]
        if rank in lineage:
            string = abbrev+'__'+lineage[rank]
            taxlist.append(string)
            #taxstring = taxstring +abbrev+'__'+lineage[rank]+'; '
        else:
            string = abbrev+'__'+'unidentified'
            taxlist.append(string)

    # now go through and delete 'unidentified from the end
    ## first reverse order of list
    taxlist.reverse() # could also use revlist = reversed(taxlist)
    #print(taxlist)

    # now go through and as soon as there is no unidentified, stop
    for i in range(0,len(taxlist)):
        # go backwards
        if 'unidentified' not in taxlist[i]:
            newlist = taxlist[i:]
            break
    newlist.reverse()
    #print(newlist)
    taxstring = '; '.join(newlist)
    qii[k]=taxstring
print(len(qii))

26
In [18]:
qii
Out[18]:
{'HELIUM_000100422_612GNAAXX:7:103:5003:24510/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae',
 'HELIUM_000100422_612GNAAXX:7:10:6016:123650/1_CONS_SUB_SUB': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae; g__Canis',
 'HELIUM_000100422_612GNAAXX:7:111:11186:21420/1_CONS_SUB_SUB': 'c__Mammalia; o__Primates; so__Haplorrhini; f__Hominidae',
 'HELIUM_000100422_612GNAAXX:7:13:6954:130390/1_CONS_SUB_SUB': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae',
 'HELIUM_000100422_612GNAAXX:7:14:17532:55460/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae; g__Canis; s__Canis lupus',
 'HELIUM_000100422_612GNAAXX:7:20:19056:140680/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae',
 'HELIUM_000100422_612GNAAXX:7:22:6019:112400/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae',
 'HELIUM_000100422_612GNAAXX:7:22:8540:147080/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Capreolus; s__Capreolus capreolus',
 'HELIUM_000100422_612GNAAXX:7:29:15520:180350/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Capreolus; s__Capreolus capreolus',
 'HELIUM_000100422_612GNAAXX:7:30:17580:19880/1_CONS_SUB_SUB': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Bovidae; g__Rupicapra',
 'HELIUM_000100422_612GNAAXX:7:34:14050:193120/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Cervus; s__Cervus elaphus',
 'HELIUM_000100422_612GNAAXX:7:45:2528:202730/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae; g__Ammospermophilus; s__Ammospermophilus harrisii',
 'HELIUM_000100422_612GNAAXX:7:46:19314:177600/1_CONS_SUB_SUB': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Capreolus; s__Capreolus capreolus',
 'HELIUM_000100422_612GNAAXX:7:52:18333:138810/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Cervus; s__Cervus elaphus',
 'HELIUM_000100422_612GNAAXX:7:54:9724:163210/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae',
 'HELIUM_000100422_612GNAAXX:7:55:3707:41580/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae; g__Marmota',
 'HELIUM_000100422_612GNAAXX:7:57:18459:161450/1_CONS_SUB_SUB': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae; g__Marmota',
 'HELIUM_000100422_612GNAAXX:7:66:13553:61370/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Cervus; s__Cervus elaphus',
 'HELIUM_000100422_612GNAAXX:7:6:9274:149510/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia; f__Cervidae; g__Cervus; s__Cervus elaphus',
 'HELIUM_000100422_612GNAAXX:7:71:9435:201700/1_CONS_SUB_SUB': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae',
 'HELIUM_000100422_612GNAAXX:7:78:15278:32850/1_CONS_SUB_SUB': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae; g__Canis',
 'HELIUM_000100422_612GNAAXX:7:78:2653:179230/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae; g__Canis; s__Canis lupus',
 'HELIUM_000100422_612GNAAXX:7:81:16334:46450/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__unidentified; so__Ruminantia',
 'HELIUM_000100422_612GNAAXX:7:89:10281:59300/1_CONS_SUB_SUB': 'c__Mammalia; o__unidentified; so__Ruminantia',
 'HELIUM_000100422_612GNAAXX:7:95:12207:167170/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Carnivora; so__Caniformia; f__Canidae; g__Canis',
 'HELIUM_000100422_612GNAAXX:7:98:10056:187480/1_CONS_SUB_SUB_CMP': 'c__Mammalia; o__Rodentia; so__Sciuromorpha; f__Sciuridae; g__Marmota'}
In [19]:
qout = open('wolf_taxonomy_qii_format4.txt','w')
In [20]:
qout.write('Feature ID\tTaxon\tConfidence\n')
Out[20]:
28
In [25]:
for k,v in qii.items():
    conf = obidict[k]['confidence']
    qout.write(k+'\t'+v+'\t'+conf+'\n')
qout.close()