所屬專題：[字符串算法](README.md)
 
## 問題
我們有一系列長度均為`$ n $`的DNA序列，它們的序列概況矩陣（profile matrix）是一個`$ 4\times n $`的矩陣`$ P $`. 其中，矩陣元素`$ P_{1,j} $`表示堿基A在各個序列的第`$ j $`位出現的次數；矩陣元素`$ P_{2,j} $`表示堿基C在各序列的第`$ j $`位出現的次數，以此類推。
一致性序列（consensus sequence）`$ c $`是一個長度同樣為`$ n $`的序列，由我們所掌握的序列集合中，各個位置出現次數最多的元素所組成。因此，序列`$ c $`的第`$ j $`位對應于序列概況矩陣第`$ j $`列中元素值最大的那個符號（堿基/氨基酸）。當然，某些位置上可能有不止一個最常見的符號，這使得存在多個可能的一致性序列。
DNA序列、序列概況矩陣（profile matrix）、一致性序列（consensus sequence）的關系如下所示：
*****
（1）多個等長的DNA序列：
:-: A T C C A G C T
G G G C A A C T
A T G G A T C T
A A G C A A C C
T T G G A A C T
A T G C C A T T
A T G G C A C T
（2）序列概況矩陣（profile matrix）
:-: **A**    5 1 0 0 5 5 0 0   
**C**    0 0 1 4 2 0 6 1   
**G**    1 1 6 3 0 1 0 0   
**T**    1 5 0 0 0 1 1 6   
（3）一致性序列（consensus sequence）
:-: A T G C A A C T
**輸入：** 包含多個等長的DNA序列的FASTA文件，序列不超過10條，長度均不超過1 kb.
**輸出：** 第一行是一致性序列，第二行至第五行是序列概況矩陣的內容。（如果存在多個一致性序列，則輸出符合條件的任一個均可）
**樣例數據：**
~~~
>Rosalind_1
ATCCAGCT
>Rosalind_2
GGGCAACT
>Rosalind_3
ATGGATCT
>Rosalind_4
AAGCAACC
>Rosalind_5
TTGGAACT
>Rosalind_6
ATGCCATT
>Rosalind_7
ATGGCACT
~~~
**樣例輸出：**
~~~
ATGCAACT
A: 5 1 0 0 5 5 0 0
C: 0 0 1 4 2 0 6 1
G: 1 1 6 3 0 1 0 0
T: 1 5 0 0 0 1 1 6
~~~
 
## 背景知識
該問題涉及序列概況矩陣與一致性序列的原理與應用。完成多序列比對后，通過構建序列概況矩陣，可以找到其中頻數較高的序列組成模式；而一致性序列則是構建了一個平均情況下出現頻率最高的序列組成，它代表了各個給定序列最相互接近的共同特征，這在序列的同源性分析中是很重要的概念。詳情請查閱ROSALIND網站上[關于該問題的背景說明](http://rosalind.info/problems/cons/)。
 
## 解答
```
def read_fa(fa):
"""讀入FASTA格式的文件對象fa，將其中的序列用字典結構存儲"""
seq = {}
for line in fa:
if line.startswith('>'):
seqid = line.replace('>', '').strip()
seq[seqid] = ''
else:
seq[seqid] += line.strip()
return seq
def profile_matrix(seq):
"""讀入字典結構的多個序列，構建其Profile matrix"""
sequences = list(seq.values())
seq_pos = list(zip(*sequences))
profile = {'A: ':[], 'C: ':[], 'G: ':[], 'T: ':[]}
for pos in seq_pos:
profile['A: '].append(pos.count('A'))
profile['C: '].append(pos.count('C'))
profile['G: '].append(pos.count('G'))
profile['T: '].append(pos.count('T'))
return profile
def consensus_seq(seq):
"""讀入字典結構的多個序列，構建一致性序列"""
sequences = list(seq.values())
seq_pos = list(zip(*sequences))
base_freq = map(lambda x:dict((base, x.count(base)) for base in "ACGT"), seq_pos)
conseq = ''.join([max(x, key=x.get) for x in base_freq])
return conseq
## --main--
with open("rosalind_cons.txt", 'r') as f1:
dnaseqs = read_fa(f1)
conseq = consensus_seq(dnaseqs)
profile_m = profile_matrix(dnaseqs)
with open("rosalind_cons_out.txt", 'w') as f2:
f2.write(conseq + '\n')
for base,profile in profile_m.items():
f2.write(base + ' '.join([str(x) for x in profile]) + '\n')
```
