改进DNA比对去间隙的代码设计
Improving code design of DNA alignment degapping
这是一个关于更高效代码设计的问题：
假设三个对齐的 DNA 序列（seq1、seq2 和 seq3；它们都是字符串）代表两个基因（gene1 和 gene2）。这些基因的起始和终止位置相对于比对的 DNA 序列是已知的。
# Input
align = {"seq1":"ATGCATGC", # In seq1, gene1 and gene2 are of equal length
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,3], "gene2":[4,7]},
"seq3":{"gene1":[0,3], "gene2":[4,7]}}
我希望从比对中去除空隙（即破折号）并保持基因起始和终止位置的相对关联。
# Desired output
align = {"seq1":"ATGCATGC",
"seq2":"ATGC",
"seq3":"ACAC"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,1], "gene2":[2,3]},
"seq3":{"gene1":[0,1], "gene2":[2,3]}}
获得所需的输出并不像看起来那么简单。下面我为这个问题写了一些（行号）伪代码，但肯定有更优雅的设计。
1 measure length of any aligned gene # take any seq, since all seqs aligned
2 list_lengths = list of gene lengths # order is important
3 for seq in alignment
4 outseq = ""
5 for each num in range(0, length(seq)) # weird for-loop is intentional
6 if seq[num] == "-"
7 current_gene = gene whose start/stop positions include num
8 subtract 1 from length of current_gene
9 subtract 1 from lengths of all genes following current_gene in list_lengths
10 else
11 append seq[num] to outseq
12 append outseq to new variable
13 convert gene lengths into start/stop positions and append ordered to new variable
任何人都可以给我 suggestions/examples 更短、更直接的代码设计吗？
以下匹配两个测试用例的示例程序的输出：
align = {"seq1":"ATGCATGC",
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,3], "gene2":[4,7]},
"seq3":{"gene1":[0,3], "gene2":[4,7]}}
(START, STOP) = (0, 1)
for alignment, sequence in align.items():
new_sequence = ''
gap = 0
for position, codon in enumerate(sequence):
if '-' == codon:
for gene in annos[alignment].values():
if gene[START] > (position - gap):
gene[START] -= 1
if gene[STOP] >= (position - gap):
gene[STOP] -= 1
gap += 1
else:
new_sequence += codon
align[alignment] = new_sequence
运行的结果：
python3 -i test.py
>>> align
{'seq2': 'ATGC', 'seq1': 'ATGCATGC', 'seq3': 'ACAC'}
>>>
>>> annos
{'seq1': {'gene1': [0, 3], 'gene2': [4, 7]}, 'seq2': {'gene1': [0, 1], 'gene2': [2, 3]}, 'seq3': {'gene1': [0, 1], 'gene2': [2, 3]}}
>>>
我希望这仍然是优雅、直接、简短和 Pythonic 的！
我自己对上述问题的解决方案既不优雅也不Pythonic，但达到了所需的输出。非常欢迎任何改进建议！
import collections
import operator
# measure length of any aligned gene # take any seq, since all seqs aligned
align_len = len(align.itervalues().next())
# initialize output
align_out, annos_out = {}, {}
# loop through annos
for seqname, anno in annos.items():
# operate on ordered sequence lengths instead on ranges
ordseqlens = collections.OrderedDict()
# generate ordered sequence lengths
for k,v in sorted(anno.items(), key=operator.itemgetter(1)):
ordseqlens[k] = v[1]-v[0]+1
# start (and later append to) sequence output
align_out[seqname] = ""
# generate R-style for-loop
for pos in range(0, len(align[seqname])):
if align[seqname][pos] == "-":
try:
current_gene = next(key for key, a in anno.items() if a[0] <= pos <= a[1])
except StopIteration:
print("No annotation provided for position", pos, "in sequence", seqname)
# subtract 1 from lengths of current_gene
ordseqlens[current_gene] = ordseqlens[current_gene]-1
# append nucleotide unless a gap
else:
align_out[seqname] += align[seqname][pos]
# convert modified ordered sequence lengths back into start/stop positions
summ = 0
tmp_dict = {}
for k,v in ordseqlens.items():
tmp_dict[k] = [summ, v+summ-1]
summ = v+summ
# save start/stop positions to correct annos
annos_out[seqname] = tmp_dict
这段代码的输出是：
>>> align_out
{'seq3': 'ACAC',
'seq2': 'ATGC',
'seq1': 'ATGCATGC'}
>>> annos_out
{'seq3': {'gene1': [0, 1], 'gene2': [2, 3]},
'seq2': {'gene1': [0, 1], 'gene2': [2, 3]},
'seq1': {'gene1': [0, 3], 'gene2': [4, 7]}}
所以，我认为尝试将每个序列分解成基因然后删除破折号的方法会导致很多不必要的 book-keeping。相反，直接查看破折号然后根据它们的相对位置更新所有索引可能更容易。这是我编写的一个似乎运行正常的函数：
from copy import copy
def rewriteGenes(align, annos):
alignments = copy(align)
annotations = copy(annos)
for sequence, alignment in alignments.items():
while alignment.find('-') > -1:
index = alignment.find('-')
for gene, (start, end) in annotations[sequence].items():
if index < start:
annotations[sequence][gene][0] -= 1
if index <= end:
annotations[sequence][gene][1] -= 1
alignment = alignment[:index] + alignment[index+1:]
alignments[sequence] = alignment
return (alignments, annotations)
这会迭代每个比对中的破折号，并在删除时更新基因索引。
请注意，这会为以下测试用例生成索引为 [2,1] 的基因：
align = {"seq1":"ATGCATGC",
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq2":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq3":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]}}
这是必要的，因为您的索引设置方式不允许空基因。例如，索引 [2,2] 将是从索引 2 开始的长度为 1 的序列。
此答案处理您更新的 annos 词典，从评论到 cdlanes 答案。该答案使 annos 字典中 seq2 gene2 的索引不正确 [2,1]。如果序列包含该区域中的所有间隙，我建议的解决方案将从字典中删除 gene 条目。还要注意，如果一个基因在最后的 align 中只包含一个字母，那么 anno[geneX] 将具有相同的开始和停止索引 --> 请参阅您的 seq3 gene1评论 annos.
align = {"seq1":"ATGCATGC",
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,3], "gene2":[4,7]},
"seq3":{"gene1":[0,3], "gene2":[4,7]}}
annos3 = {"seq1":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq2":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq3":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]}}
import re
for name,anno in annos.items():
# indices of gaps removed usinig re
removed = [(m.start(0)) for m in re.finditer(r'-', align[name])]
# removes gaps from align dictionary
align[name] = re.sub(r'-', '', align[name])
build_dna = ''
for gene,inds in anno.items():
start_ind = len(build_dna)+1
#generator to sum the num '-' removed from gene
num_gaps = sum(1 for i in removed if i >= inds[0] and i <= inds[1])
# build the de-gapped string
build_dna+= align[name][inds[0]:inds[1]+1].replace("-", "")
end_ind = len(build_dna)
if num_gaps == len(align[name][inds[0]:inds[1]+1]): #gene is all gaps
del annos[name][gene] #remove the gene entry
continue
#update the values in the annos dictionary
annos[name][gene][0] = start_ind-1
annos[name][gene][1] = end_ind-1
结果：
In [3]: annos
Out[3]: {'seq1': {'gene1': [0, 3], 'gene2': [4, 7]},
'seq2': {'gene1': [0, 1], 'gene2': [2, 3]},
'seq3': {'gene1': [0, 1], 'gene2': [2, 3]}}
上述 3 基因 annos 的结果。只需替换 annos 变量：
In [5]: annos3
Out[5]: {'seq1': {'gene1': [0, 2], 'gene2': [3, 4], 'gene3': [5, 7]},
'seq2': {'gene1': [0, 1], 'gene3': [2, 3]},
'seq3': {'gene1': [0, 0], 'gene2': [1, 2], 'gene3': [3, 3]}}
这是一个关于更高效代码设计的问题：
假设三个对齐的 DNA 序列（seq1、seq2 和 seq3；它们都是字符串）代表两个基因（gene1 和 gene2）。这些基因的起始和终止位置相对于比对的 DNA 序列是已知的。
# Input
align = {"seq1":"ATGCATGC", # In seq1, gene1 and gene2 are of equal length
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,3], "gene2":[4,7]},
"seq3":{"gene1":[0,3], "gene2":[4,7]}}
我希望从比对中去除空隙（即破折号）并保持基因起始和终止位置的相对关联。
# Desired output
align = {"seq1":"ATGCATGC",
"seq2":"ATGC",
"seq3":"ACAC"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,1], "gene2":[2,3]},
"seq3":{"gene1":[0,1], "gene2":[2,3]}}
获得所需的输出并不像看起来那么简单。下面我为这个问题写了一些（行号）伪代码，但肯定有更优雅的设计。
1 measure length of any aligned gene # take any seq, since all seqs aligned
2 list_lengths = list of gene lengths # order is important
3 for seq in alignment
4 outseq = ""
5 for each num in range(0, length(seq)) # weird for-loop is intentional
6 if seq[num] == "-"
7 current_gene = gene whose start/stop positions include num
8 subtract 1 from length of current_gene
9 subtract 1 from lengths of all genes following current_gene in list_lengths
10 else
11 append seq[num] to outseq
12 append outseq to new variable
13 convert gene lengths into start/stop positions and append ordered to new variable
任何人都可以给我 suggestions/examples 更短、更直接的代码设计吗？
以下匹配两个测试用例的示例程序的输出：
align = {"seq1":"ATGCATGC",
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,3], "gene2":[4,7]},
"seq3":{"gene1":[0,3], "gene2":[4,7]}}
(START, STOP) = (0, 1)
for alignment, sequence in align.items():
new_sequence = ''
gap = 0
for position, codon in enumerate(sequence):
if '-' == codon:
for gene in annos[alignment].values():
if gene[START] > (position - gap):
gene[START] -= 1
if gene[STOP] >= (position - gap):
gene[STOP] -= 1
gap += 1
else:
new_sequence += codon
align[alignment] = new_sequence
运行的结果：
python3 -i test.py
>>> align
{'seq2': 'ATGC', 'seq1': 'ATGCATGC', 'seq3': 'ACAC'}
>>>
>>> annos
{'seq1': {'gene1': [0, 3], 'gene2': [4, 7]}, 'seq2': {'gene1': [0, 1], 'gene2': [2, 3]}, 'seq3': {'gene1': [0, 1], 'gene2': [2, 3]}}
>>>
我希望这仍然是优雅、直接、简短和 Pythonic 的！
我自己对上述问题的解决方案既不优雅也不Pythonic，但达到了所需的输出。非常欢迎任何改进建议！
import collections
import operator
# measure length of any aligned gene # take any seq, since all seqs aligned
align_len = len(align.itervalues().next())
# initialize output
align_out, annos_out = {}, {}
# loop through annos
for seqname, anno in annos.items():
# operate on ordered sequence lengths instead on ranges
ordseqlens = collections.OrderedDict()
# generate ordered sequence lengths
for k,v in sorted(anno.items(), key=operator.itemgetter(1)):
ordseqlens[k] = v[1]-v[0]+1
# start (and later append to) sequence output
align_out[seqname] = ""
# generate R-style for-loop
for pos in range(0, len(align[seqname])):
if align[seqname][pos] == "-":
try:
current_gene = next(key for key, a in anno.items() if a[0] <= pos <= a[1])
except StopIteration:
print("No annotation provided for position", pos, "in sequence", seqname)
# subtract 1 from lengths of current_gene
ordseqlens[current_gene] = ordseqlens[current_gene]-1
# append nucleotide unless a gap
else:
align_out[seqname] += align[seqname][pos]
# convert modified ordered sequence lengths back into start/stop positions
summ = 0
tmp_dict = {}
for k,v in ordseqlens.items():
tmp_dict[k] = [summ, v+summ-1]
summ = v+summ
# save start/stop positions to correct annos
annos_out[seqname] = tmp_dict
这段代码的输出是：
>>> align_out
{'seq3': 'ACAC',
'seq2': 'ATGC',
'seq1': 'ATGCATGC'}
>>> annos_out
{'seq3': {'gene1': [0, 1], 'gene2': [2, 3]},
'seq2': {'gene1': [0, 1], 'gene2': [2, 3]},
'seq1': {'gene1': [0, 3], 'gene2': [4, 7]}}
所以，我认为尝试将每个序列分解成基因然后删除破折号的方法会导致很多不必要的 book-keeping。相反，直接查看破折号然后根据它们的相对位置更新所有索引可能更容易。这是我编写的一个似乎运行正常的函数：
from copy import copy
def rewriteGenes(align, annos):
alignments = copy(align)
annotations = copy(annos)
for sequence, alignment in alignments.items():
while alignment.find('-') > -1:
index = alignment.find('-')
for gene, (start, end) in annotations[sequence].items():
if index < start:
annotations[sequence][gene][0] -= 1
if index <= end:
annotations[sequence][gene][1] -= 1
alignment = alignment[:index] + alignment[index+1:]
alignments[sequence] = alignment
return (alignments, annotations)
这会迭代每个比对中的破折号，并在删除时更新基因索引。
请注意，这会为以下测试用例生成索引为 [2,1] 的基因：
align = {"seq1":"ATGCATGC",
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq2":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq3":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]}}
这是必要的，因为您的索引设置方式不允许空基因。例如，索引 [2,2] 将是从索引 2 开始的长度为 1 的序列。
此答案处理您更新的 annos 词典，从评论到 cdlanes 答案。该答案使 annos 字典中 seq2 gene2 的索引不正确 [2,1]。如果序列包含该区域中的所有间隙，我建议的解决方案将从字典中删除 gene 条目。还要注意，如果一个基因在最后的 align 中只包含一个字母，那么 anno[geneX] 将具有相同的开始和停止索引 --> 请参阅您的 seq3 gene1评论 annos.
align = {"seq1":"ATGCATGC",
"seq2":"AT----GC",
"seq3":"A--CA--C"}
annos = {"seq1":{"gene1":[0,3], "gene2":[4,7]},
"seq2":{"gene1":[0,3], "gene2":[4,7]},
"seq3":{"gene1":[0,3], "gene2":[4,7]}}
annos3 = {"seq1":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq2":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]},
"seq3":{"gene1":[0,2], "gene2":[3,4], "gene3":[5,7]}}
import re
for name,anno in annos.items():
# indices of gaps removed usinig re
removed = [(m.start(0)) for m in re.finditer(r'-', align[name])]
# removes gaps from align dictionary
align[name] = re.sub(r'-', '', align[name])
build_dna = ''
for gene,inds in anno.items():
start_ind = len(build_dna)+1
#generator to sum the num '-' removed from gene
num_gaps = sum(1 for i in removed if i >= inds[0] and i <= inds[1])
# build the de-gapped string
build_dna+= align[name][inds[0]:inds[1]+1].replace("-", "")
end_ind = len(build_dna)
if num_gaps == len(align[name][inds[0]:inds[1]+1]): #gene is all gaps
del annos[name][gene] #remove the gene entry
continue
#update the values in the annos dictionary
annos[name][gene][0] = start_ind-1
annos[name][gene][1] = end_ind-1
结果：
In [3]: annos
Out[3]: {'seq1': {'gene1': [0, 3], 'gene2': [4, 7]},
'seq2': {'gene1': [0, 1], 'gene2': [2, 3]},
'seq3': {'gene1': [0, 1], 'gene2': [2, 3]}}
上述 3 基因 annos 的结果。只需替换 annos 变量：
In [5]: annos3
Out[5]: {'seq1': {'gene1': [0, 2], 'gene2': [3, 4], 'gene3': [5, 7]},
'seq2': {'gene1': [0, 1], 'gene3': [2, 3]},
'seq3': {'gene1': [0, 0], 'gene2': [1, 2], 'gene3': [3, 3]}}