Biopython Tutorial and Cookbook 3.2[和訳］

今回はパソコンを使うペンギンのイラストにしてみました。
かわいい。
ただそれだけでこの上ない理由です。癒しって重要ですものね。
和んだところで本題へ行きましょう。

第3章2節　Sequences act like strings

多くの方法で、私たちはSeqオブジェクトから長さを取得したり要素を反復処理するなど通常のPython文のように扱えます。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("GATCG", IUPAC.unambiguous_dna)
>>> for index, letter in enumerate(my_seq):
        print("%i %s" % (index, letter))
0 G
1 A
2 T
3 C
4 G

>>> print(len(my_seq))
5

Pythonと同じ方法で文字列中の要素にアクセスすることができます (この時、最初の要素を0と数えることに注意)。

>>> print(my_seq[0]) #first letter
G
>>> print(my_seq[2]) #third letter
T
>>> print(my_seq[-1]) #last letter
G

Seqオブジェクトは.count()メソッドを持つ。これはPythonと同様に重複しないカウントを与える。

>>> from Bio.Seq import Seq
>>> "AAAA".count("AA")
2

>>> Seq("AAAA").count("AA")
2

生物学的用途のために、重複 (オーバーラップ)したカウントが必要になる場合があるかもしれません (上の例だと3)。一文字で検索する場合には違いはありません。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
>>> len(my_seq)
32
>>> my_seq.count("G")
9
>>> 100 * float(my_seq.count("G") + my_seq.count("C")) / len(my_seq)
46.875

上のコードによってGC%を計算可能ですが、Bio.SeqUtilsモジュールにはGC関数が既にビルトされています。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> from Bio.SeqUtils import GC
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
>>> GC(my_seq)
46.875

Bio.SeqUtils.GC()関数は自動的に配列と多義塩基“S” (G or C)が処理されることに注意してください。
また、Pythonの文字列と同じように、Seqオブジェクトは読み取り専用です。もし、点突然変異のシミュレーションなどで配列を書き換える必要がある場合は、3章12節のMutableSeqオブジェクトを見てください。

今回のまとめ

今回はstringの扱い方講座ですね。stringの扱い方を解説始めるとそれだけで一回分になってしまう気がするので割愛。

配列解析やったことないのでわからないのですが、配列中にある特定の一塩基だけをピックアップする機会ってあるのでしょうか…？それも場所指定で。それこそ一塩基置換の時くらいしか状況が想像つかないですね。別の方法で置換箇所を特定するのでしょうか。

ふと気が付いたのですが、"生物学においてオーバーラップしたカウントが必要"って言っておきながらこの中で解説してないですね。自分で調べろということか。解説書の意味。というわけで別項目にてまとめます。
以下今回のまとめ

・文字列の基本的な扱い方は普通のPythonと同じ。
・Seqオブジェクト内は.count()でカウントできる(重複を許さない)。
・GC含量はGC関数で一発で計算できる。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> from Bio.SeqUtils import GC

#.count()の使い方。
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
>>> my_seq.count("GG")
2

#GC function
>>> GC(my_seq)
46.875

すごい、3行でまとまった。
GC関数は便利ですが、ATGC以外の文字が含まれるときは注意が必要そうです。特にMとかKとかRとかYとか。この検証も下にまとめます。

文字列をオーバーラップを許してカウントする場合

説明書に説明されていなかったオーバーラップ時のカウント方法は以下の通りです。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> from Bio.SeqUtils import GC

>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)

#オーバーラップ無し
>>> my_seq.count("GG")
2

#オーバーラップ有り
>>> my_seq.count_overlap("GG")
3

.count_overlap()という、そのままなメソッドがありました。以下のURLが元になりますので疑り深い方はそちらも参照してください。
https://github.com/biopython/biopython/commit/97709cc7a4b8591e794b442796d41e4f7b855a0f

GC functionの挙動

＝クラスを参照している？＝
3.1は各クラスについて紹介されていましたが、Tutorialの中ではDNA用のクラスしか書かれていませんでした。当たり前なのですが。というわけで、GC関数にかけるSeqオブジェクトがタンパク用のクラスを持っていた場合どうなるのか試してみました。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> from Bio.SeqUtils import GC

#タンパク用クラスでGC関数を使用
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC",  IUPAC.protein)
>>> GC(my_seq)
46.875

普通に結果を吐きますね

#一応クラス確認
>>> my_seq.alphabet
IUPACProtein()

クラスはタンパク用になっていますが、関係なく計算結果吐くようです。

＝多義シンボルが含まれる場合は？＝
塩基が"S"ならGCどちらだったとしてもGC含量に影響ありませんが、MやK等々Strong組 (G or C)とWeak組 (A or T)を跨ぐ塩基の場合どうなるのかやってみました。

>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> from Bio.SeqUtils import GC

#上の配列を少し弄った
                  *                              *
                  GATCGATGGGCCTATATAGGATCGAAAATCGC
>>> my_seq = Seq("KATCGATGGGCCTATATAGGATCGAAAATCGS", IUPAC.unambiguous_dna)
>>> GC(my_seq)
43.75

下がりましたね。
挙動から察するに、GC関数はSをGまたはCに置換してから

100 * float(my_seq.count("G") + my_seq.count("C")) / len(my_seq)

を走らせてるだけな印象です。

GC関数使用時の注意事項

・GC関数はクラスに関係なく動く。
・多義シンボルが含まれる場合教えてくれるなどの優しい機能はない。

ということでした。
配列がSの時はきちんと考慮してくれているのにそれ以外のシンボルが使われている場合には
適切でないデータでも数字が出てきてしまうので注意して使いましょうということですね。

おわりに

GC含量の計算という使えるワードが出てき始めたことでオラァワクワクしてきたぞ。そして、書きながら検証してると結構な時間を使いますね。

第3章をすべてまとめ終わった段階で、和訳無しの説明書的なものを作ろうと思ってます。使う人がいるかはわからないけど。
これ全部きれいに読み終わるのいつ頃になるのだろう…（遠い目

ではまた次回