MEGAHITのPaired-endアセンブリが怪しい
MEGAHITは言わずと知れた、de Bruijn graphを利用したメタゲノム用の高速アセンブラである。
MEGAHITはFASTQをインプットとし、Paired-endとSingle-endの両方に対応している。
Input options that can be specified for multiple times (supporting plain text and gz/bz2 extensions)
-1 <pe1> comma-separated list of fasta/q paired-end #1 files, paired with files in <pe2>
-2 <pe2> comma-separated list of fasta/q paired-end #2 files, paired with files in <pe1>
--12 <pe12> comma-separated list of interleaved fasta/q paired-end files
-r/--read <se> comma-separated list of fasta/q single-end files
Paired-endの場合はファイルの形式に応じて-1, -2か-12で指定する。Single-endは-r。
ミトコンドリアの完全ゲノムを得るために、ミトコンドリアにマップされたリードのみを抽出し、いくつかのアセンブラでテストをしていたときのこと。
(MEGAHIT v1.1.1-2-g02102e1)
公開データの中には、Paired-endとSingle-endの両方があったので、Paired-endのデータは-1,-2で、Single-endは-rで指定していた。
すると、なぜかSingle-endのときにのみ完全ゲノム1本がきれいに出力された。Single-endのシーケンサはIllumina GA-IIで、長さも100bpであり、決して他に秀でてうまくアセンブルできるとは思えなかった。
そこで、Paired-endのものもForwardだけを用いて-rでアセンブルしてみたところ、(孤立した配列はいくつか混じっていたものの)完全長のコンティグが生成できた。
MEGAHITはPairの情報をそこまで利用していないという認識があったので、この結果には驚いた。他のバージョンは試していないが、Singleにしたほうが断然良い結果であった。メタゲノムデータに対する場合はわからないが、Singleにしたほうが良い結果が出る場合が他にもあるかも知れないので注意したい。
Bowtie2で許容編集距離を設定する
Bowtie2のoptionは複雑で、とりあえずデフォルトで回してからSAMをパースしている人は多いかもしれない。
マニュアルを良く読むとスコアをきちんと設定すれば色々設定できる。
Bowtie 2: Manual
特に以下の部分。
Scoring options
--mp MX,MN
Sets the maximum (MX) and minimum (MN) mismatch penalties, both integers. A number less than or equal to MX and greater than or equal to MN is subtracted from the alignment score for each position where a read character aligns to a reference character, the characters do not match, and neither is an N. If --ignore-quals is specified, the number subtracted quals MX. Otherwise, the number subtracted is MN + floor( (MX-MN)(MIN(Q, 40.0)/40.0) ) where Q is the Phred quality value. Default: MX = 6, MN = 2.--np
Sets penalty for positions where the read, reference, or both, contain an ambiguous character such as N. Default: 1.--rdg
,
Sets the read gap open () and extend ( ) penalties. A read gap of length N gets a penalty of + N * . Default: 5, 3. --rfg
,
Sets the reference gap open () and extend ( ) penalties. A reference gap of length N gets a penalty of + N * . Default: 5, 3. --score-min
Sets a function governing the minimum alignment score needed for an alignment to be considered "valid" (i.e. good enough to report). This is a function of read length. For instance, specifying L,0,-0.6 sets the minimum-score function f to f(x) = 0 + -0.6 * x, where x is the read length. See also: setting function options. The default in --end-to-end mode is L,-0.6,-0.6 and the default in --local mode is G,20,8.
日本語で要約すると
--mp: ミスマッチペナルティ(最大,最小)、リードのクオリティでスコアがバラける。
--np: N(曖昧塩基)のペナルティ、Bowtie2では他の縮退塩基もNとして扱われる。
--rdg: リードのギャップペナルティ(開始,延長それぞれについて規定)
--rfg: 参照配列のギャップペナルティ(開始,延長それぞれについて規定)
これを編集距離へを変化させるには、それぞれを1に設定すれば良い(ギャップ開始は0)。
最も煩雑なのは
--score-min
であり、これは以下のように書いてあるが、要約するとリードの長さを変数としたときの "式のタイプ, y切片,傾き" である。
Setting function options
Some Bowtie 2 options specify a function rather than an individual number or setting. In these cases the user specifies three parameters: (a) a function type F, (b) a constant term B, and (c) a coefficient A. The available function types are constant (C), linear (L), square-root (S), and natural log (G). The parameters are specified as F,B,A - that is, the function type, the constant term, and the coefficient are separated by commas with no whitespace. The constant term and coefficient may be negative and/or floating-point numbers.
For example, if the function specification is L,-0.4,-0.6, then the function defined is:
f(x) = -0.4 + -0.6 * x
If the function specification is G,1,5.4, then the function defined is:f(x) = 1.0 + 5.4 * ln(x)
See the documentation for the option in question to learn what the parameter x is for. For example, in the case if the --score-min option, the function f(x) sets the minimum alignment score necessary for an alignment to be considered valid, and x is the read length.
よって、編集距離1以内のものだけを出力するオプションは以下のようになる。
--mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "C,-1,0"
”C”については例示すら無いが一応これでちゃんと動くのは確認した。(SAMファイル内で "XM:i:" が1を超えるものがなかった。)
PythonのListとDictionaryの検索機能の速さ
"in"を使って要素の有無を判定する時、ListとDictionaryで速いのはどちらなのだろうと疑問に思ったので調べた。
with Python 2.7.11
l = [1,2,3,4,5] for i in range(10000000): if 1 in l: pass
real 0m1.017s
user 0m0.853s
sys 0m0.141s
d = {1:"", 2:"", 3:"", 4:"", 5:""}
for i in range(10000000):
if 1 in d:
pass
real 0m1.053s
user 0m0.888s
sys 0m0.143s
単純に検索するだけならList型のほうが速い。(何に使うかは知らない)
ただ、普通は検索したあとそのindexに対してカウントしたり何かしらアクションを加えるので、結局Dictionaryの方が速い。
l = [1,2,3,4,5] for i in range(10000000): if 1 in l: index = l.index(1) l[index] = l[index]
real 0m4.239s
user 0m4.074s
sys 0m0.139s
d = {1:"", 2:"", 3:"", 4:"", 5:""}
for i in range(10000000):
if 1 in d:
d[1] = d[1]
real 0m2.105s
user 0m1.938s
sys 0m0.144s
追記
要素の数が大きくなったときはDictionaryの方が早いようだった。
#d = {i:True for i in range(10000000)} d = range(10000000) for i in range(100000000): 1 in d
real 0m15.324s
user 0m13.374s
sys 0m1.961s
d = {i:True for i in range(10000000)}
#d = range(10000000)
for i in range(100000000):
1 in d
real 0m14.133s
user 0m11.922s
sys 0m2.220s
Dictionary objectはListに比べて作るのに時間がかかるため、本来ならforの回数を変化させて線形近似しなければいけないが、速さがDictionary＞Listの例が取れたのでそこは割愛した。
ちなみにDictionaryに於いて"d[1]"と"i in d"の検索速度に差はないようである。
PythonのnetworkXを使ってpathway解析をする
ある遺伝子(機能)セットと代謝ネットワークデータを用いて、上流の化合物から下流の化合物に行くpathwayがあるか調べたい。
PythonのnetworkXを使って最短経路を出すドキュメントはたくさんあったが、エッジの有無によってFlow解析についてはあまりなかったのでまとめた。
Flow解析をするには、edgeに"capacity"を設定する必要がある。
0→1→(2,3)→4 という枝分かれを1つ持つ単純なネットワークを、以下のように設定する。
>>> import networkx as nx >>> Graph = nx.DiGraph() >>> Graph.add_nodes_from([0,1,2,3,4]) >>> Graph.add_edges_from([[0,1], [1,2], [1,3], [2,4], [3,4]]) >>> Graph[0][1]["capacity"] = 1.0 >>> Graph[1][2]["capacity"] = 1.0 >>> Graph[1][3]["capacity"] = 1.0 >>> Graph[2][4]["capacity"] = 1.0 >>> Graph[3][4]["capacity"] = 1.0
edgeの"capacity"を1.0に設定することは、その遺伝子が含まれることを意味する。
今、すべてのedgeの"capacity"が1.0なので、
>>> nx.maximum_flow(Graph, 0, 4) (1.0, {0: {1: 1.0}, 1: {2: 1.0, 3: 0}, 2: {4: 1.0}, 3: {4: 0}, 4: {}})
化合物0から化合物4までの経路は1.0であり、すなわち最下流の物質を生成するpathが存在することを意味する。
(複数存在する場合、もっとも短い経路で、かつ若いインデックスを返すことに注意)
ここで、枝分かれのうちの一つの"capacity"を0.0に設定する(遺伝子を欠損させる)。
>>> Graph[2][4]["capacity"] = 0.0 >>> nx.maximum_flow(Graph, 0, 4) (1.0, {0: {1: 1.0}, 1: {2: 0, 3: 1.0}, 2: {4: 0}, 3: {4: 1.0}, 4: {}})
もう片方の枝分かれが生きているので、化合物0から化合物4までの経路は依然として1.0である。
もう一つの枝分かれを欠損させると、
>>> Graph[3][4]["capacity"] = 0.0 >>> nx.maximum_flow(Graph, 0, 4) (0, {0: {1: 0}, 1: {2: 0, 3: 0}, 2: {4: 0}, 3: {4: 0}, 4: {}})
化合物0から化合物4までの経路は0になった。
実際にはKEGGなどからPathwayを取ってきて、目的のゲノムなどの機能組成を構築し、networkXに埋め込む必要があるがここでは割愛する。
cutadaptをlocal installする
アダプタープライマートリムツールであるcutadaptをローカルでインストールしたバージョンで動かしたいという場合の方法。
cutadaptはpipでインストールすることができるが、実行時に実際に直接叩いているのは以下のスクリプト。
#!/usr/bin/python import sys try: import _preamble except ImportError: pass from cutadapt.scripts import cutadapt cutadapt.main()
すなわち、ここに書かれている/usr/bin/pythonを使って以下のスクリプトを実行しているということ。
ローカルのバージョンで使いたいという場合は、ローカルにpythonをインストールし、そこからpipでcutadaptをインストールする。
以下は例。
$ cd $WORK ##local_installするPATHへ
$ wget https://www.python.org/ftp/python/2.7.14/Python-2.7.14.tgz
$ tar zxf Python-2.7.14.tgz
$ cd Python-2.7.14
$ wget https://bootstrap.pypa.io/get-pip.py
$ ./python get-pip.py
$ ./python -m pip install cutadapt==X.X.X ##versionを指定する
$ cat > ./cutadapt.py ##コピーして書いてしまう
#!/usr/bin/python
import sys
try:
import _preamble
except ImportError:
pass
from cutadapt.scripts import cutadapt
cutadapt.main()
$ ./python ./cutadapt.py -a AACCGGTT -o output.fastq input.fastq ##実行追記
最新のバージョンではスクリプトの中身が変わっている。
# -*- coding: utf-8 -*-
import re
import sys
from cutadapt.__main__ import main
if __name__ == '__main__':
sys.argv[0] = re.sub(r'(-script\.pyw?|\.exe)?$', '', sys.argv[0])
sys.exit(main())
MAFFTでアラインメントの距離行列を得る
MAFFTのTipsに
"How to get a distance matrix, instead of alignment"
の記述があるものの、まだドキュメントが作成されていないようだったので調査。
mafftのドキュメントにはないが、binaryの中に"mafft-distance"(source)と言うものがあったので使ってみた。
mafft-distance (aa) Version 7.222
$ mafft-distance protein.faa Calculating i-j scores ... 1-2 d=0.00 l=241,240 1-3 d=0.00 l=241,240 1-4 d=0.02 l=241,240 1-5 d=0.00 l=241,240 1-6 d=0.01 l=241,240
というような感じで出てくる。各カラムは左から、i番目とj番目の配列、距離スコア、アライメントの長さ(i,j)であろう。
小数点第2位までしか出してくれないのはそれでいいのかという思いと、Protein alignmentのScoreはどうなっているのかという疑問があるが、とりあえずこれをパースして距離行列を作ることはできそう。
ただし、このツール自体はアラインメントを取っていないので、アライメントを取った後の配列で計算したほうが良さそう。
("instead of alignment" ではなくて "after alignment" になってしまっているが)
$ mafft --globalpair --maxiterate 1000 protein.faa > protein.faa.mafft $ mafft-distance protein.faa.mafft
追記
ソースコードをよく見たら(ヘルプには出てないが)"-p"というオプションを付けると小数点6位まで(8.6fで)、かつ行列に変形しやすい形で出してくれた。
$ mafft-distance -p protein.faa.mafft > protein.faa.mafft.dist
Pythonで配列のJaccard距離を計算する(scipy.spatial.distance.jaccard)
Jaccard距離とは2配列間の距離(類似性の逆)をその要素の正誤によって求める指標である。
しかし、配列の要素がNaNかNaNでないか(または0か0より大きいか)を区別したい場合と、完全に値が一致しているかしていないかを区別したい場合などがある。
scipyにはscipy.spatial.distance.jaccardとして提供されているが、挙動がよくわからなかったのでテストした。
以下のような配列を用意する。
>>> from scipy.spatial.distance import jaccard >>> import numpy as np >>> array = np.array([[np.nan, np.nan],[np.nan, 0],[1, 2],[2, 2]]) >>> array array([[ nan, nan], [ nan, 0.], [ 1., 2.], [ 2., 2.]])
普通にJaccard距離を取ると
>>> jaccard(array[:, 0], array[:, 1]) 0.75
即ち[ 2., 2.]のみが一致したとみなされ、距離は3／4となっている。
NaNかNaNでないかの違いだけを見たい場合、boolに変換する。
>>> np.array(array, dtype=bool) array([[ True, True], [ True, False], [ True, True], [ True, True]], dtype=bool) >>> jaccard(np.array(array[:, 0], dtype=bool), np.array(array[:, 1], dtype=bool)) 0.25
[ nan, 0.]のみが不一致とみなされた。
ここで注目すべきは、非常に気持ち悪いがnp.nanがTrueになるという点である。
Jaccardでは[ False, False]の比較は無いものとして扱われるので、例えば
>>> np.nan_to_num(array) array([[ 0., 0.], [ 0., 0.], [ 1., 2.], [ 2., 2.]]) >>> jaccard(np.nan_to_num(array[:, 0]), np.nan_to_num(array[:, 1])) 0.5
とすると、[ 1., 2.], [ 2., 2.]の比較が行われ距離が１／２となる。
dtype=objectだとnp.nan_to_numが適用されないのでstrを含むデータを扱う場合は頑張って置換するしか無い。
NaNをFalse、0をTrueとして扱いたい場合もあるかも知れないが、これも簡単な方法はないので頑張って置換するしか無い。
むしろPandasのDataFrameにねじ込んで、isnull(), notnull(), all(), any(), sum()を駆使するほうが楽。
tryalnigro.hatenablog.com
同様の挙動をするものにscipy.spatial.distance.hammingやscipy.spatial.distance.kulsinskiがある。
hammingは[ False, False]の比較も行うタイプのもので、kulsinskiは分子分母にoffsetが設けられたもの。
