3

以下に示すように、タンパク質配列（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でこれを行うにはどうすればよいですか?

あなたの提案をいただければ幸いです

4

4 に答える 4

6 
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
于 2012-07-23T06:51:45.807 に答える
5

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
}
于 2012-07-23T07:20:14.860 に答える
1

シーケンス ファイルは FASTA 形式のようです。などのツールを使用したカスタム解析ではなく、シーケンスの操作に特に適したツールを使用する必要がありますawk。

ファイル形式が変わったら？フォーマット仕様に従って意味のあるデータを抽出したい場合はどうすればよいでしょうか? 構築済みのパーサーは、これらの質問に最適です。

私はバイオパイソンが好きです：

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 配列を区別します)。

于 2012-07-27T11:42:02.313 に答える
0 
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())
于 2012-07-23T06:58:41.063 に答える