Python RNA.get_pr方法代码示例

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import get_pr [as 别名]
def getBPPM(sequence, structure = "", bppm_cutoff = 0.00001):
    """
        Requires ViennaRNAtools Python module
        Returns the base pair probability matrix using Vienna pf_fold, get_pr and free_pf_arrays functions.
        returns upper triangular matrix, whose entries exceed a threshold
    """
    bppm = {}
    
    
     #'--noPS', '-d 2', t, P
     
     
     
    if structure != "":
        RNA.cvar.fold_constrained = 1
    else:
        RNA.cvar.fold_constrained = 0
    #print "Before", structure
    RNA.pf_fold(sequence, structure)
    #print "After", structure
    seq_len = len(sequence)+1
    for i in xrange(1, seq_len):
        for j in xrange(1, seq_len):
            if i<j:
                bpp = RNA.get_pr(i,j)
                if bpp > bppm_cutoff:
                    bppm[str(i) + "_" + str(j)] = bpp
                else:
                    bppm[str(i) + "_" + str(j)] = 0
    RNA.free_pf_arrays()
    #print bppm
    #exit(1)
    return bppm

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import get_pr [as 别名]
def fold_probability(S, G=None):
    """Given a sequence S a secondary structure G (default mfe), we compute 
    the partition function of S given G as a constraint. The output
    is a triple (A,B,C) where A is the annotated partition folding,
    B is the energie of the ensemble A, and C a dictionary having as keys 
    a pair of positions and as value the probability of having the pair.
    """
    struct, energy = RNA.pf_fold(S, G) #Compute the partition function
    dict_probabilities = {}
    for left, right in ((x,y) for x in range(len(S)) for y in range(len(S))
                       if x < y):
        dict_probabilities[left,right] =RNA.get_pr(left + 1,right +1)
    return (struct, energy, dict_probabilities)

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import get_pr [as 别名]
def inserted_hairpin_stem(seq_five, seq_three, aptamer, shift):
    seq = terminal_hairpin(seq_five, aptamer, 0, shift)
    seq_ins = seq + seq_three
    overhang = len(seq_three)
    # print overhang
    RNA.pf_fold(seq_ins)
    bp_probability = [[RNA.get_pr(i, j) for i in range(0, len(seq_ins) + 1)] for j in range(0, len(seq_ins) + 1)]
    Iexists = 0
    IIexists = 0
    pI = 1
    pII = 1
    length = shift + overhang
    seq_length = len(seq_five)
    full_length = len(seq_ins)
    for k in range(seq_length - length, seq_length + shift - 1):
        if (bp_probability[k + 1][full_length - (k - (seq_length - length))] > 0.1) and (
            bp_probability[k + 1][full_length - (k - (seq_length - length))] < 0.2
        ):
            Iexists += 1
            pI *= bp_probability[k][full_length - (k - (seq_length - length))]
    for k in range(seq_length - overhang, seq_length + shift - 1):
        if (bp_probability[k + 1][full_length - (k - (seq_length - overhang))] > 0.1) and (
            bp_probability[k + 1][full_length - (k - (seq_length - overhang))] < 0.2
        ):
            IIexists += 1
            pII *= bp_probability[k][full_length - (k - (seq_length - overhang))]
    if length > seq_length:
        corr_shift = shift
    else:
        corr_shift = 0
    if (Iexists >= length - corr_shift) and (IIexists >= length - shift):
        # if abs(pI-pII) < 0.2:
        return seq_ins
    # else:
    #    return ""
    else:
        return ""

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import get_pr [as 别名]
def terminal_hairpin_stem(sequence, aptamer, overhang, shift):
    seq = terminal_hairpin(sequence, aptamer, overhang, shift)
    RNA.pf_fold(seq)
    bp_probability = [[RNA.get_pr(i, j) for i in range(0, len(seq) + 1)] for j in range(0, len(seq) + 1)]
    Iexists = 0
    IIexists = 0
    pI = 1
    pII = 1
    length = shift + overhang
    seq_length = len(sequence)
    full_length = len(seq)
    for k in range(seq_length - length, seq_length + shift - 1):
        if (bp_probability[k + 1][full_length - (k - (seq_length - length))] > 0.1) and (
            bp_probability[k + 1][full_length - (k - (seq_length - length))] < 0.99
        ):
            Iexists += 1
            pI *= bp_probability[k][full_length - (k - (seq_length - length))]
    for k in range(seq_length - overhang, seq_length + shift - 1):
        if (bp_probability[k + 1][full_length - (k - (seq_length - overhang))] > 0.1) and (
            bp_probability[k + 1][full_length - (k - (seq_length - overhang))] < 0.99
        ):
            IIexists += 1
            pII *= bp_probability[k][full_length - (k - (seq_length - overhang))]

    # print length
    if (Iexists >= length + shift - 2) and (IIexists >= length - 2):
        # if abs(pI-pII) < 0.2:
        print pI
        print pII
        print Iexists
        print IIexists
        return seq
    # else:
    #    return ""
    else:
        return ""

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import get_pr [as 别名]
def full_hairpin_stem(seq_five, seq_three, aptamer, shift, sec1, sec2):
    global _ACTIVE_TOLERANCE
    global _INACTIVE_TOLERANCE
    global _P_THRESHOLD
    global _MAX_ENERGY_DIFFERENCE
    global _ENERGY_THRESHOLD
    global _ENTROPY
    seq, aseq = full_hairpin(seq_five, seq_three, aptamer, shift)
    seq_ins = seq
    overhang = len(seq_three)
    # print overhang
    dump, a_free_E = RNA.pf_fold(aseq)
    bp_probability = [[RNA.get_pr(i, j) for i in range(1, len(aseq) + 1)] for j in range(1, len(aseq) + 1)]
    sec1_f, sec1_t = bond(sec1)
    sec2_f, sec2_t = bond(sec2)
    # print bond(sec2)
    Iexists = 0
    IIexists = 0
    # print len(seq_ins)
    pI = (len(seq_ins) * (len(seq_ins) - 1)) / 2 - len(sec1_f) * (len(seq_ins) - 1)
    # print pI
    pII = (len(seq_ins) * (len(seq_ins) - 1)) / 2 - len(sec2_f) * (len(seq_ins) - 1)
    # print pII
    length = shift + overhang
    seq_length = len(seq_five)
    full_length = len(seq_ins)
    for k in range(len(sec1_t)):
        # print [sec1_f[k], sec1_t[k]]
        # pI *= bp_probability[sec1_f[k]][sec1_t[k]]
        if bp_probability[sec1_f[k]][sec1_t[k]] > _P_THRESHOLD:
            # and (bp_probability[sec1_f[k]][sec1_t[k]-1] < 0.01):
            Iexists += 1
            pI *= bp_probability[sec1_f[k]][sec1_t[k]]
    # print "diff "+str(len(sec1_t)-Iexists)
    # print "1ex "+str(Iexists)+"\r"
    dump, b_free_E = RNA.pf_fold(seq_ins)
    bp_probability = [[RNA.get_pr(i, j) for i in range(1, len(seq_ins) + 1)] for j in range(1, len(seq_ins) + 1)]
    for k in range(len(sec2_t)):
        # print [sec2_f[k], sec2_t[k]]
        pII *= bp_probability[sec2_f[k]][sec2_t[k]]
        if bp_probability[sec2_f[k]][sec2_t[k]] > _P_THRESHOLD:
            IIexists += 1
            # pII *= bp_probability[sec2_f[k]][sec2_t[k]]
    # print pII, IIexists
    # print pII
    # print IIexists
    # if length > seq_length:
    #    corr_shift = shift
    # else:
    #    corr_shift = 0
    # print pI, pII
    S = sum([pair_entropy(seq, i) for i in range(len(seq))])
    if (
        (IIexists >= len(sec2_f) - _INACTIVE_TOLERANCE)
        and (Iexists >= len(sec1_f) - _ACTIVE_TOLERANCE)
        and (b_free_E < _ENERGY_THRESHOLD)
        and (S < _ENTROPY)
    ):
        # pI /= Iexists
        # pII /= IIexists
        if abs(a_free_E - b_free_E) < _MAX_ENERGY_DIFFERENCE:
            print "Free energy active: " + str(a_free_E), "Free energy inactive: " + str(
                b_free_E
            ), "Free energy difference: " + str(a_free_E - b_free_E)
            print "Entropy: " + str(S)
            print "Shift: " + str(shift)
            print "Active state base pairs: " + str(Iexists)
            print "Inactive state base pairs: " + str(IIexists)
            print "Sequence: " + str(seq_ins)
            print "good sequence"
            return seq_ins
        else:
            print "Free energy active: " + str(a_free_E), "Free energy inactive: " + str(
                b_free_E
            ), "Free energy difference: " + str(a_free_E - b_free_E)
            print "Shift: " + str(shift)
            print "Active state base pairs: " + str(Iexists)
            print "Inactive state base pairs: " + str(IIexists)
            print "Sequence: " + str(seq_ins)
            print "less-than-good sequence"
            return seq_ins
        # else:
        #    return ""
    else:
        return ""

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import get_pr [as 别名]
def pair_entropy(seq, base_idx):
    RNA.pf_fold(seq)
    prob = [[RNA.get_pr(i, j) for i in range(1, len(seq) + 1)] for j in range(1, len(seq) + 1)]
    s = sum([-elem * math.log(catch0(elem)) for elem in prob[base_idx]])
    return s