如何找到GC含量的最高百分比，包括其ID？

Question

我在Rosalind上遇到了一个问题，它涉及输出具有最高GC含量百分比和百分比的字符串的ID，我遇到了一个问题。我能够输出样本数据ID +百分比，但我不知道如何找到最高值。

#!/bin/env python3

import sys

input = input()
file = open(input, "r")
#print(input)

gc = 0
at = 0
unknown = 0

for line in file:
  if line.startswith(">"):
    if (gc + at) > 0:
     total = gc + at
     percentage = float(gc)/float(total) * 100
     result = percentage
     print(seq_id, result)
    seq_id = line.strip()
    gc = at = unknown = 0
  else:
   nuc_str = list(line.strip())
   for n in nuc_str:
    if n == "G" or n == "g" or n == "C" or n == "c":
     gc += 1.0
    elif n == "A" or n == "a" or n == "T" or n == "t":
     at += 1.0
    else:
     unknown += 1.0
total = gc + at
percentage = float(gc)/float(total) * 100
result = percentage
print(seq_id, result)

所需的输出是：

Rosalind_0808
60.919540

我得到的输出是：

>Rosalind_6404 53.75
>Rosalind_5959 53.57142857142857
>Rosalind_0808 60.91954022988506

输入文件是从rosalind获取的示例数据fasta文件：

>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT

我不想再用传记片了。

Answer 1

有一些优化可以用来使代码更好、更简洁。例如，使用list of tuples，您可以将序列头与它们的GC内容关联起来。

所以你的循环看起来是：

# Create an empty list to collect results
results = []
for line in file:
    if line.startswith(">"):
        if (gc + at) > 0:
            total = gc + at
            percentage = float(gc)/float(total) * 100
            result = percentage
            # Add results to the list here
            results.append((seq_id, result))
            seq_id = line.strip()
            gc = at = unknown = 0
    else:
        nuc_str = list(line.strip())
        for n in nuc_str:
            if n == "G" or n == "g" or n == "C" or n == "c":
                gc += 1.0
            elif n == "A" or n == "a" or n == "T" or n == "t":
                at += 1.0
            else:
                unknown += 1.0

运行循环后，results如下所示：

[
 ("Rosalind_6404", 53.75), 
 ("Rosalind_5959", 53.57142857142857), 
 ("Rosalind_0808", 60.91954022988506),
]

然后，收集具有最高GC含量的序列，我们可以使用itemgetter。

from operator import itemgetter
result = max(results, key=itemgetter(1))

使用key参数，我们告诉max查找在元组中的1位置具有最高值的列表项。

result看起来是这样的：

("Rosalind_0808", 60.91954022988506)

要以所需的格式输出结果，我们可以使用format

 print("{}\n{:.6f}".format(*result))

其中产出：

Rosalind_0808
60.919540

其他优化：

• 使用with语句打开文件，这样做的好处是在处理后正确关闭文件。
• 计数GC，并使用序列长度(来自len)作为除数
• 使用lower()代替大写和小写的比较

我写了我自己的版本，看上去如下：

from operator import itemgetter

res = []
with open("test.fa", "r") as fh:
    for line in fh:
        line = line.strip()
        if line[0] in ">":
            header = line[1:]
            seq = next(fh).strip()
            gc = sum(1 for base in seq if base.lower() in "gc")
            res.append((header, gc / len(seq) * 100))

result = max(res, key=itemgetter(1))
print("{}\n{:.6f}".format(*result))