Pythonバイオ?　Pythonバイオ/ツール?
3908　　　2019-02-21 (木) 16:56:44

配列の切出し・加工 †

↑

BioPythonのSeq型 †

Seq型はDNA、RNA、タンパクなどの配列を扱うのに便利な型です。まずは例を見てみましょう。

from Bio.Alphabet import IUPAC
from Bio.Seq import Seq

dna_seq = Seq('ATGAAACGCATTAGCACCACC', IUPAC.ambiguous_dna)
dna_seq
## Seq('ATGAAACGCATTAGCACCACC', IUPACAmbiguousDNA())

rna_seq = Seq('AUGAAACGCAUUAGCACCACC', IUPAC.ambiguous_rna)
rna_seq
## Seq('AUGAAACGCAUUAGCACCACC', IUPACAmbiguousRNA())  

aa_seq = Seq("MKRISTT", IUPAC.protein)
aa_seq
## Seq('MKRISTT', IUPACProtein())

Seq型のデータを作るには、第１引数に配列を表す文字列、第２引数に配列のタイプを指定します。配列のタイプはモジュールIUPACの説明にあるように、基本的にDNA、RNA、タンパクの３種類で、Ambiguous（アルファベットなら何でも可）かUnambiguous（GATCのみ）かが指定できます。但しSeq型変数を作るときにアルファベットをチェックしているのではなく、別途プログラムでチェックする仕組みが用意されています。

from Bio.Alphabet import _verify_alphabet
wrong_seq = Seq('PQR', IUPAC.unambiguous_dna)
_verify_alphabet(wrong_seq)
## False
↑

位置を指定して切出し †

手元システムに配列ファイルがある場合はファイル読み出しによって取り込んだのち、切り出しの操作をする。

Fasta形式参照: のファイルを読み込むには、下記のようにSeqIO.parseが便利である。

from Bio import SeqIO
for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

この例ではFastaファイル中に複数の配列データが含まれているので、forループによって１つずつ見て、そのid部分、配列自体、配列の長さを表示している。

複数の配列を含むFastaファイルからk番目の配列を取出すときは、次のようにlist関数を使って一旦リストにしたのち、k番目（ここでは2番目）の配列をseq.seqとして取出す。seq.seqとして取出したものは、Pythonのリスト（文字列）と同様に、インデックスやインデックス範囲を使って配列の一部を切り出すことができる。

from Bio import SeqIO
seq_records = list(SeqIO.parse("ls_orchid.fasta", "fasta"))
seq = seq_records[2]
print(seq.seq)
print(seq.seq[3:6])

Pythonの文字列の扱いが使える。たとえばさまざまなスライシング、２つの配列の結合

seq[::-1]

文字列と見なして、大文字・小文字を変換するupperやlowerが使える。

公共データベースから配列を取り込む場合は、取込みの手順によって取り込んだのち、切り出しの操作をする。

↑

配列の加工 †

Seq型は変更ができない型（Pythonでimmutableと呼ばれる型）になっている。たとえば、

from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
my_seq[5] = 'G'

とすると、

TypeError: 'Seq' object does not support item assignment

のようなエラーとなる。

配列を書き変えたい場合、変更可能な型（mutable sequence, MutableSeqオブジェクト）に変更するか、さもなければ一旦内容を文字列に読み出してそれを文字列として加工した上で再度Seq型に作る必要がある。

変更可能なMutableSeq型に変換するには、Seq型のオブジェクトをtomutable()でMutableSeqオブジェクトに変換するか、または直接にMutableSeqオブジェクトとして生成する。

mut_seq = my_seq.tomutable()
mut_seq
### MutableSeq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

または直接にMutableSeqオブジェクトとして生成する。

from Bio.Seq import MutableSeq
from Bio.Alphabet import IUPAC
mut_seq = MutableSeq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPAC.unambiguous_dna)
mut_seq
### MutableSeq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

このようにして得られたMutatlbeSeqオブジェクトは、要素の変更をすることができる。

mut_seq[5] = 'G'
mut_seq
### MutableSeq('GATCGGTGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

では第６塩基がAからGに書き換わっている。また

mut_seq.remove('T')
mut_seq
### MutableSeq('GACGGTGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

のようにすると最初に出現するTが取り除かれる。

なお、MutableSeqから逆に変更できないSeq型に戻すにはtoseq()を用いて

new_seq = mut_seq.toseq()
new_seq
### Seq('GACGGTGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

とすればよい。

なお、先述したようにstrを用いて文字列に書き換えてから、文字列の関数を使って変更を加え、Seqを使ってSeq型に戻すということも可能である。

my_string = str(my_seq)   # my_seqを文字列に変換する
my_string
### 'GATCGATGGGCCTATATAGGATCGAAAATCGC'
new_string = my_string.replace('T', '', 1) # 'T'を''に１回だけ置換する（最後の引数1は回数指定）
new_string
### 'GACGATGGGCCTATATAGGATCGAAAATCGC'  # 回数指定しないとすべてのTを置換するので注意
new_seq = Seq(new_string, IUPAC.unambiguous_dna)  # 再度Seq型に戻す
new_seq
Seq('GACGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())
↑

ID（染色体やdescription）を指定して切出し †

acc_name = 'Z78444.1' 
for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    if acc_name in seq_record.id:
        print(seq_record.id)
        print(repr(seq_record.seq))
        print(len(seq_record))
↑

相補鎖（complement）・逆鎖（reverse）・逆相補鎖（reverse complement）の生成 †

from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
my_seq
### Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

my_seq.complement()   # 相補鎖
### Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG', IUPACUnambiguousDNA())

my_seq[::-1]      # 逆鎖
## Seq('CGCTAAAAGCTAGGATATATCCGGGTAGCTAG', IUPACUnambiguousDNA())

my_seq.reverse_complement()    # 逆相補鎖
### Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC', IUPACUnambiguousDNA())
↑

配列の転写（transcription） †

from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
my_seq
### Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

messenger_rna = my_seq.transcribe()
messenger_rna
### Seq('GAUCGAUGGGCCUAUAUAGGAUCGAAAAUCGC', IUPACUnambiguousRNA())
↑

配列の翻訳（translation） †

from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATC", IUPAC.unambiguous_dna)
my_seq
### Seq('GATCGATGGGCCTATATAGGATCGAAAATC', IUPACUnambiguousDNA())

translated_protein = my_seq.translate()
translated_protein
### Seq('DRWAYIGSKI', IUPACProtein())

翻訳表（Translation Table）はNCBIの定義する表が用意されている。標準の表（デフォルト）はStandard Codeであるが、たとえばVertebrate Mitochondrial Codeを使いたい場合には、translateの引数にtable='Vertebrate Mitochondrial'を指定する。

translated_protein = my_seq.translate(table='Vertebrate Mitochondrial')
translated_protein
### Seq('DRWAYMGSKI', IUPACProtein())

また、終止コドン（stop codon）が現れたときに翻訳を停止するようにするには、引数にto_stop=Trueを指定する。

my_seq = Seq('GATCGATGGGCCTAAATAGGATCGAAAATC', IUPAC.unambiguous_dna)
translated_protein = my_seq.translate(to_stop=True)
translated_protein
Seq('DRWA', IUPACProtein())
↑

k-mer（k連続塩基の出現頻度解析） †

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-02-21 (木) 16:56:44 (1320d)