时间戳

首页 > 解决方案 > 分析 DNA 序列中的串联重复基序

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

问题描述

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

这是我需要的一个例子：

串联主题：

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

输入：

TTAGTTCAGGATAGTAGTTGTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCTAGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTTGCCAGGCTACCATGGAAAGACAACC

分析输入：

TTAGTTCAGGATAGTAGTTGTTGGAAGCGCAACTCTCTGAGAAACTTAGTTATTCTCTCATCTATTTAGCTACAGCAAACTTCATGTGACAAAAGCCACACCCATAACTTTTTTCCTCT AGATAGACAGATAGATGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGATAGAT ATAGATTCTCTTTCTCTGCATTCTCATCTATATTTCTGTCTTTCTCTTAATTATGGGTAACTCTTAGCCTTGCCAGGCTACCATGGAAAGACAACC

输出：

AGAT AGAC (AGAT)2 GAT (AGAT)12

等位基因：16

描述

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

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

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

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

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

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

也可能有这种情况：

串联主题：

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

输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGGCTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA

分析输入：

TTTGCATGATCTCTTCTTGATCATTTTCTTCCCCCTTTCCTAAAAAATTCTGGTCCTTTGAGGTAACTGCCATTACCATATGAGTTAGTCTGGGTTCTCCAGAGAAACAGAACCAATAGG CTATCTATCTAACTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTACTATCTCTATATTATCTATCTATCTATTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATC ATCTATCTATATCTTCTACCAAGTGATTTACTGTAATAAATTAGCTCATGCTATTATGGAGGATGAGTTCAAGATTTGTGGTCAGCAAGTTGCAGACTCA

输出：

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

等位基因：28

描述

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]

还有更多……</p>

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

标签： pythonregexdna-sequence

解决方案

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

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

概括

基序中可能有这些类型，需要定义，每个基因座都有不同的基序组合：

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

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

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

小写：示例 ctat - 表示这些不包括在最终等位基因编号中 - 正则表达式： (?: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

对于下一个示例，我只显示 regex 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)+)"

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

推荐阅读