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

# 需要导入模块: from Bio.SeqFeature import SeqFeature [as 别名]
# 或者: from Bio.SeqFeature.SeqFeature import qualifiers['note'] [as 别名]
def populateFeatures(tab_file):
    # read tab FT file
    """
FT  rRNA             1436971..1437083
FT                   /note="5s_rRNA"
FT                   /method="RNAmmer-1.2"

FT  tRNA             1614340..1614416
FT                   /note="tRNA-Arg(TCT) Cove Score 86.86"
FT                   /colour=4
FT                   /method="tRNAscan-SE"

FT   misc_feature    695000..700000
FT                   /colour=255 244 244
FT                   /algorithm="alien_hunter"
FT                   /note="threshold: 12.888"
FT                   /score=14.794
"""

    f_input = open (tab_file, 'r')
    ft_count = 0
    for line in f_input:
        values = line.strip().split()
        print values
        if values[0] == 'FT':
            if len(values) == 3:
                ft_type = values[1]
                print values[2]
                location = values[2].split('..')
                start = location[0]
                end = location[1]
                strand = 1
                if start[0:10] == 'complement':
                    start = start.split('(')[1]
                    end = end[:-1]
                    strand = -1
                if start[0] == '>':
                    start_position = AfterPosition(int(start[1:])-1)
                elif start[0] == '<':
                    start_position = BeforePosition(int(start[1:])-1)
                else:
                    start_position = ExactPosition(int(start)-1)
                if end[0] == '>':
                    end_position = AfterPosition(int(end[1:]))
                elif end[0] == '<':
                    end_position = BeforePosition(int(end[1:]))
                else:
                    end_position = ExactPosition(int(end))
                ft_count += 1
                feature = SeqFeature(FeatureLocation(start_position, end_position), strand=strand, type=ft_type)
                features.append(feature)
            if len(values) == 2:
                if values[1].startswith('/note'):
                    note = values[1].split('"')[1]
                    feature.qualifiers['note'] = [note]
                if values[1].startswith('/method'):
                    method = values[1].split('"')[1]
                    feature.qualifiers['method'] = [method]
                if values[1].startswith('/algorithm'):
                    method = values[1].split('"')[1]
                    feature.qualifiers['method'] = [method]
    f_input.close()
    print "INFO: %s features added from %s" % (ft_count, method)

# 需要导入模块: from Bio.SeqFeature import SeqFeature [as 别名]
# 或者: from Bio.SeqFeature.SeqFeature import qualifiers['note'] [as 别名]
def gene2features(r, gene, gene2position, gene2product, start, end, gcode, partialyes, verbose):
    """
    """
    contig, CDSs, gffstrand, function, frames = gene2position[gene]
    if gffstrand in ('1','+'):
        strand = +1
    else:
        strand = -1
        CDSs.reverse()
    '''#add stop codon if not partial seq
    if strand==1 and CDSs[-1][1]+3 <= len(r.seq):
            CDSs[-1][1] += 3
    elif strand==-1 and CDSs[0][0]-3 > 0:
        CDSs[0][0] -= 3'''
    cdsloc, mrnaloc = get_locations(CDSs, start, end, strand)
    #add gene
    geneid = gene #".".join(gene.split('.')[:-1])
    #get product
    product = "hypothetical protein"
    if geneid in gene2product:
        product = gene2product[geneid]    
    if gene.endswith('.t1'):
        sf = SeqFeature(FeatureLocation(BeforePosition(start-1),AfterPosition(end)), strand=strand, type='gene', id=geneid)
        sf.qualifiers={"locus_tag": geneid, "gene": geneid, "product": product}
        r.features.append(sf)
    #get mRNA sf
    sf = SeqFeature(mrnaloc, type='mRNA', id=gene)
    sf.qualifiers={"locus_tag": geneid, "gene": geneid, "product": product} #"protein_id": gene
    r.features.append(sf)
    #get CDS sf
    sf = SeqFeature(cdsloc, type='CDS', id=gene)
    #get translation
    seq = sf.extract(r.seq)
    aa = str(seq.translate(table=gcode))
    #solve non-triplets issue
    if len(seq) % 3:
        if strand==1:
            end   -= len(seq) % 3
        else:
            start += len(seq) % 3
    ##check for partial sequence - no M as first or no * as last aa
    partial = 0
    #both ends partial
    if aa[0]!="M" and aa[-1]!="*":
        partial = 1
        sf.location = FeatureLocation(BeforePosition(start-1),AfterPosition(end))
    #left end partial
    elif aa[0]!="M" and strand==1 or aa[-1]!="*" and strand==-1:
        partial = 1                
        sf.location = FeatureLocation(BeforePosition(start-1),end)
    #right end partial
    elif aa[-1]!="*" and strand==1 or aa[0]!="M" and strand==-1:
        partial = 1
        sf.location = FeatureLocation(start-1,AfterPosition(end))
    #strip stop codon
    aa = aa.strip("*")
    #replace internal stop codons by X
    if "*" in aa:
        if verbose:
            sys.stderr.write("[Warning] Stop codon(s) in: %s. Skipped!\n" % gene)
        return r
        #aa = aa.replace("*","X")
    sf.qualifiers = {'transl_table': gcode, "locus_tag": geneid, "gene": geneid, "product": product, "translation": aa} #"protein_id": gene,
    if function:
        sf.qualifiers['note'] = function
    #inform about partial entries
    if partial:
        #skip if not partial are allowed
        if not partialyes:
            return r
        if aa[0]!="M":
            sf.qualifiers['codon_start'] = 1
        sf.qualifiers['product']    += ", partial cds"
        if verbose:
            sys.stderr.write("[Warning] Partial sequence: %s\n" % (gene,))
            #sys.stderr.write("[Warning] Partial sequence: %s %s\n" % (gene,sf))
    #add to features
    r.features.append(sf)
    return r