ERRORこの広告は、90日以上更新していないブログに表示しています。graphを利用したメタゲノム用の高速アセンブラである。MEGAHITはFASTQをインプットとし、Paired-endとSingle-endの両方に対応している。-2か-12で指定する。Single-endは-r。ミトコンドリアの完全ゲノムを得るために、ミトコンドリアにマップされたリードのみを抽出し、いくつかのアセンブラでテストをしていたときのこと。v1.1.1-2-g02102e1)公開データの中には、Paired-endとSingle-endの両方があったので、Paired-endのデータは-1,-2で、Single-endは-rで指定していた。すると、なぜかSingle-endのときにのみ完全ゲノム1本がきれいに出力された。GA-IIで、長さも100bpであり、決して他に秀でてうまくアセンブルできるとは思えなかった。そこで、Paired-endのものもForwardだけを用いて-rでアセンブルしてみたところ、(孤立した配列はいくつか混じっていたものの)完全長のコンティグが生成できた。MEGAHITはPairの情報をそこまで利用していないという認識があったので、この結果には驚いた。他のバージョンは試していないが、Singleにしたほうが断然良い結果であった。メタゲノムデータに対する場合はわからないが、Singleにしたほうが良い結果が出る場合が他にもあるかも知れないので注意したい。18:00MEGAHITのPaired-endアセンブリが怪しいTweet広告を非表示にするもっと読むコメントを書く2017-12-06Bowtie2で許容編集距離を設定するBowtie2のoptionは複雑で、とりあえずデフォルトで回してからSAMをパースしている人は多いかもしれない。マニュアルを良く読むとスコアをきちんと設定すれば色々設定できる。Manual特に以下の部分。ミスマッチペナルティ(最大,最小)、リードのクオリティでスコアがバラける。N(曖昧塩基)のペナルティ、Bowtie2では他の縮退塩基もNとして扱われる。参照配列のギャップペナルティ(開始,延長それぞれについて規定)これを編集距離へを変化させるには、それぞれを1に設定すれば良い(ギャップ開始は0)。length.よって、編集距離1以内のものだけを出力するオプションは以下のようになる。"C,-1,0"”C”については例示すら無いが一応これでちゃんと動くのは確認した。1を超えるものがなかった。04:04Bowtie2で許容編集距離を設定するTweet広告を非表示にするもっと読むコメントを書く2017-11-24PythonのListとDictionaryの検索機能の速さ"in"を使って要素の有無を判定する時、ListとDictionaryで速いのはどちらなのだろうと疑問に思ったので調べた。0m0.143s単純に検索するだけならList型のほうが速い。(何に使うかは知らない)ただ、普通は検索したあとそのindexに対してカウントしたり何かしらアクションを加えるので、結局Dictionaryの方が速い。要素の数が大きくなったときはDictionaryの方が早いようだった。objectはListに比べて作るのに時間がかかるため、本来ならforの回数を変化させて線形近似しなければいけないが、速さがDictionary>Listの例が取れたのでそこは割愛した。d"の検索速度に差はないようである。15:01PythonのListとDictionaryの検索機能の速さTweet広告を非表示にするもっと読むコメントを書く2017-11-13PythonのnetworkXを使ってpathway解析をするある遺伝子(機能)セットと代謝ネットワークデータを用いて、上流の化合物から下流の化合物に行くpathwayがあるか調べたい。PythonのnetworkXを使って最短経路を出すドキュメントはたくさんあったが、エッジの有無によってFlow解析についてはあまりなかったのでまとめた。Flow解析をするには、edgeに"capacity"を設定する必要がある。いう枝分かれを1つ持つ単純なネットワークを、以下のように設定する。1.0edgeの"capacity"を1.0に設定することは、その遺伝子が含まれることを意味する。{}})化合物0から化合物4までの経路は1.0であり、すなわち最下流の物質を生成するpathが存在することを意味する。(複数存在する場合、もっとも短い経路で、かつ若いインデックスを返すことに注意)ここで、枝分かれのうちの一つの"capacity"を0.0に設定する(遺伝子を欠損させる)。{}})もう片方の枝分かれが生きているので、化合物0から化合物4までの経路は依然として1.0である。{}})化合物0から化合物4までの経路は0になった。実際にはKEGGなどからPathwayを取ってきて、目的のゲノムなどの機能組成を構築し、networkXに埋め込む必要があるがここでは割愛する。installするアダプタープライマートリムツールであるcutadaptをローカルでインストールしたバージョンで動かしたいという場合の方法。cutadaptはpipでインストールすることができるが、実行時に実際に直接叩いているのは以下のスクリプト。cutadaptcutadapt.main()すなわち、ここに書かれている/usr/bin/pythonを使って以下のスクリプトを実行しているということ。ローカルのバージョンで使いたいという場合は、ローカルにpythonをインストールし、そこからpipでcutadaptをインストールする。最新のバージョンではスクリプトの中身が変わっている。alignment"の記述があるものの、まだドキュメントが作成されていないようだったので調査。mafftのドキュメントにはないが、binaryの中に"mafft-distance"(source)と言うものがあったので使ってみた。l=241,240というような感じで出てくる。各カラムは左から、i番目とj番目の配列、距離スコア、アライメントの長さ(i,j)であろう。alignmentのScoreはどうなっているのかという疑問があるが、とりあえずこれをパースして距離行列を作ることはできそう。このツール自体はアラインメントを取っていないので、アライメントを取った後の配列で計算したほうが良さそう。protein.faa.mafft追記ソースコードをよく見たら(ヘルプには出てないが)"-p"というオプションを付けると小数点6位まで(8.6fで)、かつ行列に変形しやすい形で出してくれた。17:14MAFFTでアラインメントの距離行列を得るTweet広告を非表示にするもっと読むコメントを書く2017-10-23Pythonで配列のJaccard距離を計算する(scipy.spatial.distance.jaccard)Jaccard距離とは2配列間の距離(類似性の逆)をその要素の正誤によって求める指標である。配列の要素がNaNかNaNでないか(または0か0より大きいか)を区別したい場合と、完全に値が一致しているかしていないかを区別したい場合などがある。scipyにはscipy.spatial.distance.jaccardとして提供されているが、挙動がよくわからなかったのでテストした。以下のような配列を用意する。2.]のみが一致したとみなされ、距離は3/4となっている。NaNかNaNでないかの違いだけを見たい場合、boolに変換する。0.]のみが不一致とみなされた。ここで注目すべきは、非常に気持ち悪いがnp.nanがTrueになるという点である。2.]の比較が行われ距離が1/2となる。dtype=objectだとnp.nan_to_numが適用されないのでstrを含むデータを扱う場合は頑張って置換するしか無い。NaNをFalse、0をTrueとして扱いたい場合もあるかも知れないが、これも簡単な方法はないので頑張って置換するしか無い。sum()を駆使するほうが楽。tryalnigro.hatenablog.com同様の挙動をするものにscipy.spatial.distance.hammingやscipy.spatial.distance.kulsinskiがある。False]の比較も行うタイプのもので、kulsinskiは分子分母にoffsetが設けられたもの。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。