obiaddtaxids.py 16.2 KB
Newer Older
1
#!/usr/local/bin/python
2
'''
Aurélie Bonin committed
3 4
:py:mod:`obiaddtaxids`: adds *taxids* to sequence records using an ecopcr database
==================================================================================
5 6 7

.. codeauthor:: Celine Mercier <celine.mercier@metabarcoding.org>

Aurélie Bonin committed
8
The :py:mod:`obiaddtaxids` command annotates sequence records with a *taxid* based on 
Eric Coissac committed
9
a taxon scientific name stored in the sequence record header.
10

Aurélie Bonin committed
11 12
Taxonomic information linking a *taxid* to a taxon scientific name is stored in a 
database formatted as an ecoPCR database (see :doc:`obitaxonomy <obitaxonomy>`) or 
Eric Coissac committed
13
a NCBI taxdump (see NCBI ftp site).
14

Eric Coissac committed
15
The way to extract the taxon scientific name from the sequence record header can be
16
specified by two options:
17

18
    - By default, the sequence identifier is used. Underscore characters (``_``) are substituted
Eric Coissac committed
19 20 21
      by spaces before looking for the taxon scientific name into the taxonomic
      database.

Aurélie Bonin committed
22
    - If the input file is an ``OBITools`` extended :doc:`fasta <../fasta>` format, the ``-k`` option
Eric Coissac committed
23 24
      specifies the attribute containing the taxon scientific name.

Aurélie Bonin committed
25
    - If the input file is a :doc:`fasta <../fasta>` file imported from the UNITE or from the SILVA web sites,
Eric Coissac committed
26 27 28 29 30 31 32
      the ``-f`` option allows specifying this source and parsing correctly the associated 
      taxonomic information.
      
  
For each sequence record, :py:mod:`obiaddtaxids` tries to match the extracted taxon scientific name 
with those stored in the taxonomic database.

Aurélie Bonin committed
33
    - If a match is found, the sequence record is annotated with the corresponding *taxid*.
Eric Coissac committed
34

35
Otherwise,
36
    
Eric Coissac committed
37
    - If the ``-g`` option is set and the taxon name is composed of two words and only the 
38 39 40
      first one is found in the taxonomic database at the 'genus' rank, :py:mod:`obiaddtaxids` 
      considers that it found the genus associated with this sequence record and it stores this 
      sequence record in the file specified by the ``-g`` option.
41
    
Eric Coissac committed
42 43 44
    - If the ``-u`` option is set and no taxonomic information is retrieved from the 
      scientific taxon name, the sequence record is stored in the file specified by the 
      ``-u`` option.
45

Aurélie Bonin committed
46 47 48 49 50
    *Example*
    
    
    .. code-block:: bash
    
51
       > obiaddtaxids -k species_name -g genus_identified.fasta \\
Aurélie Bonin committed
52 53 54 55
                      -u unidentified.fasta -d my_ecopcr_database \\
                      my_sequences.fasta > identified.fasta

    Tries to match the value associated with the ``species_name`` key of each sequence record 
Aurélie Bonin committed
56
    from the ``my_sequences.fasta`` file with a taxon name from the ecoPCR database ``my_ecopcr_database``. 
Aurélie Bonin committed
57 58 59
        
            - If there is an exact match, the sequence record is stored in the ``identified.fasta`` file. 
        
Aurélie Bonin committed
60
            - If not and the ``species_name`` value is composed of two words, :py:mod:`obiaddtaxids` 
Aurélie Bonin committed
61 62 63 64 65 66 67
              considers the first word as a genus name and tries to find it into the taxonomic database. 
        
                - If a genus is found, the sequence record is stored in the ``genus_identified.fasta``
                  file. 
                  
                - Otherwise the sequence record is stored in the ``unidentified.fasta`` file.

Celine Mercier committed
68
'''
69

70

71
import re
72 73 74 75 76 77 78 79 80

from obitools.fasta import fastaIterator,formatFasta
from obitools.options import getOptionManager
from obitools.options.taxonomyfilter import addTaxonomyDBOptions
from obitools.options.taxonomyfilter import loadTaxonomyDatabase
from obitools.format.genericparser import genericEntryIteratorGenerator
from obitools import NucSequence


Celine Mercier committed
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
def addObiaddtaxidsOptions(optionManager):
    
    optionManager.add_option('-g','--genus_found',
                             action="store", dest="genus_found",
                             metavar="<FILENAME>",
                             type="string",
                             default=None,
                             help="(not with UNITE databases) file used to store sequences with the genus found.")

    optionManager.add_option('-u','--unidentified',
                             action="store", dest="unidentified",
                             metavar="<FILENAME>",
                             type="string",
                             default=None,
                             help="file used to store completely unidentified sequences.")

    optionManager.add_option('-s','--dirty',
                             action='store', dest="dirty",
                             metavar="<FILENAME>",
                             type="str",
                             default=None,
                             help="(not with UNITE databases) if chosen, ALL the words in the name used to identify the sequences will be searched"
                                  " when neither the exact name nor the genus have been found."
                                  " Only use if the sequences in your database are badly named with useless words or numbers"
                                  " in the name etc."
                                  " The sequences identified this way will be written in <FILENAME>.")
    
    optionManager.add_option('-f','--format',
                             action="store", dest="db_type",
                             metavar="<FORMAT>",
                             type="string",
                             default='raw',
113 114 115
                             help="type of the database with the taxa to be added. Possibilities : 'raw', 'UNITE_FULL', 'UNITE_GENERAL' or 'SILVA'."
                                  "The UNITE_FULL format is the one used for the 'Full UNITE+INSD dataset', and the UNITE_GENERAL format is the "
                                  "one used for the 'General FASTA release'."
116
                                  " Default : raw.")
Celine Mercier committed
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
    
    optionManager.add_option('-k','--key-name',
                             action="store", dest="tagname",
                             metavar="<KEYNAME>",
                             type="string",
                             default='',
                             help="name of the key attribute containing the taxon name in databases of 'raw' type. Default : the taxon name is the id "
                             "of the sequence. The taxon name MUST have '_' between the words of the name when it's the id, and "
                             "CAN be of this form when it's in a field.")

    optionManager.add_option('-a','--restricting_ancestor',
                             action="store", dest="res_anc",
                             type="str",
                             metavar="<ANCESTOR>",
                             default='',
                             help="can be a word or a taxid (number). Enables to restrict the search of taxids under a "
                                  "specified ancestor. If it's a word, it's the field containing the ancestor's taxid "
                                  "in each sequence's header (can be different for each sequence). If it's a number, "
                                  "it's the taxid of the ancestor (in which case it's the same for all the sequences)")



139 140 141 142 143 144 145 146
def numberInStr(s) :
    containsNumber = False
    for c in s :
        if c.isdigit() :
            containsNumber = True
    return containsNumber


147
def UNITEIterator_FULL(f):
148 149 150 151 152
    
    fastaEntryIterator = genericEntryIteratorGenerator(startEntry='>')
    for entry in fastaEntryIterator(f) :
        all = entry.split('\n')
        header = all[0]
153 154
        fields = header.split('|')        
        seq_id = fields[0][1:]
155
        seq = all[1]
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
        s = NucSequence(seq_id, seq)
        
        path = fields[1]

        species_name_loc = path.index('s__')
        species_name_loc+=3
        s['species_name'] = path[species_name_loc:]
        
        genus_name_loc = path.index('g__')
        genus_name_loc+=3
        s['genus_name'] = path[genus_name_loc:species_name_loc-4]
        
        path = re.sub('[a-z]__', '', path)
        s['path'] = path.replace(';', ',')

171 172 173
        yield s


174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
def UNITEIterator_GENERAL(f):
    
    fastaEntryIterator = genericEntryIteratorGenerator(startEntry='>')
    for entry in fastaEntryIterator(f) :
        all = entry.split('\n')
        header = all[0]
        fields = header.split('|')  
        
        seq_id = fields[0][1:]
        seq = all[1]
        s = NucSequence(seq_id, seq)
        
        s['species_name'] = seq_id.replace("_", " ")
        
        path = fields[4]
        path = re.sub('[a-z]__', '', path)
        path = path.replace(';', ',')
        s['path'] = path.replace(',,', ',')
        
        yield s


196
def SILVAIterator(f, tax):
197 198 199 200 201 202 203 204 205 206
    
    fastaEntryIterator = genericEntryIteratorGenerator(startEntry='>')
    for entry in fastaEntryIterator(f) :
        all = entry.split('\n')
        header = all[0]
        fields = header.split(' | ')
        id = fields[0][1:]
        seq = all[1]
        s = NucSequence(id, seq)
        
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238
        if (
            '(' in fields[1] 
            and len(fields[1].split('(')[1][:-1]) > 2 
            and ')' not in fields[1].split('(')[1][:-1] 
            and not numberInStr(fields[1].split('(')[1][:-1])
            ) :
            species_name = fields[1].split('(')[0][:-1]
            other_name = fields[1].split('(')[1][:-1]
            
            ancestor = None
            notAnAncestor = False
            
            if (len(other_name.split(' ')) == 1 and other_name[0].isupper()):
                
                try:
                    ancestor = tax.findTaxonByName(other_name)
                except KeyError :
                    notAnAncestor = True
            
            if (ancestor == None and notAnAncestor == False):
                s['common_name'] = other_name
                s['original_silva_name'] = fields[1]
                s['species_name'] = species_name
            
            elif (ancestor != None and notAnAncestor == False) :
                s['ancestor_name'] = other_name
                s['ancestor'] = ancestor[0]
                s['original_silva_name'] = fields[1]
                s['species_name'] = species_name
                
            elif notAnAncestor == True :
                s['species_name'] = fields[1]
Celine Mercier committed
239
                        
240
        else :
241 242
            s['species_name'] = fields[1]
        
243 244
        yield s

245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
    
def dirtyLookForSimilarNames(name, tax, ancestor):
    
    similar_name = ''
    taxid = None
    
    try :
        t = tax.findTaxonByName(name)
        taxid = t[0]
        similar_name = t[3]
    
    except KeyError :
        taxid = None
        
    if ancestor != None and not tax.isAncestor(ancestor, taxid) :
        taxid = None
    
    return similar_name, taxid
Celine Mercier committed
263
    
264 265 266

def getGenusTaxid(tax, species_name, ancestor):
    genus_sp = species_name.split(' ')
267 268 269 270
    genus_taxid = getTaxid(tax, genus_sp[0], ancestor)
    if tax.getRank(genus_taxid) != 'genus' :
        raise KeyError()
    return genus_taxid
271 272 273


def getTaxid(tax, name, ancestor):
274

275
    taxid = tax.findTaxonByName(name)[0][0]
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
    if ancestor != None and not tax.isAncestor(ancestor, taxid) :
        raise KeyError()
    return taxid


def get_species_name(s, options) :
    
    species_name = None
    if options.tagname == '' or options.tagname in s :
        if options.tagname == '' :
            species_name = s.id
        else :
            species_name = s[options.tagname]
                
        if "_" in species_name :
            species_name = species_name.replace('_', ' ')
          
        if len(species_name.split(' ')) == 2 and (species_name.split(' ')[1] == 'sp' or species_name.split(' ')[1] == 'sp.' or species_name.split(' ')[1] == 'unknown') :
            species_name = species_name.split(' ')[0]
          
        if options.tagname == '' :
            s['species_name'] = species_name
    
    return species_name


def getVaguelySimilarNames(species_name, tax, restricting_ancestor) :
    
    kindOfFound = False              
    uselessWords = ['sp', 'sp.', 'fungus', 'fungal', 'unknown', 'strain', 'associated', 'uncultured']
    for word in species_name.split(' ') :
        if word not in uselessWords :
            similar_name, taxid = dirtyLookForSimilarNames(word, tax, restricting_ancestor)
            if taxid != None :
                if len(similar_name) > len(s['species_name']) or kindOfFound == False :
                    s['species_name'] = similar_name
                    kindOfFound = True
    return kindOfFound


def openFiles(options) :
    
    if options.unidentified is not None:
        options.unidentified=open(options.unidentified,'w')
    
    if options.genus_found is not None:
        options.genus_found=open(options.genus_found,'w')
        
    if options.dirty is not None:
        options.dirty = open(options.dirty, 'w')


Celine Mercier committed
328
################################################################################################
329 330 331

if __name__=='__main__':

Celine Mercier committed
332
    optionParser = getOptionManager([addObiaddtaxidsOptions, addTaxonomyDBOptions], progdoc=__doc__)
333 334 335
    
    (options,entries) = optionParser()
    
336 337
    tax=loadTaxonomyDatabase(options)
    
338
    if options.db_type == 'raw' :
339
        entryIterator = fastaIterator
340
        entries = entryIterator(entries)
341 342 343 344
    elif options.db_type == 'UNITE_FULL' :
        entryIterator = UNITEIterator_FULL
    elif options.db_type == 'UNITE_GENERAL' :
        entryIterator = UNITEIterator_GENERAL
345
        entries = entryIterator(entries)
346
    elif options.db_type == 'SILVA' :
347
        entryIterator = SILVAIterator
348
        entries = entryIterator(entries, tax)
349
        options.tagname = 'species_name'
Celine Mercier committed
350

351 352
    openFiles(options)
    
353
    if (options.db_type == 'raw') or (options.db_type == 'SILVA') :
354 355 356 357 358 359 360 361 362
        
        if options.res_anc == '' :
            restricting_ancestor = None
        elif options.res_anc.isdigit() :
            restricting_ancestor = int(options.res_anc)
        
        for s in entries:
            
            if options.res_anc != '' and not options.res_anc.isdigit():
363
                restricting_ancestor = int(s[options.res_anc])
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
            
            species_name = get_species_name(s, options)
            
            if species_name != None :    
                try:
                    taxid = getTaxid(tax, species_name, restricting_ancestor)
                    s['taxid'] = taxid
                    print formatFasta(s)
                
                except KeyError:
                    
                    genusFound = False
                    if options.genus_found is not None and len(species_name.split(' ')) >= 2 :
                        try:
                            genusTaxid = getGenusTaxid(tax, species_name, restricting_ancestor)
379
                            s['genus_taxid'] = genusTaxid
380 381 382 383 384 385 386 387 388 389 390 391 392 393 394
                            print>>options.genus_found, formatFasta(s)
                            genusFound = True
                        except KeyError :
                            pass
                    
                    kindOfFound = False
                    if options.dirty is not None and not genusFound :
                        kindOfFound = getVaguelySimilarNames(species_name, tax, restricting_ancestor)
                        if kindOfFound == True :
                            print>>options.dirty, formatFasta(s)
                    
                    if options.unidentified is not None and not genusFound and not kindOfFound :
                        print>>options.unidentified,formatFasta(s)


395
    elif ((options.db_type =='UNITE_FULL') or (options.db_type =='UNITE_GENERAL')) :
396 397
        
        restricting_ancestor = tax.findTaxonByName('Fungi')[0][0]
398
                
399
        for s in entries :
400 401 402
            
            try :
                species_name = s['species_name']
403
                taxid = getTaxid(tax, species_name, restricting_ancestor)
404 405
                s['taxid'] = taxid
                s['rank'] = tax.getRank(taxid)
406
                print formatFasta(s)
407
                
408 409 410
            
            except KeyError:
                
411 412 413
                genusFound = False
                if options.genus_found is not None :
                    try:
414
                        genusTaxid = getGenusTaxid(tax, species_name, restricting_ancestor)
415 416 417
                        s['genus_taxid'] = genusTaxid
                        print>>options.genus_found, formatFasta(s)
                        genusFound = True
418
                    
419 420 421 422 423
                    except KeyError:
                        pass
                
                if options.unidentified is not None and not genusFound :
                    print>>options.unidentified,formatFasta(s)