ryought / dbgphmm Goto Github PK
View Code? Open in Web Editor NEWProbabilistic genome assembler based on de Bruijn Graph and profile HMM
Home Page: https://ryought.github.io/dbgphmm/dbgphmm/
Probabilistic genome assembler based on de Bruijn Graph and profile HMM
Home Page: https://ryought.github.io/dbgphmm/dbgphmm/
add_seq
とかでkの値を指定しないといけないようになっている
#28 で言われていることのもう一度整理
stat
read-stat
compare
dbg
cdbgをcompressする
→
基本的には足りないk-merを足していくような行為
optimalに近づくにつれて良くなるはず
trueから、全basisを上下に動かした時のスコアをチェックするみたいなのも面白い。
(摂動を加えてスコアがどう変化するか?)
単体でcircularなものができた時、端が同じkmerが有ればinjectできる
wikiにまとめる?
ref.faを作る
ランダムに
cmd generate > ref.fa
リファレンスから切り出す
seqkit faidx ...
ref.faからリードを生成する
cmd sample ref.fa > reads.fa
dbgを作った時の統計を作る
cmd stat -k 8 ref.fa
cmd stat -k 8 reads.fa
dotを作る
dbgを作って、readの出力確率を計算する
cmd forward dbg.fa reads.fa > reads.tsv
readからdbgを作って、一番良いコピー数を推定する
(答えを与えると、正解のコピー数と比較できる)
cmd optimize reads.fa --answer ref.fa > dbgs.tsv
正解のkmer、エラーkmerの割合
-strict
ピッタリ長さLのリードを出力する
-end-prob
ending probabilityに従って止める
ひょっとしたらend-prob=0にすれば同じことかもしれない
2つのコピー数の間の距離
(ベンチマーク等で、正解のコピー数にどの程度近いか?みたいなことを計測したい)
dbg間の距離
(テストで、serializeして元に戻るか?などを確かめたい)
とりあえず実装してためす。
ボトルネックを特定して、高速化する
今はまずナイーブに実装する段階。
ランダムの場合は何個入れるかも選べるようにする
clapのstructを渡した方がミスが少ない
ひょっとしたら
各cycle edgeを検出してその個数だけcycleを足し引きする
だけで良いかもしれない
実験スケジュール全体の概観。
kが固定されれば、faから作れるdbgは一意
dbgから作る時は、euler circuitを適当に選んで作る
色々なコピー数で、リードの出力確率を計算できるようにする
1/2のk-merコピー数推定は動いていそう。
・高速化
Forward計算にかかる時間(readの塩基、kmerの個数あたり)
MMWC計算にかかる時間
・精度
今のところ失敗する例が見つけられない…
・他手法との比較
比較できる手法があるか?dbgベースの
完了
emitableでないノードもグラフ中には許すが、
transition probabilityは0になるようにする。
dbgphmmを作っていく意味があることを示したい
7/26まであと2週間
文字の複雑さ
k-mer dbgのkを変えた時の性質で測れるのでは
逆にそれをinputにすると、似たような複雑さ(=アセンブリの難しさ)のゲノムが生成できるのでは?
経緯:EMアルゴリズムがメインの推定手段になりそう。その際のForward/Backwardアルゴリズムはn
=(dBG中の全ノード数) m
=(リードの総塩基数)としてnm
の表を埋める必要があるので重い。
4^k
。しかしだんだん後になると(k=40等)ほとんどがユニークなk-merになり、またreadのエラーも除外されていると想定できるので、ノード数はゲノムサイズに近づく。速度を測定できるようにする
(配列をランダムに作る。)
PHMMのモジュールだけ考えれば良い
理論的な保証
EMアルゴリズムの正当性
Q関数が一個前よりも下がっているか Q(H_t,H_t-1) > Q(H_t-1,H_t-1)
を確かめる
実験的な保証
リファレンスゲノムを分断する
・reference, readのfastaを読めるようにする
・実験結果もcsv等で(parseできる形で)出力するようにする
read.fa, dbg.faを与える
read.faからdbgを作って、周辺をうろつく
リファレンス生成から自動でやる
シードを適当に振って、リファレンスを生成して、周辺をうろついた結果を報告する
dbgのコピー数が最適化されていく様子を動画にできるようにする
余力があれば3dでinteractiveにすることも可能
少なくともうまく最適化が回るかを確認してからで良い…
(デバッグに便利という説もあるが)
optimalなD_k
からD_k+1
を計算する方法
k+1に変換したときにどこがambiguousか?
2以上in-2以上outのノード(intersection)
P(D_k+1)
を試して一番スコアが高いものを採用する
C
:choice数, I
:intersection数としてC^I
になるk+1
に対してもk
と同じEMアルゴリズムを動かす
k
がかなり大きくなって、intersection同士が離れ始めたら有効そうk→k+1に行けるようにする
まずはとりあえず例によって小さいデータで実験して、うまく行きそうかをチェックする
read中の出現回数と比較
(この情報があると最適化が楽になるのでは?)
各リードについて、
read_id, prob
を出力する
最後にコメントで合計の確率を書いておく
HPCの考慮ありのprofile hmmモデルを作るのもあり
ループの個数がどういう感じで分布するのだろうか…?
そして平均ループサイズは…?
kのサイズ
8~40
ゲノムサイズ
error profile
端っこ(NNNNN)を入れる入れない
7/20の調査
y=(freqs/mean_coverage) x=(copy_nums) としたときに ||x-y||^2
を最小化したい。以前の定式化では cx-y
を使っていたが、割り算すれば良いことに気づいた。カバレッジを引数に与えなくて良いので若干キレイ
一般のinteger least square(ILS) min(||Ax-y||^2, x)
はNP-hard。xをサイクル基底上の係数ベクトルとすると、Aはサイクル上に出てくるk-merに+1/-1が立っている行列で、ILSとして解ける。
逆はAが整数行列だったら成り立つ。Aが整数行列・binary行列の場合/x,yが非負の場合は特別早く解ける、みたいなことがありえるか?
ILSはGPUの位置推定で使われるらしく、Aを上三角行列にして、力任せで探索する方針が取られるらしい。
meta heuristics的な方法。
SAStateをcdbg/phmmむけに作る
optimize後の結果をfastaで出力する
変化させるパラメーターは
sigmaからかなり離れると桁落ちする
分散で良い説もある。
log(Forward) + lambda log(Prior)
大きなスケール
小さなスケール
・dbg上でやるのは比較的簡単
全nodeを列挙して、各nodeに対して、in-edges-sum = out-edges-sumかを判定すれば良い
・cdbg上でできるようにしたい
kは答えとして与える
k+1を早く推定できるか?
・答えとして与えたkから、可能なk+1を全部残したk+1のdbgを作る
・その上で最適化する
大きいstructにまとめて、serdeでjson出力する
indentつければ、人間にもみやすいはず
複雑にしようと思えばいくらでも出来て、沼にハマってしまう。
細かくできる作業に分割してどんどんやっていこう
直近で出来ていないといけないこと
10*4=40
ステート考える)、またその上位10ステートのみを保存して次にいくという感じ
すぐにやらないといけない感じではないが、いずれやるとメンテナンス性が上がりそうなもの
u128
等で管理するdbgの実装を入れたい。smerとしてファイルは作った2回目以降のforwardを高速化する
変化してないところはそのままにする
どのベースが出てきたかも覚えておく
1Mbpくらいの領域でも動くようにしたい
kmerの判定力をチェックする
リファレンスを切り刻む
セントロメア等でためす
どのくらいエラーkmerが減るのかをテストする
長さに対する計算時間を測定する
k k+1をテストする
dbgの実装を変更する
・可変k、small k、large kでも動くようにする
可変k(simple pathを潰す。nodeがkmer、edgeがk-1merの表記の時はnodeが>kmerになる場合があるということ。)
small k(k-merをuint64等で持っておく)
large k(各ノードのk-merを持っておくのはもったいない。ただし分岐点のところだけわかれば良い…?)
Forwardの軽量化
・0-copy k-merは落とす(cdbg→phmmに変換する時に)
・sparsevecを使うようにする(sparse vecを使って良いかどうかの検証)
大文字の変数とかがおこられている
あとむかし作ったいらないファイルは消す
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.