Python RNA.pf_fold方法代码示例

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import pf_fold [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 pf_fold [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 pf_fold [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 pf_fold [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 pf_fold [as 别名]
def main():
    usage = """
    python scripts/dotplus.py sequence

    Create a file for displaying as a dotplot using the provided sequence.

    If the provided sequence is '-', then read the sequence from stdin.
    """
    num_args= 1
    parser = OptionParser(usage=usage)

    parser.add_option('-p', '--probability', dest='probability', default=0.01, help="The probability cutoff for displaying base pair points", type='float')
    #parser.add_option('-u', '--useless', dest='uselesss', default=False, action='store_true', help='Another useless option')

    (options, args) = parser.parse_args()

    if len(args) < num_args:
        parser.print_help()
        sys.exit(1)

    if args[0] == '-':
        seq = sys.stdin.read()
    else:
        seq = args[0].strip()

    # pf_fold returns the MFE as well as the partition function energy
    mfe, pfe = RNA.pf_fold(seq)
    prob_matrix = RNA.export_bppm()

    bp_to_seq = {}

    # get all of the suboptimal structures
    subopts = RNA.zukersubopt(seq)
    for i in range(0, subopts.size()):
        s = subopts.get(i)
        pt = fus.dotbracket_to_pairtable(s.structure)

        # go through each base pair and store the sequence
        # it came from so that it can be output later
        for j in range(1, pt[0]+1):
            bp = tuple(sorted([j, pt[j]]))

            if bp in bp_to_seq:
                continue

            bp_to_seq[bp] = (s.structure, s.energy)

        #print s.structure, s.energy

    bps = []
    structs = []

    counter = 0
    struct_dict = {}
    base_probs = col.defaultdict(float)

    print >>sys.stderr, "probability:", options.probability
    for i,j in it.combinations(range(1, len(seq)+1), 2):
        prob = RNA.doubleP_getitem(prob_matrix, 
                                   RNA.intP_getitem(RNA.cvar.iindx, i) - j)

        base_probs[i] += prob
        base_probs[j] += prob

        if prob > options.probability:
            struct, energy = bp_to_seq[(i,j)]
            pp = math.exp((pfe - (energy / 100.)) / .616310776)

            if struct not in struct_dict:
                struct_dict[struct] = (struct, pp, counter)
                index = counter
                counter += 1
            else:
                index = struct_dict[struct][2]

            print >>sys.stderr, "ix:", index, "struct:", struct, "pp:", pp
            bps += [{"i": i, "j": j, "p": "{:.3f}".format(math.sqrt(prob)), "ix": index}]


    structs = struct_dict.values()
    structs.sort(key=lambda x: -x[1])
    base_probs = dict([(i, "{:.3f}".format(j)) for (i,j) in base_probs.items()])

    structs = [{"struct": st[0], "sprob": st[1], "ix": st[2]} for st in structs]
    print json.dumps({"seq": seq, "structs": structs, "bps": bps, "baseProbs": base_probs.items()}, indent=2)

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import pf_fold [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 pf_fold [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