读取 ATCG DNA 序列，计算第三位 ATCG 的数量

Posted 2023-02-25

【中文标题】读取 ATCG DNA 序列，计算第三位 ATCG 的数量

【英文标题】：Reading ATCG DNA sequences, and calculates the numbers of ATCG in third place

【发布时间】：2016-05-21 14:32:26

【问题描述】：

我想读取 ATCG DNA 序列，并计算第三位的 ATCG 数量。

例如1：

DNA = AAATTTCCCGGG

第三位ATCG是这样的：AA'A'TT'T'CC'C'GG'G'

所以在这个序列中 A=1 T=1 C=1 G=1。

例如2：

DNA = ATGGTATTTAAA

AT"G"GT"A"TT"T"AA"A"

我想计算 3、6、9、12 位 ATCG 号码。所以在 DNA 中 A=2 T=1 C=0 G=1

我的txt文件是这样的：

>seq1
ATGGTATTTAAA
ATCGTTTTTAAA
>seq2
ATGGTATTTAAA
ATCGTTTTTAAA
ATCGTTTTTAAA
>seq3
ATGGTATTTAAA

我的代码是这样的：

f = open("a.txt","r")
seqlist = []
for line in f.readlines():
  line = line.strip("\n")
  if line.startswith(">"):
    print(line)
  elif line.startswith("A") or line.startswith("T") or line.startswith("C") or line.startswith("G"):
    seq = line
    y = 0
    for y in range(2, len(seq), 3):
      x = seq[y]
      print(x)

现在我可以拿到第三名的ATCG了，我想把它放在一个列表中。

然后我可以数 ATCG。

但我不知道如何将它放在一个列表中。并得到以下结果。

seq1 A=3 T=3 C=1 G=1
seq2 A=? T=? C=? G=?
seq3 A=? T=? C=? G=?

非常感谢您对我的帮助。

【问题讨论】：

“计算第三位 ATCG 的数量”对我来说并不完全清楚。

我在我的问题中举了一个例子。我要算seq1~3中的第三名ATCG。非常感谢。

例如在 seq3(AT"G"GT"A"TT"T"AA"A") 中。我想数 3、6、9、12 位 ATCG 号码。所以在 seq3 A=2 T=1 C=0 G=1.

【参考方案1】：

这是一个尽可能少地修改代码的选项：

from collections import Counter

counter = None
for line in f.readlines():
    line = line.strip("\n")
    if line.startswith(">"):
        if counter is not None:
            print(counter)
        print(line)
        counter = Counter()
    elif line.startswith("A") or line.startswith("T") or line.startswith("C") or line.startswith("G"):
        seq = line
        y = 0
        for y in range(2, len(seq), 3):
            x = seq[y]
            counter[x] += 1
print(counter)

输出：

>seq1
Counter('A': 3, 'T': 3, 'C': 1, 'G': 1)
>seq2
Counter('T': 5, 'A': 4, 'C': 2, 'G': 1)
>seq3
Counter('A': 2, 'T': 1, 'G': 1)

这是同样的事情，但整体上改进了你的代码，并更好地格式化了输出：

from collections import Counter

counter = None
bases = 'ATCG'

def print_counter():
    print(' '.join('%s=%s' % (k, counter[k]) for k in bases))

with open("a.txt", "r") as f:  # Always open files like this
    for line in f:  # no need for readlines
        line = line.strip("\n")
        if line.startswith(">"):
            if counter is not None:
                print_counter()
            print(line)
            counter = Counter()
        elif line and line[0] in bases:
            counter.update(line[2::3])
print_counter()

输出：

>seq1
A=3 T=3 C=1 G=1
>seq2
A=4 T=5 C=2 G=1
>seq3
A=2 T=1 C=0 G=1

【讨论】：

非常感谢！我想问我可以打印这样的表格A =吗？ T=? C=? G=?。谢谢！