如何从序列中删除像'>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\n'
这样的id?
我有这个代码:
with open('sequence.fasta', 'r') as f :
while True:
line1=f.readline()
line2=f.readline()
line3=f.readline()
if not line3:
break
fct([line1[i:i+100] for i in range(0, len(line1), 100)])
fct([line2[i:i+100] for i in range(0, len(line2), 100)])
fct([line3[i:i+100] for i in range(0, len(line3), 100)])
输出:
['>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\n']
['CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG\n']
['AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG\n']
['CCGCCTCGGGAGCGTCCATGGCGGGTTTGAACCTCTAGCCCGGCGCAGTTTGGGCGCCAAGCCATATGAA\n']
['AGCATCACCGGCGAATGGCATTGTCTTCCCCAAAACCCGGAGCGGCGGCGTGCTGTCGCGTGCCCAATGA\n']
['ATTTTGATGACTCTCGCAAACGGGAATCTTGGCTCTTTGCATCGGATGGAAGGACGCAGCGAAATGCGAT\n']
['AAGTGGTGTGAATTGCAAGATCCCGTGAACCATCGAGTCTTTTGAACGCAAGTTGCGCCCGAGGCCATCA\n']
['GGCTAAGGGCACGCCTGCTTGGGCGTCGCGCTTCGTCTCTCTCCTGCCAATGCTTGCCCGGCATACAGCC\n']
['AGGCCGGCGTGGTGCGGATGTGAAAGATTGGCCCCTTGTGCCTAGGTGCGGCGGGTCCAAGAGCTGGTGT\n']
['TTTGATGGCCCGGAACCCGGCAAGAGGTGGACGGATGCTGGCAGCAGCTGCCGTGCGAATCCCCCATGTT\n']
['GTCGTGCTTGTCGGACAGGCAGGAGAACCCTTCCGAACCCCAATGGAGGGCGGTTGACCGCCATTCGGAT\n']
['GTGACCCCAGGTCAGGCGGGGGCACCCGCTGAGTTTACGC\n']
['\n']
...
我的功能是:
def fct(input_string):
code={"a":0,"c":1,"g":2,"t":3}
p=[code[i] for i in input_string]
n=len(input_string)
c=0
for i, n in enumerate(range(n, 0, -1)):
c +=p[i]*(4**(n-1))
return c+1
fct()
从字符串返回一个整数。例如,ACT
给出8
,即:我的函数必须采用输入字符串序列只包含以下基数A,C,G,T
但是当我使用我的功能时它会给出:
KeyError: '>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\n'
我尝试通过剥离行来删除id,从>
开始并将其余部分写入文本文件中,因此,我的文本文件output.txt
只包含没有id的序列,但是当我使用我的函数fct时,我发现了同样的错误:
KeyError: 'CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG\n'
我能做什么?
我在您的代码中看到两个主要问题:您在解析FASTA序列时遇到问题,并且您的函数没有正确迭代每个序列。
我建议使用优秀的Biopython包吗?它内置了极好的FASTA支持(读写)(参见Sequences中的Tutorial)。
要解析FASTA文件中的序列:
for seq_record in SeqIO.parse("seqs.fasta", "fasta"):
print record.description # gi|2765658|emb|Z78533.1...
print record.seq # a Seq object, call str() to get a simple string
>>> print record.id
'gi|2765658|emb|Z78533.1|CIZ78533'
>>> print record.description
'gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA'
>>> print record.seq
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet())
>>> print str(record.seq)
'CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACC' #(truncated)
在你的代码中,你有一个传递给fct()
的字符串列表(input_string
实际上不是一个字符串,而是一个字符串列表)。解决方案只是构建一个输入字符串,并迭代它。
fct
:c
立即返回。p
时,为什么还要构造code
呢?n
循环中使用它作为变量名称来编写序列的长度(for
)修改后的代码(使用适当的PEP 8格式化),并将变量重命名为更清楚它们的意思(仍然不知道c
应该是什么):
from Bio import SeqIO
def dna_seq_score(dna_seq):
nucleotide_code = {"A": 0, "C": 1, "G": 2, "T": 3}
c = 0
for i, k in enumerate(range(len(dna_seq), 0, -1)):
nucleotide = dna_seq[i]
code_num = nucleotide_code[nucleotide]
c += code_num * (4 ** (k - 1))
return c + 1
for record in SeqIO.parse("test.fasta", "fasta"):
dna_seq_score(record.seq)