分析 DNA 序列中的串联重复基序答案

【问题标题】：Analyze tandem repeat motifs in DNA sequences分析 DNA 序列中的串联重复基序
【发布时间】：2018-12-19 03:53:42
【问题描述】：

Hy Py-guys :)。由于我是编码世界和 Python 的新手，因此我没有太多编码经验，因此将不胜感激。我正在处理 DNA 序列中的短串联重复序列，我想要一个代码，可以根据指定基因座的串联基序读取和计算重复核苷酸。

这是我需要的示例：

串联主题：

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

输入：

TTAGTTCAGGATAGTAGTTGTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTGCCAGGCTACCATGAAAGAp>

分析输入：

TTAGTTCAGGATAGTAGTTGTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTTAGAGp

输出：

AGAT AGAC (AGAT)2 GAT (AGAT)12

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

等位基因：16个

每个主题的总副本数 (1 + 1 + 2 + 12)

说明

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

所以在这种情况下，整个主题以AGAT 开头，以AGAT 的最后一个副本结束

在串联基序中指定的那些之间没有未知核苷酸（A/C/T/G），并且应该忽略该定义的基序之前和之后的所有内容

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

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

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

也可能有这种情况：

串联主题：

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

输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGGCTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA P>

分析输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGG的 CTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATC 强> ATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA P>

输出：

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

等位基因：28

(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]

还有更多……

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

【问题讨论】：

您需要使用您在tandem motif 中指定的格式（例如：您正在从其他地方读取它们）或者您可以更改格式？你用的是python3吗？
可以更改。它通常像我介绍的那样使用，但是对于编码，可以根据需要对其进行修改:)。是的，我有 Python3。
一个有 4 个字母的等位基因（例如：ACAT）的值为 1？
是的，主题中的每个“项目”都计为 1，除了小写字母、带有 N## 的那些以及添加 0,1 - 0,3 的那些，我猜这些将是最有问题的。通常他们有 4 个字母，但也有案例，有 3 个字母和 5 个字母。所以有时 3 个字母的图案算作 1，有时它们会加上 0,3 :/
..... 使其更简单：当图案有 4 个字母的“项目”计为 1 时，如果该图案中还有 3 个字母的项目则计为 0， 3. （例如：[TCTA],TCA,[TCTA]）但是当主题有 3 个字母的“项目”时，这些都算作 1（例如：[CCT],CTT,[TCT],CCT,[TCT]）。我希望这是有道理的：D

标签： python regex dna-sequence

【解决方案1】：

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

让我们看看正则表达式的工作原理：
下面是您的摘要在正则表达式模式中的样子：

总结

在主题中可以有这些类型，需要定义，并且每个基因座会有不同的基序组合：

方括号： 示例 [CTAT] – CTAT 的多个副本 - RegEx： (CTAT)+（一个或多个）或 (CTAT)*（零个或多个）

无括号： 示例 CTAT - 只有一份 CTAT - RegEx： (CTAT){1}

N#: 示例 N30 - 表示 30 个未指定的核苷酸 (A/C/G/T) - RegEx: .{30}

小写： 示例 ctat - 表示这些不包含在最终等位基因编号中 - 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))

输出：等位基因总数为：16

对于接下来的示例，我只显示正则表达式 tandem，其余代码相同

示例 2：

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

输出：等位基因总数为：32.2

示例 3：

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

输出：等位基因总数为：31.0

示例 4：

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

输出：等位基因总数为：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)+)"

开发一个完整的工作框架需要一点时间，具体取决于您想要实现的灵活性级别、输入类型、自动化...

【讨论】：

顺便说一句，问题写得很好，即使至少要解决 3 个不同的问题（通常每个问题都应该有 1 个问题）。见minimal reproducible example
哇，这看起来很棒，谢谢。所以我假设我可以创建像字典或元组之类的东西，在那里我可以用它们的正则表达式来定义所有基因座名称？那么我可以将这个字典应用于原始输入数据，然后分析它吗？
是的。如果我正在开发那部分，我可能会编写一个小函数，其中[ATAC] 被(ATAC)+) 取代，ata 被(?:ATA) 取代，N42 被(.{42}) 取代等等。请记住，您还可以使用 + 符号连接 RegEx 字符串。我建议您开发一些代码并在遇到困难时询问。如果您愿意，请通知我，因为我已经知道您的要求

【解决方案2】：

很抱歉，我没有时间完成所有案例，但希望这能给您一些帮助。

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

我们有一个位置 (pos) 变量来标记我们正在处理的序列中的索引。

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

现在我们已经初始化了设置，我们可以开始对串联主题字符串中的每个主题进行迭代。这是通过用逗号分隔字符串来完成的。

然后在 for 循环中，我们需要确定这是哪种主题情况（例如方括号重复、无括号、N# 等）。由于时间关系，我只为重复的方括号实现了这个，因为它很容易地演示了这个过程。

一旦测试用例通过，您需要处理需要完成的步骤。

例如，在这种带有方括号的情况下，主题是重复的，所以我将初始 count 变量初始化为 0，然后将 pos 跳转到 sequence 中主题的第一次出现处，如果 @ 987654329@ 是0 - 即，如果这是我们的第一个主题，我们需要将我们的位置跳到第一次出现的末尾。我还增加了count，因为我们找到了一个主题。

从这里开始，虽然sequence 中的下一个字符等于我们的主题字符串，但我们将位置增加主题的长度（因此它位于下一个主题的末尾）并增加 count .

最后，我们将格式化字符串 ((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}

以及我编造的一个非常基本的案例的测试：

>>> dna('[TCTA]', 'TGCAGCATTCTATCTATCTAGCTAAGCC')
{'alleles': 3, 'output': '(TCTA)3'}

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

【讨论】：