分析DNA序列中的串联重复模体

Question

帅哥们：）。我对Python编程没有太多的帮助，因为我在编写Python方面没有太多的经验。我正在研究DNA序列中的短串联重复序列，我想有一个代码，可以读取和计算基于特定位点串联基序的重复核苷酸。在

下面是我需要的一个例子：

---

串联图案：

AGAT,AGAC,[AGAT],gat,[AGAT]

输入：

TTAGTTCAGGATAGTAGTTGTTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTGCCAGGCTACCATGGAAAGACAACCTTTAT

分析输入：

TTAGTTCAGGATAGTAGTTGTTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTGCCAGGCTACCATGGAAAGACAACCTTTAT

输出：

^{pr2}$

• 份数。（在输出中，即使不计算viz描述，GAT也是大写的）

等位基因：16

• 各图案总印数（1+1+2+12）

---

说明

每个基因座的串联基序是不同的，所以我需要为每个基因座（总共约130个基因座）手动指定它。在

所以在本例中，整个motif以AGAT开头，以AGAT的最后一个副本结束

在串联基序中指定的核苷酸之间没有未知的核苷酸（A/C/T/G），在这个定义的基序之前和之后的一切都应该被忽略

如你所见，当串联基序中有以小写字母（gat）书写的核苷酸时，它们不包括在最终的等位基因值中

那些在括号里的图案，可以重复多次

那些不在括号中的-在序列中只有一个副本

---

也可能出现这种情况：

---

串联图案：

[CTAT],CTAA,[CTAT],N30,[TATC]

输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGGCTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA

分析输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGGCTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA

输出：

(CTAT)2 CTAA (CTAT)12 (TATC)13

等位基因

• （2+1+12+13）

  ---

说明

N30意味着在最后的串联重复之前有30个未指定的核苷酸

---

---

摘要

母题中可以有这些类型，需要加以界定，每个位点都有不同的母题组合：

括号：示例[CTAT]–CTAT的多个副本

无括号：示例CTAT–只有一份CTAT副本

N#：示例N30-表示30个未指定核苷酸（A/C/G/T）

小写：例如ctat-表示这些不包括在最终等位基因数中

---

真实母题的例子：

[CTTT],TT,CT,[CTTT]

[TCTA],[TCTG],[TCTA],ta,[TCTA],tca,[TCTA],tccata,[TCTA],TA,[TCTA]

[TAGA],[CAGA],N48,[TAGA],[CAGA]

[AAGA],[AAGG],[AAGA]

还有更多

---

提前谢谢大家。任何帮助和想法将不胜感激！：）

Answer 1

很抱歉，我没有时间完成所有的案件，但希望这会给你一些开始的东西。在

设计是一个线性自动机，其磁带是核苷酸序列。在

我们有一个position（pos）变量，它标记了我们正在使用的序列中的索引。在

还有两个正在运行的累积变量：一个output字符串和一个alleles的整数计数。在

现在我们已经初始化了设置，我们可以开始迭代串联motif字符串中的每个motif。这是通过在逗号上拆分字符串来完成的。在

然后在for循环中，我们需要确定这是哪个motif大小写（例如，方括号重复，无括号，N#等）。为了节省时间，我只对重复的方括号实现了这一点，因为它很容易演示过程。在

一旦测试用例通过了，您需要处理需要执行的步骤。在

例如，在这个方括号中，motif是重复的，所以我将一个初始的count变量初始化为0，然后跳pos到sequence中第一个出现的motif，如果pos是{}-也就是说，如果这是我们的第一个motif，我们需要将位置跳到第一个出现的末尾。我还增加了count，因为我们发现了一个母题。在

从这里开始，虽然sequence中的下一个字符等于我们的motif字符串，但是我们将位置增加为motif的长度（因此它位于下一个的末尾），同时还增加count。在

最后，我们将格式化字符串（(motif)#）附加到输出字符串，并将motif（等位基因）的数量添加到主alleles计数器中。在

然后我们将输出作为字典返回（如果需要，可以使用元组）。在

def dna(tandem_motif, sequence):
    pos = 0
    output = ''
    alleles = 0
    for motif in tandem_motif.split(','):
        if motif[0] == '[' and motif[-1] == ']':
            motif = motif.replace('[', ''). replace(']', '')
            count = 0
            if pos == 0:
                pos = sequence.index(motif) + len(motif)
                count += 1
            while sequence[pos:pos+len(motif)] == motif:
                pos += len(motif)
                count += 1
            output += '({}){}'.format(motif, count)
            alleles += count
        #elif ... :    <  where you add the criteria for the other motif test cases
    return {'alleles': alleles, 'output': output}

我做了一个很基本的测试：

^{pr2}$

这是正确的，因为：'TGCAGCAT|TCTATCTATCTA|GCTAAGCC'

Answer 2

解决问题的一个好方法是使用regex。正则表达式是编程中解析strings的常用方法。
使用regex，您可以定义要在字符串中搜索的模式（就像您所做的那样），这是问题的核心。
这意味着regex有自己的格式，与您的格式相似，但不完全相同。
您也可以编写一些代码来将格式转换为regex格式，但您可能应该编写另一个问题，避免所有的DNA内容。在

让我们看看regex是如何工作的：
以下是您的摘要在regex模式中的外观：

Summary

There can be these types in motifs, which need to be defined, and each locus would have different combination of motifs:

Brackets: example [CTAT] – multiple copies of CTAT - RegEx:(CTAT)+ (one or more) or (CTAT)* (zero or more)

No brackets: example CTAT – only one copy of CTAT - RegEx:(CTAT){1}

N#: example N30 - means 30 unspecified nucleotides (A/C/G/T) - RegEx:.{30}

Lower case: example ctat - means that these are not included in final allele number - RegEx:(?:CTAT)

有了这些知识，我们可以将正则表达式应用于我们的输入：
例1：

import re # import regex module

tandem = r"((AGAT){1}(AGAC){1}(AGAT)+(?:GAT){1}(AGAT)+)"

mystring = "TTAGTTCAGGATAGTAGTTGTTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTGCCAGGCTACCATGGAAAGACAACCTTTAT" #input string

analyzed_input = re.findall(tandem, mystring)[0]

print(analyzed_input) #see the match found

tot = 0
max_len = max((len(al) for al in analyzed_input[1:] if len(al) <= 4)) # longest allele, maximum 4
remaining_string = analyzed_input[0] #string to analyzed. will be cutted in for loop
for allele in analyzed_input[1:]: #for each allele
    section = re.findall(r"((" + re.escape(allele) + ")+)", remaining_string)[0][0] # section where the allele is repeated
    value = section.count(allele) if len(allele) == max_len else section.count(allele)*(len(allele)/10.0) # get the value of the alleles. /10.0 if allele is shorter than the longest allele found
    remaining_string = remaining_string[remaining_string.index(section)+len(section):] # cut away from remaining string the current section
    print("The value of allele {0} is {1}\n".format(allele, value))
    if len(allele) <= 4: #add the allele value if his length is less than 5
        tot += value

print("total allele number is: {0}".format(tot))

OUTPUT: total allele number is: 16

对于下一个示例，我只展示regex tandem，其余代码都是相同的

例2：

^{pr2}$

OUTPUT: total allele number is: 32.2

例3：

tandem3 = r"((TCTA)+(TCTG)+(TCTA)+(?:TA){1}(TCTA)+(?:TCA){1}(TCTA)+(?:TCCATA){1}(TCTA)+(TA)*(TCTA)*)"

OUTPUT: total allele number is: 31.0

例4：

tandem4 = r"((CTAT)+(CTAA){1}(CTAT)+(.{30})(TATC)+)"

OUTPUT: total allele number is: 28.0

另一个例子是：

[CTTT],TT,CT,[CTTT] r"((CTTT)+(TT){1}(CT){1}(CTTT)+)"

[TCTA],[TCTG],[TCTA],ta,[TCTA],tca,[TCTA],tccata,[TCTA],TA,[TCTA] r"((TCTA)+(TCTG)+(TCTA)+(?:TA){1}(TCTA)+(?:TCA){1}(TCTA)+(?:TCCATA){1}(TCTA)+(TA){1}(TCTA)+)"

[TAGA],[CAGA],N48,[TAGA],[CAGA] r"((TAGA)+(CAGA)+(.{48})(TAGA)+(CAGA)+)"

[AAGA],[AAGG],[AAGA] r"((AAGA)+(AAGG)+(AAGA)+)"

开发一个完整的工作框架需要一点时间，这取决于你想要达到的灵活性水平、输入类型、自动化程度。。。在