在给定结构中寻找RNA碱基配对

Question

我有一个RNA序列，看起来像这样。每个字符代表一个核苷酸（程序员可以忽略这一点，您可以将它们视为元素）：

         (((((((..((((.....(..)))).((((.........)))).....(((((..)....))))))))))))....

我将用这个术语来回答这个问题：

^{pr2}$

所以本质上，具有lhb的元素与具有rhb的元素相连，所有的点都是自由的地区。The它们的联系方式很复杂。很难用文字来表达，为了方便起见，我将把数字放在一些连接的元素下面：

^{3}$

我想这会让你知道他们是如何联系在一起的。我感兴趣的是找出连接的元素和自由区域的位置。在

（例如，元件1连接到元件72，元件8到9是自由的）。在

我 我选择了C来编写代码，但我没有接近逻辑的地方。在

而且，用C语言编程变得越来越困难，我觉得在python中使用正则表达式或perl可以很容易地做到这一点，但我对这些语言没有太多的经验。所以，如果有人能提供一个更简单的方法，那将是一个巨大的帮助。改进C代码的想法也很受欢迎。这是我的C代码：

#include <stdio.h>

int main() {

char dot[500];
int i = 0, j = 0;
int count = 0, count1 = 0, count2 = 0;
int lhb[100];
int rhb[100];
int dots[100];
int pair_1[100];
int pair_2[100];
int pair_3[100];
FILE * fp;

fp = fopen("structure.txt", "r");

while (fscanf(fp, "%c", & dot[i]) != EOF) {

  i++;
}

fclose(fp);

for (i = 0; dot[i] != '\0'; i++) {

  if (dot[i] == '(') {
    lhb[count] = dot[i];
    pair_1[count] = i;
    count++;
    } 
  else if (dot[i] == '.') {
    rhb[count1] = dot[i];
    pair_2[count1] = i;
    count1++;
    }
  else {
    dots[count2] = dot[i];
    pair_3[count2] = i;
    count2++;
    }
}


printf("Base-pair details :\n");

for (j = 0; j < count; j++)

  printf("%d--%d\n", pair_1[j] + 1, pair_3[count - j - 1] + 1);

printf("Loop details :\n");

// for(j=0;j<=count;j++)

// printf("--%d-",pair_2[j]+1);

return 0;

}

Answer 1

真正解决这个问题的最佳方法是维护堆栈。对于遇到的每一个(，您可以将索引值推送到堆栈中，对于每一个)，您需要使用插入最后一个(索引的数字弹出堆栈。这意味着(的索引和{}的索引形成一对。在

这可以通过这样做来实现

seq = '(((((((..((((.....(..)))).((((.........)))).....(((((..)....))))))))))))....'
xStack = []
for i, x in enumerate(seq):
    if x == '(':
        xStack.append(i)
    if x == ')':
        o = xStack.pop()

现在有了基本的步骤，除了维护括号的索引之外，还需要一些其他的东西。在pop操作之后，您需要存储匹配对，为此，我们引入另一个变量，当遇到.时，基本上什么都不做

^{pr2}$

现在我们得到了如下所示的结果对

^{3}$

我们需要找出所有的空闲空间在哪里，这可以很容易地通过以下操作来完成

spacesList = [i for i in range(len(seq)) if seq.startswith('.', i)]

结果是

[7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59, 72, 73, 74, 75]

现在可以很容易地编写一个函数，其中传递spacesList和{}，并获得每个可能对之间的可用空间数。有很多可能的优化，但这应该能够让你朝着正确的方向开始。在

def getSpacesCount(spacesList, resultingPairs):
    for pair in resultingPairs:
        a = pair[0]
        b = pair[1]
        spacesCount = 0
        for val in spacesList:
            if a < val < b:
                spacesCount+=1
        print a,b,spacesCount

seq = '(((((((..((((.....(..)))).((((.........)))).....(((((..)....))))))))))))....'
xStack = []
resultingPairs = []
for i, x in enumerate(seq):
    if x == '(':
        xStack.append(i)
    if x == ')':
        o = xStack.pop()
        tempPair = [o, i]
        resultingPairs.append(tempPair)
    if x == '.':
        pass

spacesList = [i for i in range(len(seq)) if seq.startswith('.', i)]

getSpacesCount(spacesList, resultingPairs)

有左括号的位置值，右括号的位置值，以及它们之间的自由空间数。在

>>> getSpacesCount(spacesList, resultingPairs)
18 21 2
12 22 7
11 23 7
10 24 7
29 39 9
28 40 9
27 41 9
26 42 9
52 55 2
51 60 6
50 61 6
49 62 6
48 63 6
9 64 28
6 65 30
5 66 30
4 67 30
3 68 30
2 69 30
1 70 30
0 71 30

编辑 似乎不知道如何编辑stackoverflow的注释，请将函数更新为

^{8}$

这将给你计数和位置。你喜欢哪一个你都可以留着。在

>>> getSpacesCount(spacesList, resultingPairs)
18 21 2 [19, 20]
12 22 7 [13, 14, 15, 16, 17, 19, 20]
11 23 7 [13, 14, 15, 16, 17, 19, 20]
10 24 7 [13, 14, 15, 16, 17, 19, 20]
29 39 9 [30, 31, 32, 33, 34, 35, 36, 37, 38]
28 40 9 [30, 31, 32, 33, 34, 35, 36, 37, 38]
27 41 9 [30, 31, 32, 33, 34, 35, 36, 37, 38]
26 42 9 [30, 31, 32, 33, 34, 35, 36, 37, 38]
52 55 2 [53, 54]
51 60 6 [53, 54, 56, 57, 58, 59]
50 61 6 [53, 54, 56, 57, 58, 59]
49 62 6 [53, 54, 56, 57, 58, 59]
48 63 6 [53, 54, 56, 57, 58, 59]
9 64 28 [13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
6 65 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
5 66 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
4 67 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
3 68 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
2 69 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
1 70 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]
0 71 30 [7, 8, 13, 14, 15, 16, 17, 19, 20, 25, 30, 31, 32, 33, 34, 35, 36, 37, 38, 43, 44, 45, 46, 47, 53, 54, 56, 57, 58, 59]

Answer 2

这是一个潜在的解决方案。Free包含一个索引每个点的列表，pairs包含一个元组的列表，这些元组包含每个对的索引。这是在假设您的数据是完美的（也就是说，左对和右对的数量相同，在左对之前没有右对）然而，可以修改它来检查边缘情况。根据你的问题，我从1开始索引，而不是0。在

data = '(((((((..((((.....(..)))).((((.........)))).....(((((..)....))))))))))))....'

left = []
group = []
prev = ''

free = []
pairs = []

for index, elem in enumerate(data, 1):
    if elem == '.' and prev == '.':
        group.append(index)
    elif elem == '.':
        group = [index]
    else:
        if len(group) >= 1:
            free.append(group)
            group = []
        if elem == '(':
            left.append(index)
        elif elem == ')':
            pairs.append([left.pop(), index])
    prev = elem
if len(group) > 0:
    free.append(group)
pairs.sort()

输出：

^{pr2}$

代码通过在每个元素上迭代一次来工作。如果遇到左括号，则将索引追加到列表中。当它遇到一个右paren时，它从最后一个看到的左paren的索引中弹出，创建一个对。至于自由点，跟踪前面的元素会告诉您如何处理每个点。如果已经看到一个，继续添加到当前列表，否则启动一个新列表。在

Answer 3

我不确定您到底想得到什么，但是这个perl生成的与Navidad20相同，它依赖于正则表达式。它使用Regexp::Common模块来获取平衡括号的位置，并使用一个简单的while循环来获取'free'元素的位置。从1开始计数，而不是从0开始计数。在

我不知道它在处理无效数据时会有什么反应。在

#!/usr/bin/perl
use strict;
use warnings;
use Regexp::Common qw /balanced/;

my $seq = '(((((((..((((.....(..)))).((((.........)))).....(((((..)....))))))))))))....';

while ($seq =~ /(?=($RE{balanced}{-parens=>'()'}))/g) {
    print 1 + $-[1], ' ',  $+[1], " $1\n";
}

my @free;
push @free, [ 1 + $-[0] .. $+[0] ] while $seq =~ /\.+/g;

use Data::Dumper; print Dumper \@free;

输出为：

^{pr2}$