提问人：RichardAzolla 提问时间：11/3/2015 最后编辑：Timur ShtatlandRichardAzolla 更新时间：11/18/2023 访问量：562
如何在fasta文件中查找序列长度？
How to find sequence lengths in a fasta file?
问：
我目前正在处理一个文件（一个文本文件），其中包含 DNA 提取序列（重叠群）列表，每个序列都有一个标题，后跟核苷酸行，这是该重叠群的核苷酸长度。有 120 个重叠群，每个条目都用以“>”开头的行标记，以表示序列信息。在这条线之后，给出该序列的核苷酸长度。fasta
例：
>gi|571136972|ref|XM_006625214.1| Plasmodium chabaudi chabaudi small subunit ribosomal protein 5 (Rps5) (rps5) mRNA, complete cds
ATGAGAAATATTTTATTAAAGAAAAAATTATATAATAGTAAAAATATTTATATTTTATATTATATTTTAATAATATTTAAAAGTATTTTTATTATTTTATTTAATAGTAAATATAATGTGAATTATTATTTATATAATAAAATTTATAATTTATTTATTATATATATAAAATTATATTATATTATAAATAATATATATTATAATAATAATTATTATTATATATATAATATGAATTATATA
TATTTTTATATTTATAAATATAATAGTTTAAATAATA
>gi|571136996|ref|XM_006625226.1| Plasmodium chabaudi chabaudi small subunit ribosomal protein 2 (Rps2) (rps2) mRNA, complete cds
ATGTTTATTACATTTAAAGATTTATTAAAATCTAAAATATATATAGGAAATAATTATAAAAATATTTATATTAATAATTATAAATTTATATATAAAATAAAATATAATTATTGTATTTTAAATTTTACATTAATTATATTATATTTATATAAATTATATTTATATATTTATAATATATCTATATTTAATAATAAAATTTTATTTATTATTAATAATAATTTAATTACAAATTTAATTATT
AATATATGTAATTTAACTAATAATTTTTATATTATTA
我想做的是列出每个重叠群。我的问题是，我不知道告诉 Python 所需的语法：
- 找到以“>”开头的行之后的行
- 计算该序列各行中的所有字符
- 将值返回到所有重叠群值的列表（给出每个重叠群长度列表的列表，即 126、300、25...）
- 确保最后一个重叠群（没有“>”表示其结束）被计算在内。
我想要一个整数列表，这样我就可以计算重叠群的平均长度、标准差、酷基因方程等。
答：
0赞
Robᵩ
11/3/2015
#1
这可能对您有用：它会在包含以下内容的行后面的行中打印 ACGT 的数量：>
import re
with open("input.txt") as input_file:
data = input_file.read()
data = re.split(r">.*", data)[1:]
data = [sum(1 for ch in datum if ch in 'ACGT') for datum in data]
print(data)
评论
0赞
RichardAzolla
11/4/2015
这对我有用，我已将其纳入我的程序，谢谢。
3赞
jpalmer
11/3/2015
#2
不要重新发明轮子，按照 Martin 的建议使用 biopython。这是您的开始，它将序列 ID 和长度打印到终端。您可以使用 pip 安装 biopython，即pip install biopython
from Bio import SeqIO
import sys
FileIn = sys.argv[1]
handle = open(FileIn, 'rU')
SeqRecords = SeqIO.parse(handle, 'fasta')
for record in SeqRecords: #loop through each fasta entry
length = len(record.seq) #get sequence length
print "%s: %i bp" % (record.id, length) #print sequence ID: seq length
或者，您可以将结果存储在字典中：
handle = open(FileIn, 'rU')
sequence_lengths = {}
SeqRecords = SeqIO.parse(handle, 'fasta')
for record in SeqRecords: #loop through each fasta entry
length = len(record.seq) #get sequence length
sequence_lengths[record.id] = length
#access dictionary outside of loop
print sequence_lengths
0赞
RichardAzolla
11/4/2015
#3
感谢所有的帮助。我看过生物蟒蛇的东西，很高兴能理解它并融入它。这个作业的总体目标是教我如何理解 python，而不是直接找到解决方案，或者至少如果我找到解决方案，我必须能够用我自己的话来解释它。
无论如何，我已经创建了一个包含该元素以及其他元素的代码。我还有几件事要做，如果我感到困惑，我会回来问。
这是我的第一个工作代码，除了直接与我的主管合作或我制作和理解的教程（呜呜！
import re
with open("COPYFORTESTINGplastid.1.rna.fna") as fasta:
contigs = 0
for line in fasta:
if line.strip().startswith('>'):
contigs = contigs + 1
with open("COPYFORTESTINGplastid.1.rna.fna") as fasta:
data = fasta.read()
data = re.split(r">.*", data)[1:]
data = [sum(1 for ch in datum if ch in 'ACGT') for datum in data]
print "Total number of contigs: %s" %contigs
total_contigs = sum(data)
N50 = sum(data)/2
print "number used to determine N50 = %s" %N50
average = 0
total = 0
for n in data:
total = total + n
mean = total / len(data)
print "mean length of contigs: %s" %mean
print "total nucleotides in fasta = %s" %total_contigs
#print "list of contigs by length: %s" %sorted([data])
l = data
l.sort(reverse = True)
print "list of contigs by length: %s" %l
这符合我的要求，但如果您有任何意见或建议，我很想听听。
接下来，用这个甜蜜的清单确定N50。再次感谢！
评论
0赞
BioGeek
11/6/2015
代码审查并不是 Stackoverflow 的真正目的，但如果你把你的代码发布在 codereview.stackexchange.com 上，我会在那里给你我的代码评论。
0赞
plantpathbreeding
11/10/2020
#4
我创建了一个函数来计算 N50，它似乎工作得很好。我可以解析命令行并通过程序运行任何 .fa 文件
def calc_n50(array):
array.sort(reverse = True)
n50 = 0 #sums lengths
n = 0 #n50 sequence
half = sum(array)/2
for val in array:
n50 += val
if n50 >= half:
n = val
break #breaks loop when condition is met
print "N50 is",n
评论