如何在python中填充矩阵

Question

我试图编写一个需要输出为矩阵的代码，但由于是新手，我没有正确理解它。基本上我想为每列的A，C，G，T生成一个计数矩阵。我能够做到这一点，但很难为其他专栏做。

输入文件

>Rosalind_1 
ATCCAGCT 
>Rosalind_2 
GGGCAACT 
>Rosalind_3 
ATGGATCT 
>Rosalind_4 
AAGCAACC 
>Rosalind_5 
TTGGAACT 
>Rosalind_6 
ATGCCATT 
>Rosalind_7 
ATGGCACT 

到目前为止我的代码

fh_in = open("consensus_seq.txt", 'r') 

A_count = 0 
C_count = 0 
G_count = 0 
T_count = 0 

result = [] 
for line in fh_in: 
    line = line.strip() 
    if not line.startswith(">"): 
     for nuc in line[0]: 
      if nuc == "A": 
       A_count += 1 
      if nuc == "C": 
       C_count += 1 
      if nuc == "G": 
       G_count += 1 
      if nuc == "T": 
       T_count += 1 
result.append(A_count) 
result.append(C_count) 
result.append(G_count) 
result.append(T_count) 
print result 

输出

[5, 0, 1, 1] 

我想要的实际产量

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 

任何帮助/提示表示赞赏。

Answer 1

你可以使用numpy加载文本文件。由于格式是有点古怪，很难加载，但在这之后的总和变得微不足道：

import numpy as np 
data = np.loadtxt("raw.txt", comments=">", 
        converters = {0: lambda s: [x for x in s]}, dtype=str) 

print (data=="A").sum(axis=0) 
print (data=="T").sum(axis=0) 
print (data=="C").sum(axis=0) 
print (data=="G").sum(axis=0) 

输出：

[5 1 0 0 5 5 0 0] 
[1 5 0 0 0 1 1 6] 
[0 0 1 4 2 0 6 1] 
[1 1 6 3 0 1 0 0] 

的真正的好处，这是你所构建的numpy的阵列可以做其他事情。例如，假设我想知道的是我们在“Rosalinds”列中发现的A的平均次数，而不是总和：

print (data=="A").mean(axis=0) 
[ 0.71428571 0.14285714 0. 0. 0.71428571 0.71428571 0. 0.] 

Answer 2

import collections 
answer = [] 
with open('blah') as infile: 
    rows = [line.strip() for _,line in zip(infile, infile)] 
    cols = zip(*rows) 
    for col in cols: 
    d = collections.Counter(col) 
    answer.append([d[i] for i in "ATCG"]) 
    answer = [list(i) for i in zip(*answer)] 
for line in answer: 
    print(' '.join([str(i) for i in line])) 

输出：

5 1 0 1 0 
0 0 1 6 6 
5 0 2 0 1 
0 1 6 0 0 

Answer 3

首先使行的列表，剥出开始的行>。然后你可以zip这把它变成列表的列表。然后你可以列出每个字母的列数。

rows = [line.strip() for line in infile if not line.startswith('>')] 
columns = zip(*rows) 
for letter in 'ACGT': 
    print letter, [column.count(letter) for column in columns] 

但是，如果您的文件非常大，这可能会占用大量内存。另一种方法是逐行排列字母。

counts = {letter: [0] * 8 for letter in 'ACGT'} 
for line in infile: 
    if not line.startswith('>'): 
     for i, letter in enumerate(line.strip()): 
      counts[letter][i] += 1 
for letter, columns in counts.items(): 
    print letter, columns 

你也可以使用一个Counter，特别是如果你不知道提前多少列有如何将是：

from collections import Counter 
# ... 
counts = Counter() 
for line in infile: 
    if not line.startswith('>'): 
     counts.update(enumerate(line.strip())) 

columns = range(max(counts.keys())[0]) 
for letter in 'ACGT': 
    print letter, [counts[column, letter] for column in columns]