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

Question

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

这是我需要的示例：

串联基序：

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

输入:

TTAGTTCAGGATAGTAGTTGTTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTGCCAGGCTACCATGGAAAGACAACCTTTAT

分析输入：

TTAGTTCAGGATAGTAGTTGTTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTGCCAGGCTACCATGGAAAGACAACCTTTAT

输出:

AGAT AGAC (AGAT)2 GAT (AGAT)12

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

等位基因: 16

每个主题的总拷贝数 (1 + 1 + 2 + 12)

描述

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

所以在这种情况下，整个主题以 AGAT 开始，以 AGAT

的最后一个副本结束

在串联基序中指定的核苷酸之间没有未知核苷酸 (A/C/T/G)，应忽略此定义基序前后的所有内容

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

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

不在括号内的——序列中只有一份

也可以有这种情况：

串联基序：

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

输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGGCTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA

分析输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGGCTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA

输出：

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

等位基因: 28

(2+1+12+13)

描述

N30 表示，在最终串联重复之前有 30 个未指定的核苷酸

总结

motif可以有这些类型，需要定义，每个locus会有不同的motif组合：

方括号：示例[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

很抱歉我没有时间完成所有案例，但希望这能给你一些开始。

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

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

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

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

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

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

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

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

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

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

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'

Answer 2

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

让我们看看正则表达式是如何工作的：
以下是您的摘要在正则表达式模式中的样子：

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

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

示例 2：

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

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)+)"

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

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

Analyze tandem repeat motifs in DNA sequences

python

regex

dna-sequence