如何找到最高百分比的 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

我尽量不使用 biopython。

Answer 1

您可以使用一些优化来使您的代码更好、更简洁。例如，使用 tuples 的 list 可以使序列头与其 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))

如何找到最高百分比的 GC 内容及其 ID？

How to find the highest percentage of GC content including its ID?

python

bioinformatics