Python RNA.export_bppm方法代码示例

# 需要导入模块: import RNA [as 别名]
# 或者: from RNA import export_bppm [as 别名]
    def add(self, seq):
        def bp_pr(i, j):
            if i > j:
                i, j = j, i
            return bppm[iindx[i + 1] - (j + 1)]

        def pos_bp_pr_iter():
            """Yield probability for base i being paired in each structure."""
            struct_0, struct_1 = target_structs
            for i in xrange(len(seq)):
                a = struct_0.basepairs[i]
                b = struct_1.basepairs[i]
                if a is not None:
                    if b is not None:
                        yield bp_pr(i, a), bp_pr(i, b)
                    else:
                        sum_ = sum((bp_pr(i, j) for j in xrange(size)))
                        yield bp_pr(i, a), 1 - sum_
                elif b is not None:
                    sum_ = sum((bp_pr(i, j) for j in xrange(size)))
                    yield 1 - sum_, bp_pr(i, b)
                else:
                    sum_ = 1 - sum((bp_pr(i, j) for j in xrange(size)))
                    yield sum_, sum_

        seq_str = str(seq)
        if seq_str in self.container:
            return
        size = len(seq)
        self.container[seq_str] = {}

        # with LOCK_VIENNA_RNA:
        energy_ensemble = vienna_rna.pf_fold_par(seq_str, None, None, 1, 0,
                                                 0)
        bppm = vienna_rna.doubleArray_frompointer(vienna_rna.export_bppm())
        i = vienna_rna.get_iindx(len(seq_str))
        iindx = vienna_rna.intArray_frompointer(i)
        energy_structs = (
            vienna_rna.energy_of_struct(seq_str, str(target_structs[0])),
            vienna_rna.energy_of_struct(seq_str, str(target_structs[1])))
        struct_mfe = " " * len(seq_str)
        energy_mfe = vienna_rna.fold_par(seq_str, struct_mfe, None, 0, 0)
        self.container[seq_str]["pos_bp_pr"] = list(pos_bp_pr_iter())
        self.container[seq_str]["energy_ensemble"] = energy_ensemble
        self.container[seq_str]["energy_mfe"] = energy_mfe
        self.container[seq_str]["struct_mfe"] = struct_mfe
        self.container[seq_str]["seq_pr"] = (
            math.exp((energy_ensemble - energy_structs[0]) / kT),
            math.exp((energy_ensemble - energy_structs[1]) / kT))

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