Python SeqFeature.qualifiers['estimated_length']方法代码示例

# 需要导入模块: from Bio.SeqFeature import SeqFeature [as 别名]
# 或者: from Bio.SeqFeature.SeqFeature import qualifiers['estimated_length'] [as 别名]
def doConvert(embl_file, dep_file, contact, project, genome_project_id, organism_name, strain, locus_tag, taxon_id, dna_source, authors, comment, ac, clean=False):
    record = SeqIO.read(open(embl_file), "embl")

    # ----------------------------------------
    # HEADER
    # ----------------------------------------
    # remove accession
    if 'accession' in record.annotations.keys():
        del record.annotations['accession']
    record.annotations['accession'] = [ac]
    # ID line
    record.id = "XXX"
    record.name = "XXX"
    record.annotations['data_file_division'] = 'PRO'
    record.annotations['data_file_class'] = 'WGS'
    # PR line
    record.dbxrefs = ["Project:%s" % genome_project_id]
    # OS line
    record.annotations["organism"] = "%s %s" % (organism_name, strain)
    # DE line
    if project == 'metahit':
        record.description = "%s %s draft genome." % (organism_name, strain)
    else:
        record.description = "%s %s genome." % (organism_name, strain)
    # RN & RL lines
    if dna_source == 'GHP':
        dna_source = 'Rowett Institute of Nutrition and Health, University of Aberdeen -- http://www.rowett.ac.uk/divisions/ghp/'
        authors = 'Pajon A., Turner K., Parkhill J., Duncan S., Flint H.'
    elif dna_source == 'INRA':
        dna_source = 'INRA Clermont-Ferrand-Theix -- http://www.clermont.inra.fr/'
        authors = 'Pajon A., Turner K., Parkhill J., Bernalier A.'
    elif dna_source == 'HCIR':
        dna_source = 'Helmholtz Centre for Infection Research -- http://www.helmholtz-hzi.de/'
        authors = 'Pajon A., Turner K., Parkhill J., Timmis K., Oxley A., Wurdemann D.'
    elif dna_source == 'DSMZ':
        dna_source = 'German Collection of Microorganisms and Cell Cultures -- http://www.dsmz.de/'
        authors = 'Pajon A., Turner K., Parkhill J.'
    elif dna_source == 'NCTC':
        dna_source = 'Health Protection Agency\'s National Collection of Type Cultures -- http://www.hpacultures.org.uk/'
        authors = 'Pajon A., Turner K., Parkhill J.'
    elif dna_source == 'DPM':
        dna_source = 'Departments of Periodontology and Microbiology, King\'s College London -- http://www.kcl.ac.uk/'
        authors = 'Pajon A., Turner K., Parkhill J., Wade W., Vartoukian S.'
    else:
        dna_source = dna_source
        authors = authors
    ref_journal = Reference()
    ref_journal.journal = 'Unpublished.'
    if project == 'metahit':
        ref_journal.consrtm = "metaHIT consortium -- http://www.metahit.eu/"
    ref_journal.title = 'The genome sequence of %s %s' % (organism_name, strain)
    ref_journal.authors = authors
    ref_dep = Reference()
    ref_dep.authors = CONTACTS[contact]['author']
    today = date.today()
    ref_dep.journal = "Submitted (%s) to the EMBL/GenBank/DDBJ databases. Sanger Institute, Wellcome Trust Genome Campus, Hinxton, Cambridge CB10 1SA, United Kingdom." % today.strftime("%d-%b-%Y")
    ref_dep.title = 'Direct submission'
    record.annotations['references'] = [ref_journal, ref_dep]
    # CC line
    record.annotations['comment'] = ['Data release policy http://www.sanger.ac.uk/legal/#t_2',
                                     'DNA source: %s' % dna_source,
                                     '%s' % comment]
    
    # ----------------------------------------
    # GAP FEATURE (only with clean option)
    # ----------------------------------------
    # Add FT gap 
    seq = record.seq
    in_N = False
    gap_features = []
    if clean:
        # TODO - Cope with a sequence which ends with N
        if seq[-1] != "N":
            print "WARNING: sequence ends with N"
        for i in range(len(seq)):
            if seq[i] == 'N' and not in_N:
                start_N = i
                in_N = True
            if in_N and not seq[i+1] == 'N':
                end_N = i + 1
                length = end_N - start_N
                assert length > 0
                assert str(seq[start_N:end_N]) == "N"*length
                # do not create FT for 1bp gap
                if length > 1:
                    gap_feature = SeqFeature(FeatureLocation(start_N,end_N), strand=1, type="gap")
                    gap_feature.qualifiers['estimated_length'] = [length]
                    gap_features.append(gap_feature)
                in_N = False
    
    # ----------------------------------------
    # OTHER FEATURE (only with clean option)
    # ----------------------------------------
    new_features = []
    first_source = True
    has_source = False
    removed_cds = 0
    for i in range(len(record.features)):
        feature = record.features[i]

#.........这里部分代码省略.........