我有一个包含蛋白质序列(200 个序列)的文本文件,如下所示。
>ptn1
AAGHM
>ptn2
MGLKKRR
我需要为序列的每个字符提供以下值,并且必须找到每个序列的平均值。
A= 0.2, G= 0.5, L=0.14, M= 0.70, R= 0.55, C=0.48, H= 1.00 , K=0.4
期望的输出
ptn1 - 0.52
ptn2 - 0.462
如何使用 awk 或 python 做到这一点?
您的建议将不胜感激
我有一个包含蛋白质序列(200 个序列)的文本文件,如下所示。
>ptn1
AAGHM
>ptn2
MGLKKRR
我需要为序列的每个字符提供以下值,并且必须找到每个序列的平均值。
A= 0.2, G= 0.5, L=0.14, M= 0.70, R= 0.55, C=0.48, H= 1.00 , K=0.4
期望的输出
ptn1 - 0.52
ptn2 - 0.462
如何使用 awk 或 python 做到这一点?
您的建议将不胜感激
def avg(sequence):
v= {'A': 0.2, 'C': 0.48, 'R': 0.55, 'G': 0.5, 'H': 1.0,
'K': 0.4, 'M': 0.7, 'L': 0.14}
return sum(v[x] for x in sequence) / len(sequence)
avg("AAGHM") # => 0.5199999999999999
avg("MGLKKRR" # => 0.46285714285714274
FS=""
http://www.gnu.org/software/gawk/manual/html_node/Single-Character-Fields.html#Single-Character-Fields需要 gawk
用法:
awk -f foo.awk foo.txt
BEGIN {
FS=""
k["A"]=0.2; k["G"]=0.5; k["L"]=0.14; k["M"]=0.70
k["R"]=0.55; k["C"]=0.48; k["H"]=1.00; k["K"]=0.4
}
/^>/{
$1=""
name=$0
next
}
{
s=0
for (i=1; i<=NF; i++) {
s+=k[$(i)]
}
printf "%s - %.3f\n", name, s/NF
}
您的序列文件似乎是 FASTA 格式。您应该使用专门适用于处理序列的工具,而不是使用awk
.
如果您的文件格式发生变化怎么办?如果您想根据格式规范提取有意义的数据怎么办?已经构建的解析器非常适合这些问题。
我喜欢Biopython:
from Bio import SeqIO
records = SeqIO.parse("sequences.fasta", "fasta")
# Function borrowed from Tichodrama's answer
def avg(sequence):
v= {'A': 0.2, 'C': 0.48, 'R': 0.55, 'G': 0.5, 'H': 1.0,
'K': 0.4, 'M': 0.7, 'L': 0.14}
return sum(v[x] for x in sequence) / len(sequence)
for record in records:
print("Score for {}: {:.2f}".format(record.id, avg(record.seq)))
请注意,record
上述每个SeqRecord
对象都是一个包含有用信息的对象,例如解析的定义信息和序列字母表(例如,您可以区分蛋白质序列和 DNA 序列)。
mapping = {
'A': 0.2,
'G': 0.5,
'L': 0.14,
'M': 0.7,
'R': 0.55,
'C': 0.48,
'H': 1.0,
'K': 0.4}
def process(seqs):
i = 0
while i < len(seqs) - 1:
print seqs[i].strip()[1:], '-', avg(seqs[i+1])
i += 2
def avg(seq):
return sum(mapping[v] for v in seq) / len(seq)
# assume the file is data.txt
process(open('data.txt').readlines())