獣医師が必死にバイオインフォマティクスPythonによる実践レシピをwin10で動かす（その6）

Biopythonで配列を扱う

ようやく時間が取れました
では前回取り出した犬のNGFの配列を基にあれこれいじってみる
まずは配列をNCBIから呼び出す

from Bio import Entrez, Seq, SeqIO
from Bio.Alphabet import IUPAC

Entrez.email = "put@your_email.here" 
hdl = Entrez.efetch(db='nucleotide', id=['XM_038690419'], rettype='fasta') 
seq = SeqIO.read(hdl, 'fasta')

のっけからここを調べていて、Bio.Alphabetパッケージが廃止されることを知る
本が発売されたのが2018年だけど翻訳されたのは最近だからわかるけども、日本語で勉強をする弊害をひしひしと感じる
もはやこの時点で暗澹たる気持ちなのだけど、気にせず突貫
とりあえず取った情報を確認

print("sequence_length: "+ str(len(seq)) + " bp")
seq

意外と小さいmRNAだなーと思う
とは言え非モデル動物の場合RNA-seqによって出てきたものから推定した情報が多いので、しょっちゅう覆されるのですが
とりあえず使いやすくするために情報を格納

w_hdl = open('example.fasta', 'w')
SeqIO.write([seq], w_hdl, 'fasta')
w_hdl.close()

繰り返し開け閉めしているとメモリが圧迫されるという基礎知識はあるけど、実際にそんな状態に遭遇したことは無い
誰か実感する方法を教えてプリーズ

recs = SeqIO.parse('example.fasta', 'fasta')
for rec in recs:
    seq = rec.seq
    print(rec.description)
    print(seq[:10])
    print(seq.alphabet)

これで情報は見やすくなった
しかし、このSingleLetterAlphabet()って情報はいるのだろうか？わからん

続いて、IUPACのATCGだけの塩基を表示させてみる

seq = Seq.Seq(str(seq), IUPAC.unambiguous_dna)
seq

基本的には同じ情報が表示されているだけだけど、ATGC以外が入っていると表示されないようになった
IUPAC.unambiguous_dnaとIUPAC.ambiguous_dnaを使って多型のある塩基を処理しようと思ったこともあったけど、
じつはprimer設計をするならば全てNに変換してしまった方が早く楽
seq = str.maketrans("URYMKSWHBVD", "NNNNNNNNNNN")
こんな感じに書くとIUPACコードを全部Nに変換可能

続いて転写したRNAを取得してみる

print((seq[:12], seq[-50:]))
rna = seq.transcribe()
rna

情報もIUPACUnambiguousRNA()となっているのでAUCGになっていると思われる
更にこれをタンパク質に変換してみる

prot = seq.translate()
prot

怒られた
C:\Users\anaconda3\envs\bioinfomatics\lib\site-packages\Bio\Seq.py:2309: BiopythonWarning: Partial codon, len(sequence) not a multiple of three. Explicitly trim the sequence or add trailing N before translation. This may become an error in future. BiopythonWarning)
3の倍数じゃねーよ、ちゃんとORFだけにしたのかと怒られていて、開始コドンを考えていなかったことに思い当たる
というわけで開始コドンの情報を取りに行く

recs = list(SeqIO.parse(hdl, 'gb'))
for feature in rec.features:
    if feature.type == 'CDS':
        print("CDS_Location: " + str(feature.location))

CDS_Location: [220:946](+)なので、これを組み込んで再度挑戦

cds_seq = seq[220:946]
prot = cds_seq.translate()
prot

無事に怒られずに終始コドンまで到達！
今更思ったけど、この場合のCDSの位置は-1とかしなくても良いみたい
これがBiopythonを使うメリットなのかもと思い、今日を終える