add_seqとかでkの値を指定しないといけないようになっている

#28 で言われていることのもう一度整理
stat
read-stat
compare

• 正しいk-merの存在比(kが大きくなるとだんだん少なくなるはず)
• 正しいコピー数との差

dbg

• k=32以下ではkmerをusizeとして表す
• unitigなノードを端折る(隣に行くときに1文字しか追加できない制約をやめる)
• k可変を許す

cdbgをcompressする
→

基本的には足りないk-merを足していくような行為
optimalに近づくにつれて良くなるはず

trueから、全basisを上下に動かした時のスコアをチェックするみたいなのも面白い。
(摂動を加えてスコアがどう変化するか？)

単体でcircularなものができた時、端が同じkmerが有ればinjectできる

wikiにまとめる？

option(params)

• dbg
  • k
• phmm
  • p_mis, p_open, p_ext等
• prior
  • ave
  • std
• sa
  • init_temp
  • temp_deg

commands

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して元に戻るか？などを確かめたい)

とりあえず実装してためす。
ボトルネックを特定して、高速化する
今はまずナイーブに実装する段階。

forward algorithm

dbg construction

• zero
• true
• read
• random

ランダムの場合は何個入れるかも選べるようにする

• hmmのpathの情報を使う(やはりBaumWelch的に？)
• 全体にではなくk-merに温度を割り振る　まだ変化しうるk-mer、みたいに

clapのstructを渡した方がミスが少ない

ひょっとしたら
各cycle edgeを検出してその個数だけcycleを足し引きする
だけで良いかもしれない

実験スケジュール全体の概観。

A 計算する意味があるか

• A1 真のコピー数が確率最大になるか？(エラー込みのリードを作る。事前分布を適当に設定。kやゲノムを変えて、真のコピー数の周りを徘徊して、最大かどうかを検討する)
  • 問題になる可能性がありそうな場所：事前分布の形状、kが大きい時の状況、リードの端っこの扱い方
• A2 simulated annealingで問題ないか？近傍の設定方法はcycle basisで問題ないか？
  • 問題になる可能性がありそうな場所：cycleのサイズが大きすぎるとbasisだけだと近傍が遠すぎる

B 計算を高速化できるか、コードがもう少し簡潔になるか

• B1 Forward algorithmm
  • differential update コピー数が変わったところだけ変更する
  • chopped reads kのdbgのforwardを動かすのに、大体2kぐらいの長さしか関係なさそう。
• B2 cdbg/dbg construction
  • readの端っこの処理は今のままで問題ないか？(明らかにここで終わるわけがない、みたいな場所は残しておく？)
• B3 kからk+1が高速に作れるか？
  • この場合にも、dbgは作れなくてもcdbgは作れる、などがありそう。

C ゲノムスケールで動くか？

• dbg作成やサイクル検出に必要な時間、メモリの見積もり。KIR領域程度のサイズからリードを生成して、dbgを作って、どのぐらいのサイズになるか検証(雑には今のリードセットからdbgを作るのでもよい) kに応じて、kmerの個数、次数分布、サイクル長分布
• forwardを高速に回せるようになったか？

kが固定されれば、faから作れるdbgは一意
dbgから作る時は、euler circuitを適当に選んで作る

• optimize、コマンドとして実行できるようにする
  • 答えありでコピー数周辺を彷徨く
• copy数が0があるときに、hmmの中のinit/trans-probがバグったりしないか？
  • init_prob/trans_probを表示させてみて変な値がないか検証
  • dbg_faから作ったcdbgと、read_faから作ってcopy_numをdbg_faから作ったcdbgで確率が同じになることを検証
• 答えのcopynumをcdbgから計算する処理を作る
  • true k-merが存在しないことが結構ある？
  • リードの長さ、kを変更してtrue k-merがある状況を考える
• kmerの統計値を充実させる
  • リードとの相関: 答えがある場合、正しい/誤りkmerの割合
  • リファレンスの複雑性: kを変えて、何merでユニークになるかを検証できる
• annealerのlogを作る(historyを出力するように)
  • annealerのlogの種類とは？各サイクル基底ごとの選択率
• サブコマンドdotを作る(statは標準出力に書き出すようにする)
• コンパイル時間を早める(並列？自動で裏でビルドする？)
• sampler
  • リード長を全部同じにするモード、全長吐き出すモードも作りたい
  • 終わった後に統計値(エラー率、ノードカバー率)を出力したい

色々なコピー数で、リードの出力確率を計算できるようにする

1/2のk-merコピー数推定は動いていそう。

・高速化
Forward計算にかかる時間(readの塩基、kmerの個数あたり)
MMWC計算にかかる時間
・精度
今のところ失敗する例が見つけられない…
・他手法との比較
比較できる手法があるか？dbgベースの

完了

• cdbgとdbgで同じ値が出るか
• cdbgのemissionがNのノードを無視する

emitableでないノードもグラフ中には許すが、
transition probabilityは0になるようにする。

dbgphmmを作っていく意味があることを示したい

• 真のコピー数がスコア最大になること
  • リードから作る
  • 色々なコピー数をためす。
• cycle basisのsimulated annealingで、スコアを最大化できること
  • 何割ぐらいが訪問された？

7/26まであと2週間

• 最適化の方法を色々試す
  • 勾配法
  • k-mer温度、基底温度
• 最適化が難しい理由を探る
  • 地形の可視化。(どうやって？)
  • ランダムなシードから勾配法を登る。何割ぐらいの確率で正解に行き着くか？
• モデルとしての良さを示す図を作る
  • 正解をシードにして、勾配法を登る。

文字の複雑さ
k-mer dbgのkを変えた時の性質で測れるのでは

逆にそれをinputにすると、似たような複雑さ(=アセンブリの難しさ)のゲノムが生成できるのでは？

経緯：EMアルゴリズムがメインの推定手段になりそう。その際のForward/Backwardアルゴリズムはn=(dBG中の全ノード数) m=(リードの総塩基数)としてnmの表を埋める必要があるので重い。

考える方向性

• sum→maxの場合はviterbiになり、基本的にはアライメントのDPと同じ。アライメントの高速化手法が参考になる。
• アライメントと違う点(つまりこの問題固有の性質)としては
  • 各ステートをemissionに使った回数f=(forwardとbackwardの足し算?)のみが必要
  • 何回も似たモデルに対して実行する。
  • モデルのパラメーターも、そこまで急激に変化するわけではない

高速化のアイデア

• nを少なくする
  • とりあえずコピー数0のノードは除外できる
  • 最初(例えばk=8等)は全k-merが出現しそうで、ノード数は4^k。しかしだんだん後になると(k=40等)ほとんどがユニークなk-merになり、またreadのエラーも除外されていると想定できるので、ノード数はゲノムサイズに近づく。
  • しかしノード数がゲノムサイズになったとて、毎回リードをリファレンスにmapしているようなもので、かなり重い。もちろん領域は限定するが。
• mを少なくする
  • readをkmer単位に分解する
• 埋め方を工夫する
  • baseが前回どこから出力されたかを覚えておく。その近辺だけで動かす。

速度を測定できるようにする
(配列をランダムに作る。)
PHMMのモジュールだけ考えれば良い

そのほか

• ちょっとした高速化テク
  • 2次元行列を1次元行列にする
  • メモリ管理をもうちょっと頑張る(コピーしなくて済むところなど)
  • logsumexp等、何度も動かしている部分を見直す
• またはsum特有の方向性があるかもしれない
• その高速化方法がどの程度良いかをどう測るか？
  • exactな方法との相関を見る

• dbg.fa(正解) reads.fa(リード)を受け取る
• reads.faからdbgを作る→statと同じ内容を出力
• それに加えてdbg.faの正解があるから、以下の補足情報を作る
  • 正しいkmerの数
  • 正しいkmerのdbg
  • 正しいゲノムサイズ
• optimizeモード
  • 正解のcopy数から初めて、SAでランダムに遷移する。各ステートでの確率を出力する
  • priorのゲノムサイズは真の値を使う(分散は適当に指定する)。

理論的な保証
EMアルゴリズムの正当性
Q関数が一個前よりも下がっているか Q(H_t,H_t-1) > Q(H_t-1,H_t-1) を確かめる

実験的な保証
リファレンスゲノムを分断する

・reference, readのfastaを読めるようにする
・実験結果もcsv等で(parseできる形で)出力するようにする

(1)

read.fa, dbg.faを与える
read.faからdbgを作って、周辺をうろつく

(2)

リファレンス生成から自動でやる
シードを適当に振って、リファレンスを生成して、周辺をうろついた結果を報告する

dbgのコピー数が最適化されていく様子を動画にできるようにする

• rustでdbgの点の位置計算をする。2dで
• annealerのcopy数のhistoryで時系列の画像を出力する

余力があれば3dでinteractiveにすることも可能

少なくともうまく最適化が回るかを確認してからで良い…
(デバッグに便利という説もあるが)

optimalなD_kからD_k+1を計算する方法

k+1に変換したときにどこがambiguousか？
2以上in-2以上outのノード(intersection)

方針

• 色々なP(D_k+1)を試して一番スコアが高いものを採用する
  • 全choiceで動かす: C:choice数, I:intersection数としてC^Iになる
  • 各choiceで確率最大のものに割り当てる
• k+1に対してもkと同じEMアルゴリズムを動かす
  • ゲノムサイズを変える必要はない…
  • k+1として可能なk-merのみを入れてk+1を作る。コピー数の初期値は、ambiguousなところは等分する。
  • リードk-merを全部扱わなくて良くなる反面、結局最初からk+1で実行しているのと同じでは？
• intersectionにのみ注目する
  • 多分kがかなり大きくなって、intersection同士が離れ始めたら有効そう
  • 実験的に確かめてみたら？

実装

k→k+1に行けるようにする

スケジュール

まずはとりあえず例によって小さいデータで実験して、うまく行きそうかをチェックする

read中の出現回数と比較
(この情報があると最適化が楽になるのでは？)

各リードについて、
read_id, prob
を出力する
最後にコメントで合計の確率を書いておく

• LASTのパラメタファイルが使える？
• 自分で作る場合、pair-hmmを作ることになる？

HPCの考慮ありのprofile hmmモデルを作るのもあり

• 探索空間(PCAして、軌跡が線で表示されているようなやつ)
• 停留点のアセンブリの状況(fa)

ループの個数がどういう感じで分布するのだろうか…？
そして平均ループサイズは…？

kのサイズ
8~40

ゲノムサイズ

• 100bp
• 1000bp(1kbp)
• over 20kbp

error profile

• no error
• 1%(hifi)

端っこ(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的な方法。

  • simulated-annealingで選択にスコアによってbiasを入れたもの。
  • 基底を増やして、色々なサイクルを考慮するもの。
  • 各ステップでMST的な方法を使って、一番必要なサイクルを追加するもの。
  • STをたくさん計算する。STを固定すれば、freqs→copy_nums (in f)は一意で、四捨五入でuにできる。そのコンセンサスを取る？

• optimalなスコアからためす(どこにも行かない？)
  • →違うところに行ってしまう
• forwardの部分を並列化する
  • →Forwardサブコマンドでまず実験
• 全体のどのあたりを探索しているのか？
  • 提示された状態を列挙する
• 可能性のある方向にだけ移動するようにする

• 全部2倍,3倍にしてもスコアは同じか？
• 勾配はどんな感じになっているか？

SAStateをcdbg/phmmむけに作る

optimize後の結果をfastaで出力する

• パラメーターを変化させた時の
  • annealer・grad (途中で止まる)
  • freq-em、full-em constant (kが大きいとおかしい)
  • full-em linearGradient (いい感じ)
• kが大きいと真のk-merのfreqが小さくなる。ヒストグラムみたいな感じで出力する？

変化させるパラメーターは

• データセット
  • length
  • error_rate
  • depth
• 最適化
  • k
  • method

sigmaからかなり離れると桁落ちする

分散で良い説もある。

log(Forward) + lambda log(Prior)

大きなスケール

• 計算する意味があることを示す図を完成させる
• 小さいゲノムについて最適化を走らせる。1000bpのゲノムに対して、真のコピー数が見つかるか確認

小さなスケール

• Forward, Optimizeコマンドを完成させる
• 適当にリファレンスから生成して、1塩基変えて、得点が下がるかチェック
• Optimizeで、正解のデータがある時の統計値の出力を考える
• SimulatedAnnealingで正解のコピー数の周りをフラつくのを事前分布のみで(リードの出力確率を考えずに)やってみる
• (ここまでできたら、真のゲノムの確率が最大になることを示せる？)
• 温度戦略等を変えて、実際に最適化を走らせる

・dbg上でやるのは比較的簡単
全nodeを列挙して、各nodeに対して、in-edges-sum = out-edges-sumかを判定すれば良い

・cdbg上でできるようにしたい

kは答えとして与える
k+1を早く推定できるか？

・答えとして与えたkから、可能なk+1を全部残したk+1のdbgを作る
・その上で最適化する

大きいstructにまとめて、serdeでjson出力する
indentつければ、人間にもみやすいはず

複雑にしようと思えばいくらでも出来て、沼にハマってしまう。
細かくできる作業に分割してどんどんやっていこう

優先度高め

直近で出来ていないといけないこと

• Forward/Backwardを高速化する: 目標としてはk=32等で1Mぐらいの領域のリードサブセットに対して動くようにしたい。
  • 何列ぐらいを平均で使っているのか？
    • 10/27,28 可視化したい。一番スコアの高いものが対角線に来るようにソートしてあげる。どのぐらいの順位になるのだろう？
  • phmmをsparsevecを使って書き直す。各塩基に対して、スコア(確率)が高かった上位10ステートのみを保存しておく。その次の塩基は、その子供に対しては全部考えて(なので10*4=40ステート考える)、またその上位10ステートのみを保存して次にいくという感じ
    • 10/27 layerをtraitにすれば良い？
    • 10/28 deletion stateもひとつ前の上位陣だけを切り取ればよさそう(連続でdeletionした場合なので)。
  • 0-copyのkmerを落とすようにする
  • 2度目、スコアが変わっていないところは修正しなくて良い。(分割してもスコアが変わらないような気がするのは何故だろう？)
• 大きなkでも動くようにする
  • 全node/edgeにkmerを持っておく必要はなくて、親をたどっていけば自明に定まる
  • nodeをcontractできるようにする
• k→k+1のアルゴリズムを作る

すぐにやらないといけない感じではないが、いずれやるとメンテナンス性が上がりそうなもの

• u128等で管理するdbgの実装を入れたい。smerとしてファイルは作った
• graph, deBruijnGraph等、データ構造と手法を分離したい
  • algorithmから、具体的にstructに依存する部分を分離して、trait(インターフェース)を介してアクセスするようにする

下書き

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を使って良いかどうかの検証)

大文字の変数とかがおこられている
あとむかし作ったいらないファイルは消す