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

嗨，伙计们：）。由于我在编码领域和python领域都是新手，所以我在编码方面没有太多经验，因此任何帮助都是值得赞赏的。我正在研究DNA序列中的短串联重复序列，我想有一个基于特定位点串联基序的重复核苷酸读取和计数代码。
下面是我需要的一个例子：
串联图案：

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

输入：
附加可接受可接受可接受可接受可接受可接受可接受可接受可接受可接受
分析输入：
附加可接受可接受可接受可接受可接受可接受可接受可接受可接受可接受
输出：

AGAT AGAC (AGAT)2 GAT (AGAT)12

份数。（在输出中，GAT为大写，即使它不计算viz描述）
等位基因：16
每个图案的副本总数（1+1+2+12）
说明
每个位点的串联基序是不同的，所以我需要手动为每个位点指定它（总共约130个位点）。
因此在本例中，所有motif都以AGAT开头，以AGAT的最后一个副本结尾。
在串联基序中所指定的序列和在这个定义基序之前和之后的所有序列之间没有未知的核苷酸（A/C/T/G）
如你所见，当串联基序中有以小写字母（gat）书写的核苷酸时，它们不包括在最终等位基因值中。
那些在括号里的图案，可以重复多次
不在括号中的，序列中只有一个副本
也可能存在这种情况：
串联图案：

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

输入：
状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态附加
分析输入：
状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态状态附加
输出：

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

等位基因：28
（2+1+12+13）
说明
N30意味着在最后串联重复之前有30个未指定的核苷酸
摘要
母题中可以有这些类型，需要定义，每个轨迹都有不同的母题组合：
括号：示例[CTAT]：CTAT的多个副本
无括号：示例CTAT只有CTAT的一个副本
N'35；：示例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]

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

最佳答案

解决问题的一个好方法是使用regex。正则表达式是编程中解析strings的常用方法。
使用regex，您可以定义要在字符串中搜索的模式（几乎与您一样），这是问题的核心。
这意味着regex有自己的格式，与您的格式相似，但不完全相同。
您也可以编写一些代码来将格式转换为regex格式，但是您可能应该编写另一个问题，避免所有的dna内容。
让我们看看regex是如何工作的：
以下是您的摘要在regex模式中的外观：
摘要
在主题中可以有这些类型，需要定义，并且每个
轨迹将有不同的主题组合：
括号：示例[CTAT]–CTAT的多个副本-RegEx:(CTAT)+（一个或多个）或(CTAT)*（零个或多个）
无括号：示例CTAT–只有CTAT的一个副本-regex:(CTAT){1}
n：示例n30-表示30个未指定核苷酸（a/c/g/t）-正则表达式：.{30}
小写：示例ctat-表示这些不包括在最终等位基因数中-RegEx:(?:CTAT)
有了这些知识，我们可以将regex应用于我们的输入：
例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
对于下一个示例，我只显示regextandem，其余代码是相同的
例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)+)"

开发一个完整的工作框架需要花费一点时间，这取决于你想要达到什么样的灵活性，输入类型，自动化…