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

Question

嗨，伙计们 :)。由于我是编码世界和 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 可以有这些类型，需要定义，每个位点会有不同的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

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

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