ここまでのあらすじ
非モデル生物のゲノムをショットガンシーケンスで読み、リードのフィルタリングを行いました。あとはアセンブリをすればドラフトゲノム完成！な訳でしたが…
アセンブリが終わらなーーーーーい！
まず初めに試みたのは有名なde novo アセンブラのPlatanus（Kajitani et al.,2014）。これは特にパラメータ設定とか必要なくて楽ちんなアセンブラです。
ショットガンシーケンスの強み「短い配列をたくさん読む」を活かすためにk-mer をデフォルトより短く設定しました。
結果、計算が止まってしまいました。
計算にはスパコンを使っているのでメモリ不足などではなく、アセンブラの限界だったと見えます。KmerGenie (Chikhi and Medvedev, 2015) という最適k-mer を求めてくれるソフトを使ってみると最適k-mer は101 と出たので、試しにそれも試してみましたが、短すぎて全然使い物にならなかったのでした。
これさえあればオッケー、と思っていたPlatanus がダメぽだったので電子の海をさまよって色々試しています。以下、日記です（現在進行形）。
ちなみに、ドラフトゲノムの使い道はRAD-seq のリファレンスです。
ゲノムアセンブリの記録
202018までの進捗
Platanus→k=32 以下では収束せず。Kmergenieを用いて出したk=101 は大きすぎて捨てる部分が多くなってしまうため不可。
Minimap2 で近縁種ゲノムにフィルタリング後データをマッピングしたものをhypo でポリッシング→RAD-seq には使えなかった
MaSuRCA (Zimin et al., 2013) →MaSuRCA アセンブラ - macでインフォマティクス (hatenablog.com)
名古屋大学の上坂一馬先生のサイトで発見。日本語でバイオインフォマティクスのソフトの説明が書いてあるサイトはめずらしいので非常に助かる。
Masuruca →アセンブリ自体はリファレンスなしでできるものの、インサートサイズとインサートサイズのSDが必要であるため、結局ゲノムへのマッピングが必要。
※これは後日勘違いであると発覚
インストールはanaconda3 を使った。スパコンではいつもの[conda activate 環境] が何故か効かなくて困ったが、[source activate 環境] で仮想環境移動ができた。よかった。
source activate masurca
masurca config.txt -o assemble.sh
nohup ./assemble.sh
20201219
学校に来たらアセンブリが止まっていた
エラーログがいくつか出ていた
エラーその１ WARNING! CA_PARAMETERS = cgwErrorRate=0.25 and LIMIT_JUMP_COVERAGE = 60 in config file should only be used for bacterial genomes; set cgwErrorRate=0.15 and LIMIT_JUMP_COVERAGE=300 for eukaryotes and plants!
パラメータがいくつかバクテリアのものになっていたらしい（めんどうくさがってHP のコードをマルパクした影響）
→対策：config.txt の中身を指示通りに書き直した。
エラーその２
ERROR: the runCA option ovlMemory is not recognized.
This runCA functionality changed in CABOG version 7.
Here are suggestions for new runCA option values:
ovlHashBits=23
ovlHashBlockLength= 30000000
(replaces ovlMemory=2GB)
ovlHashBits=24
ovlHashBlockLength=110000000
(replaces ovlMemory=4GB)
ovlHashBits=25
ovlHashBlockLength=180000000
(replaces ovlMemory=8GB)
ギャップクローズ段階で出るエラーのよう。ovlHashBitが何だかわからないし、公式サイトで検索しても引っかからず困ったが、とりあえず、エラー表示の中で勧められていたovlHashBits=23 を試すことにした。ovlHashBlockLength= 30000000 はconfig.txt の中に見当たらなかったし、何だかもわからなかったので一先ずスルーした。
20201220 またアセンブリが止まっていた
今回もギャップクローズ段階で問題が生じたらしい
ovlHashBits=24
ovlHashBlockLength=110000000
(replaces ovlMemory=4GB)
を試すことにした。やけっぱちである。
21：29 またアセンブリが止まっていた。
今、mazurcaではk=85 （オート）でやっているようだが、25~127 までサポートされているらしいので最小に設定してみようかな。というかk=14でやったらどうなるんだろう。
やってみた。
masurca config.txt -o assemble.sh
Error line 22 of configuration file 'config.txt':
bad value for GRAPH_KMER_SIZE, enter auto or number >= 15 and <= 151
15からいけるって！とりま15にしよ！
あと、CA_PARAMETERS = cgwErrorRate=0.15 に戻した。
20201221
また止まった。
いままでエラーメッセージが書かれていたrunCA0.outファイルの中を見てみる。
usage: runCA -d <dir> -p <prefix> [options] <frg> ...
-d <dir> Use <dir> as the working directory. Required
-p <prefix> Use <prefix> as the output prefix. Required
-s <specFile> Read options from the specifications file <specfile>.
<specfile> can also be one of the following key words:
[no]OBT - run with[out] OBT
noVec - run with OBT but without Vector
-version Version information
-help This information
-options Describe specFile options, and show default values
<frg> CA formatted fragment file
Complete documentation at http://wgs-assembler.sourceforge.net/
File not found or invalid command line option 'superReadSequences_shr.frg'
File not found or invalid command line option 'pe.linking.frg'
今までに見たことがないパターン。
superReadSequences_shr.frg
pe.linking.frg
の二つがないぞと怒られる。いままでのトライで出ていたのかわからないのが痛い。MaSuRCAは試行のたびにものすごく重いファイルをたくさん吐くので毎回消してしまっていた。一先ず今回のフォルダを確認する。
たしかにないけど、superReadSequences_shr.frg.tmpというかなり惜しい名前のファイルがある。中を見ると、途中で処理が止まってしまっていた。
原因不明。今度はnohup.out ファイルの中身を見てみる。
[2020年 12月 20日 日曜日 21:47:23 JST] Processing pe library reads
[2020年 12月 20日 日曜日 21:55:48 JST] Average PE read length 127
[2020年 12月 20日 日曜日 21:55:49 JST] MIN_Q_CHAR: 33
[2020年 12月 20日 日曜日 21:55:49 JST] Creating mer database for Quorum
[2020年 12月 20日 日曜日 22:05:44 JST] Error correct PE
[2020年 12月 20日 日曜日 23:03:08 JST] Estimating genome size
[2020年 12月 20日 日曜日 23:19:55 JST] Estimated genome size: 747917875
[2020年 12月 20日 日曜日 23:19:55 JST] Creating k-unitigs with k=15
[2020年 12月 21日 月曜日 00:46:39 JST] Computing super reads from PE
Unable to flush stdout: ディスク使用量制限を超過しました[2020年 12月 21日 月曜日 16:33:06 JST] Using linking mates
Unable to flush stdout: ディスク使用量制限を超過しました[2020年 12月 21日 月曜日 16:39:10 JST] Celera Assembler
gkStore::gkStore()-- GateKeeper Store 'CA/genome.gkpStore' doesn't exist.
Search pattern not terminated at -e line 1.
Search pattern not terminated at -e line 1.
[2020年 12月 21日 月曜日 16:39:12 JST] Overlap/unitig failed, check output under CA/ and runCA1.out
普通にディスク使用量制限を超過したからのような気がしてきた。
えー・・・これどうしたらええの。
とりあえず、ディスクの中を泣く泣く整理した。消していいかわからないけどデカすぎて自分のPCに移せないものがちらほら（涙）。
あと、config txtから以下の記述を削除した。ジャンピングライブラリーってメイトペアシーケンスに使われるものっぽい（？？？）ので。一応今回外注先が使ったライブラリ作製のプロトコルも見てみたけどそれらしき用語は出ていなかった。
#this parameter is useful if you have too many jumping library mates. Typically set it to 60 for bacteria and something large (300) for mammals
LIMIT_JUMP_COVERAGE = 300
masurca/README.md at master · alekseyzimin/masurca · GitHub
README を見直して無駄なパラメータがほかにもないか再度確認。
#whether to attempt to close gaps in scaffolds with Illumina data
CLOSE_GAPS=1
あとはこれを入れてみました
それから、CAのHPを見てようやくovlHashBlockLengthが何だかわかりました。
RunCA - wgs-assembler (sourceforge.net)
ovlHashBits=integer (デフォルト 22) (CA 7.0 の新機能)
（直訳）ハッシュテーブルのサイズをビット単位で指定します。これは固定サイズの割り当てであり、ovlHashBlockLengthやovlRefBlockSizeに基づいて変更されることはありません。
ovlHashBlockLength
（直訳）ハッシュテーブルにロードするシーケンスの量を塩基単位で指定します。デフォルトの100,000,000,000塩基をロードすると、 ovlHashBitsが使用するメモリに加えて1GBのメモリを消費します。
ovlRefBlockSize
（直訳）重複するジョブの数と各ジョブの実行時間を直接制御します。値を小さくすると、より多くのジョブが実行され、それぞれのジョブが終了するまでの時間が短くなります。この値が小さすぎるとオーバーヘッドが総時間の大半を占め、大きすぎると同時実行性が低下します。
これらのパラメータの選択は、オーバーラップを計算しているハードウェアに基づいて行われます。
ジョブ数やジョブが消費するメモしサイズを決めるタグだったようです。
ovlHashBits=23
ovlHashBlockLength= 30000000
(replaces ovlMemory=2GB)
大きいジョブにするとまた容量超過しそうでこわいので小さめに入れてみる。
今回メモリオーバーしたのはこのタグを消したせいもあるのかもしれない…。
という訳で、再びアセンブリスタート。
20201222
またもやアセンブリが止まる。今までで最速（スタートして３時間）。
nohup.out の中を見る。
[2020年 12月 21日 月曜日 23:09:23 JST] Processing pe library reads
[2020年 12月 21日 月曜日 23:17:17 JST] Average PE read length 127
[2020年 12月 21日 月曜日 23:17:17 JST] MIN_Q_CHAR: 33
[2020年 12月 21日 月曜日 23:17:17 JST] Creating mer database for Quorum
[2020年 12月 21日 月曜日 23:26:35 JST] Error correct PE
[2020年 12月 22日 火曜日 00:20:01 JST] Estimating genome size
[2020年 12月 22日 火曜日 00:35:45 JST] Estimated genome size: 747917875
[2020年 12月 22日 火曜日 00:35:45 JST] Creating k-unitigs with k=15
[2020年 12月 22日 火曜日 01:52:01 JST] Computing super reads from PE
[2020年 12月 22日 火曜日 01:52:01 JST] Super reads failed, check super1.err and files in ./work1/
super1.err and files in ./work1/ と言われたが、どっちもないんですけど！
という訳で原因はわかりません。
今までの試行と変えた部分は
LIMIT_JUMP_COVERAGE = 300
これ取って
CLOSE_GAPS=1
これいれたことくらいしかない。
ギャプクローズ段階より前に止まったから心当たりがあるとするなら前者かなあ…
ひとまず前者を入れて再試行。
16:35 またアセンブリが止まる。
[2020年 12月 22日 火曜日 15:48:38 JST] Super reads failed, check super1.err and files in ./work1/
またお前か…。
今回はwork1フォルダもsuper1.errもあるので中を見る。
super1.errの中身
terminate called after throwing an instance of 'std::runtime_error'
what(): Error allocating memory
child died with signal 6, with coredump
mv work1/overlap.overlaps work1/createKUnitigMaxOverlaps.Failed
mv work1/overlap.coords work1/createKUnitigMaxOverlaps.Failed
mv: `work1/overlap.coords' を stat できません: そのようなファイルやディレクトリはありません
terminate called after throwing an instance of 'std::runtime_error'はファイルを開くのに失敗したときに出るエラーらしい。
シグナル６は前段階がダメでkill されたのではなく自然消滅的にプロセスが終了したことを示すっぽい。
前段階のファイルが作られていない。それがoverlap.coordsってことなのかな。
どうしたものかな。
実はひそかに気になっていたのがJF_SIZEの設定で、今までの解析でゲノムサイズが747917875と推定されていたので、いままで設定していた値（JF_SIZE=3096134359：ソフトウェアにおすすめされた）が大きすぎるのではと思っていました。
#this is mandatory jellyfish hash size.-- a safe value is estimated_genome_size*20
次はREADMEにあったgenomesize*20に設定して見ようと思います。
JF_SIZE=15000000000
それから、ハッシュサイズをビビッて一番小さい値で設定していましたが、このままでは埒が明かないので最大化してみます。
ovlHashBits=25
ovlHashBlockLength=180000000
ドン！
masurca config.txt -o assemble.sh
エラーが出た。
Error line 38 of configuration file 'config.txt':
Invalid line 'ovlHashBits=25'
あびゃー
ovlHashBits=25を消してもう一度
Error line 38 of configuration file 'config.txt':
Invalid line 'ovlHashBlockLength=180000000'
あびゃー
CA_PARAMETERS = ovlHashBits=25
CA_PARAMETERS = ovlHashBlockLength=180000000
"CA_PARAMETERS ="の入れ忘れが原因だったみたいです。
さて、もう一周スタートしよう。
2020/12/23 今回は前回よりは先（Celera Assembler）まではいけたのだけど、やっぱり止まってしまいました
nohupの中身
[2020年 12月 23日 水曜日 17:26:06 JST] Recomputing A-stat for super-reads
/~i/software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 行 14: CA/5-consensus-coverage-stat/markRepeatUnique.err: そのようなファイルやディレクトリはありません
[2020年 12月 23日 水曜日 18:06:07 JST] Filtering overlaps
mv: `5-consensus-coverage-stat' を stat できません: そのようなファイルやディレクトリはありません
mv: `5-consensus-insert-sizes' を stat できません: そのようなファイルやディレクトリはありません
[2020年 12月 23日 水曜日 19:39:01 JST] Recomputing A-stat for super-reads
MultiAlignStore::MultiAlignStore()-- ERROR, didn't find any unitigs or contigs in the store.
MultiAlignStore::MultiAlignStore()-- asked for store 'CA/genome.tigStore', correct?
MultiAlignStore::MultiAlignStore()-- asked for version '5', correct?
MultiAlignStore::MultiAlignStore()-- asked for partition unitig=0 contig=0, correct?
MultiAlignStore::MultiAlignStore()-- asked for writable=0 inplace=0 append=0, correct?
/~/software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 行 14: CA/5-consensus-coverage-stat/markRepeatUnique.err: そのようなファイルやディレクトリはありません
[2020年 12月 23日 水曜日 19:48:04 JST] CA failed, check output under CA/ and runCA3.out
「recompute_astat_superreads_CA8.sh: 行 14: CA/5-consensus-coverage-stat/markRepeatUnique.err: そのようなファイルやディレクトリはありません」
これ以降がダメになっているくさい。でもそんなファイルもディレクトリもない！
せっかくCelera Assemblerまでいけたから、config.txtファイルからやりなおすんじゃなくてCA工程もう一周やってみようかな、と思ってassemble.shをCA部分だけ残して保存し直してもう一回かけてみた。
20201224 はい、ダメでした。
nohupファイルの中身
./assemble_CA.sh: 行 1: log: コマンドが見つかりません
./assemble_CA.sh: 行 1: log: コマンドが見つかりません
./assemble_CA.sh: 行 33: save: コマンドが見つかりません
./assemble_CA.sh: 行 35: save: コマンドが見つかりません
./assemble_CA.sh: 行 37: save: コマンドが見つかりません
./assemble_CA.sh: 行 39: save: コマンドが見つかりません
./assemble_CA.sh: 行 51: log: コマンドが見つかりません
./assemble_CA.sh: 行 57: log: コマンドが見つかりません
/~/software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 行 14: CA/5-consensus-coverage-stat/markRepeatUnique.err: そのようなファイルやディレクトリはありません
./assemble_CA.sh: 行 69: 予期しないトークン `(' 周辺に構文エラーがあります
./assemble_CA.sh: 行 69: `cat <(overlapStore -b 1 -e $NUM_SUPER_READS -d genome.ovlStore | awk '{if($1<$2 && ($1<'$NUM_SUPER_READS' && $2<'$NUM_SUPER_READS')) print $0}'|filter_overlap_file -t 20 <(gatekeeper -dumpfragments -withsequence genome.gkpStore| grep -P '^fragmentIdent|^fragmentSequence' | awk 'BEGIN{flag=1}{if(flag){print ">"$3}else{ print $3;} flag=1-flag; }') unitig_mers /dev/fd/0) <(overlapStore -d genome.ovlStore | awk '{if($1<$2 && ($1>='$NUM_SUPER_READS' || $2>='$NUM_SUPER_READS')) print $1" "$2" "$3" "$4" "$5" "$6" "$7}') |convertOverlap -ovl |gzip > overlaps_dedup.ovb.gz && '
あわわ。
消した部分にlogやsaveコマンドへのパスがあったっぽい。
しかし、多分一番の問題はやはり昨日と同じ。
5-consensus-coverage-statっていつ作られるファイルなんだ・・・？
assemble.shの中を見ると5-consensus-coverage-statという文字が一度だけ出てくる。
mv 4-unitigger 5-consensus 5-consensus-coverage-stat 5-consensus-insert-sizes genome.tigStore genome.ovlStore ovlStoreBackup &&
つまり、この前には既にこのフォルダはあるということだ。
nohupファイルによれば
Overlap/unitig success
ここまではいけて、
if [ ! -e CA/recompute_astat.success ];then
log "Recomputing A-stat for super-reads"
recompute_astat_superreads_CA8.sh genome CA $PE_AVG_READ_LENGTH work1/readPlacementsInSuperReads.final.read.superRead.offset.ori.txt superReadSequences_shr.frg
fi
ここが乗り越えられなかったので、この前の部分でトラブっているはず
で、しかもCAのまえにrm-rf CAしているので一旦CA内は空になっているはず
ということで犯人のプロセスは
rm -f CA/0-overlaptrim-overlap/overlap.sh CA/1-overlapper/overlap.sh CA/3-overlapcorrection/frgcorr.sh CA/3-overlapcorrection/ovlcorr.sh CA/5-consensus/consensus.sh
これになってしまうのだが、おかしい。これらが消えても5-consensus-coverage-statは消えなさそう
もう一回runCAの中を見ると
And finally, compute the coverage stat for all unitigs
#
my $l = getGlobal("utgGenomeSize");
if (! -e "$wrk/5-consensus-coverage-stat/computeCoverageStat.err") {
system("mkdir $wrk/5-consensus-coverage-stat") if (! -d "$wrk/5-consensus-coverage-stat");
$cmd = "$bin/computeCoverageStat \\\n";
$cmd .= " -g $wrk/$asm.gkpStore \\\n";
$cmd .= " -t $wrk/$asm.tigStore 5 \\\n";
$cmd .= " -s $l \\\n" if defined($l);
$cmd .= " -o $wrk/5-consensus-coverage-stat/$asm \\\n";
$cmd .= "> $wrk/5-consensus-coverage-stat/computeCoverageStat.err 2>&1";
このフォルダはすべてのunitigの計算が終わらないと作られない。
Cellera Assemblerの公式サイト(RunCA - wgs-assembler (sourceforge.net))を見てパラメータをいじってみることにした。
Unitigger
ユニティガーモジュールは、アセンブリパイプラインのステップ4として実行されます。このモジュールは、リードデータとメイトデータに加えて、パイプラインの上流で計算されたペアワイズオーバーラップを使用します。ユニティグは初期の高信頼度コンティグです。真のユニティグは、入力データの矛盾をゼロにします。ノイズの多いデータの中で大規模なユニティグを構築するために、Celera Assembler は入力データの矛盾が少ないユニティグのみを構築しようとします。
cnsMaxCoverage=7→3
#cnsMaxCovergae=float (default=0) (New in CA 8.2)
Limit read coverage during unitig consensus. Contained reads are randomly removed until the coverage limit is met.
Default 0 means unlimited coverage.
cnsReuseUnitigs=1→0
#cnsReuseUnitigs=integer (default=0) (New in CA 8.0)
For contigs that are formed from a single unitig, do not recompute the consensus sequence in stage 8-consensus. Instead, reuse the 5-consensus result. Note that VAR records will not be computed for these contigs.
utgGraphErrorLimit=1000→500
#utgGraphErrorLimit=float (default=3.25) (new in CA 7.0) Overlaps used for Best Overlap Graph construction in the bogart unitigger.
utgMergeErrorLimit=1000→500
#utgMergeErrorLimit=float (default=5.25) (new in CA 7.0)Overlaps below this threshold are Overlaps below this threshold are used for bubble popping and repeat detection in the bogart unitigger.
utgGraphErrorRate=0.015→0.030
#utgGraphErrorRate=float (default=0.030) (new in CA 7.0)
Overlaps below this threshold are used for Best Overlap Graph construction in the bogart unitigger.
utgMergeErrorRate=0.025→0.045
#utgMergeErrorRate=float (default=0.045) (new in CA 7.0)
Overlaps below this threshold are used for bubble popping and repeat detection in the bogart unitigger.
もういい加減これでダメならあきらめたいよ～～～
20201225
ブチギレてうっかり必要なファイルまで消してしまったので全部最初からやり直した
生成し直したassenble.shの上記数値を張り付け直したら、一部数値がconfig.txtの部分とかみ合わなくてエラー吐いて止まった
cgwErrorRate=0.1→0.15
これで動き出した
20201226
Cellera Assemblerまではいけたが、また止まった
nohup.out
[2020年 12月 26日 土曜日 16:15:11 JST] Celera Assembler
gkStore::gkStore()-- GateKeeper Store 'CA/genome.gkpStore' doesn't exist.
Search pattern not terminated at -e line 1.
Search pattern not terminated at -e line 1.
で、runCA.outを見たら、
runCA -s runCA.spec stopAfter=initialStoreBuilding -p genome -d CA doFragmentCorrection=0 ovlHashBlockLength=18000000
doOverlapBasedTrimming=0 doExtendClearRanges=0 ovlMerSize=30 superReadSequences_shr.frg pe.linking.frg 1> runCA0.out 2>&1
のovlHashBlockLength=18000000がダメと言われた
で、assemble.shの中をよく見たら
ovlHashBlockLength=10000000
save ovlHashBlockLength
ovlCorrBatchSize=`perl -e '$s=int('$TOTAL_READS'/100); if($s>10000){print $s}else{print "10000"}'`
save ovlCorrBatchSize
と指定している箇所があった
これと数値がかみ合ってなかったのと、これが先の文より後にあったのがまずかった？
ovlHashBlockLength=180000000
save ovlHashBlockLength
ovlCorrBatchSize=`perl -e '$s=int('$TOTAL_READS'/100); if($s>10000){print $s}else{print "10000"}'`
save ovlCorrBatchSize
echo "ovlRefBlockSize=$ovlRefBlockSize
ovlHashBlockLength=$ovlHashBlockLength
ovlCorrBatchSize=$ovlCorrBatchSize
を先の一文の前にかえて、あとはCellera assemplerより前の処理は消して別名でスクリプトを保存→再スタート
止まった・・・
File not found or invalid command line option 'ovlHashBlockLength'
File not found or invalid command line option '180000000'
'TOTAL_READS' is not a valid option; see 'runCA -options' for a list of valid options.
File not found or unknown specFile option line 'save TOTAL_READS'.
File not found or unknown specFile option line 'save ovlRefBlockSize'.
'' is not a valid option; see 'runCA -options' for a list of valid options.
元の位置に戻して再スタート
23：45 また止まった
File not found or invalid command line option 'ovlHashBlockLength'
File not found or invalid command line option '180000000'
'' is not a valid option; see 'runCA -options' for a list of valid options.
ovlHashBlockLengthを消すしかなさそう。
実行条件
ovlHashBlockLength=18000000
save ovlHashBlockLength
以外の部分を消して再スタート。
1：56
止まった
ovlRefBlockSize=24842374
ovlHashBlockLength=18000000
echo ovlRefBlockSize=24842374
./assemble1226.sh: 行 203: 一致する `"' を探索中に予期しないファイル終了 (EOF) です
./assemble1226.sh: 行 206: 構文エラー: 予期しないファイル終了 (EOF) です
たしかに206行目に無意味な”あったけど、スクリプトの一番最後の行で、解析はその段階まで行ってないはずなのでめちゃ謎。
20201228
また止まったし、またovlHashBlockLengthがダメらしい
usage: runCA -d <dir> -p <prefix> [options] <frg> ...
-d <dir> Use <dir> as the working directory. Required
-p <prefix> Use <prefix> as the output prefix. Required
-s <specFile> Read options from the specifications file <specfile>.
<specfile> can also be one of the following key words:
[no]OBT - run with[out] OBT
noVec - run with OBT but without Vector
-version Version information
-help This information
-options Describe specFile options, and show default values
<frg> CA formatted fragment file
Complete documentation at http://wgs-assembler.sourceforge.net/
File not found or invalid command line option '10000000'
ovlHashBlockLengthに関する部分全部消してもう一回・・・
20201228
今までにないくらい進んだ が、止まった
nohup.out の内容
[2020年 12月 28日 月曜日 01:14:43 JST] Celera Assembler
[2020年 12月 28日 月曜日 09:24:35 JST] Overlap/unitig success
[2020年 12月 28日 月曜日 09:24:35 JST] Recomputing A-stat for super-reads
/~/software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 12 行: 12187 中止 (コアダンプ) tigStore -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 4 -E unitig_cov.txt > tigStore.err 2>&1
//software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 13 行: 13422 中止 (コアダンプ) tigStore -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 5 -E unitig_cov.txt > tigStore.err 2>&1
/~/software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 14 行: 14225 中止 (コアダンプ) markRepeatUnique -short 300 -span 2 -lowcov 1 1 -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 5 -j 1 -k 5 -o ${ASM_DIR}/5-consensus-coverage-stat/genome.markRepeatUnique > ${ASM_DIR}/5-consensus-coverage-stat/markRepeatUnique.err 2>&1
[2020年 12月 28日 月曜日 10:29:37 JST] Filtering overlaps
Invalid uint64_t '-t' for [-s, --size=uint64]: Negative value
Could not malloc memory (4846228621574095687 bytes)
tigStore: src/AS_UTL/AS_UTL_alloc.C:65: void* safe_malloc(size_t): Assertion `p != __null' failed.
Failed with 'Aborted'
Backtrace (mangled):
tigStore[0x4140f2]
/lib64/libpthread.so.0(+0xf630)[0x7fa3f764c630]
/lib64/libc.so.6(gsignal+0x37)[0x7fa3f6f76387]
/lib64/libc.so.6(abort+0x148)[0x7fa3f6f77a78]
/lib64/libc.so.6(+0x2f1a6)[0x7fa3f6f6f1a6]
/lib64/libc.so.6(+0x2f252)[0x7fa3f6f6f252]
tigStore[0x4115b5]
tigStore[0x423798]
tigStore[0x423948]
tigStore[0x455551]
tigStore[0x40b418]
/lib64/libc.so.6(__libc_start_main+0xf5)[0x7fa3f6f62555]
tigStore[0x402759]
Backtrace (demangled):
[0] tigStore() [0x4140f2]
[1] /lib64/libpthread.so.0::(null) + 0xf630 [0x7fa3f764c630]
[2] /lib64/libc.so.6::(null) + 0x37 [0x7fa3f6f76387]
[3] /lib64/libc.so.6::(null) + 0x148 [0x7fa3f6f77a78]
[4] /lib64/libc.so.6::(null) + 0x2f1a6 [0x7fa3f6f6f1a6]
[5] /lib64/libc.so.6::(null) + 0x2f252 [0x7fa3f6f6f252]
[6] tigStore() [0x4115b5]
[7] tigStore() [0x423798]
[8] tigStore() [0x423948]
[9] tigStore() [0x455551]
[10] tigStore() [0x40b418]
[11] /lib64/libc.so.6::(null) + 0xf5 [0x7fa3f6f62555]
[12] tigStore() [0x402759]
GDB:
[2020年 12月 28日 月曜日 10:37:30 JST] Recomputing A-stat for super-reads
Could not malloc memory (4846228621574095687 bytes)
tigStore: src/AS_UTL/AS_UTL_alloc.C:65: void* safe_malloc(size_t): Assertion `p != __null' failed.
Failed with 'Aborted'
Backtrace (mangled):
tigStore[0x4140f2]
/lib64/libpthread.so.0(+0xf630)[0x7f7d59dd8630]
/lib64/libc.so.6(gsignal+0x37)[0x7f7d59702387]
/lib64/libc.so.6(abort+0x148)[0x7f7d59703a78]
/lib64/libc.so.6(+0x2f1a6)[0x7f7d596fb1a6]
/lib64/libc.so.6(+0x2f252)[0x7f7d596fb252]
tigStore[0x4115b5]
tigStore[0x423798]
tigStore[0x423948]
tigStore[0x455551]
tigStore[0x40b418]
/lib64/libc.so.6(__libc_start_main+0xf5)[0x7f7d596ee555]
tigStore[0x402759]
Backtrace (demangled):
[0] tigStore() [0x4140f2]
[1] /lib64/libpthread.so.0::(null) + 0xf630 [0x7f7d59dd8630]
[2] /lib64/libc.so.6::(null) + 0x37 [0x7f7d59702387]
[3] /lib64/libc.so.6::(null) + 0x148 [0x7f7d59703a78]
[4] /lib64/libc.so.6::(null) + 0x2f1a6 [0x7f7d596fb1a6]
[5] /lib64/libc.so.6::(null) + 0x2f252 [0x7f7d596fb252]
[6] tigStore() [0x4115b5]
[7] tigStore() [0x423798]
[8] tigStore() [0x423948]
[9] tigStore() [0x455551]
[10] tigStore() [0x40b418]
[11] /lib64/libc.so.6::(null) + 0xf5 [0x7f7d596ee555]
[12] tigStore() [0x402759]
GDB:
/~/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 9 行: 1954 中止 (コアダンプ) tigStore -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 5 -U -d layout > unitig_layout.txt
/~/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 12 行: 32782 中止 (コアダンプ) tigStore -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 4 -E unitig_cov.txt > tigStore.err 2>&1
/~/software/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 13 行: 33356 中止 (コアダンプ) tigStore -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 5 -E unitig_cov.txt > tigStore.err 2>&1
/~/anaconda3/envs/masurca/bin/recompute_astat_superreads_CA8.sh: 14 行: 34209 中止 (コアダンプ) markRepeatUnique -short 300 -span 2 -lowcov 1 1 -g ${ASM_DIR}/${PREFIX}.gkpStore -t ${ASM_DIR}/${PREFIX}.tigStore 5 -j 1 -k 5 -o ${ASM_DIR}/5-consensus-coverage-stat/genome.markRepeatUnique > ${ASM_DIR}/5-consensus-coverage-stat/markRepeatUnique.err 2>&1
[2020年 12月 28日 月曜日 11:47:54 JST] CA failed, check output under CA/ and runCA3.out
これによると
「Recomputing A-stat for super-reads」とやらからおかしくなっている
Overlap/unitig successとあるので今回はOverlap/unitigまではいけたっぽい
ただ、このあとおかしくなっているということはこの部分のアウトプットファイルが変である可能性が高い
ので、見てみた
unitigger.successファイルが空だった
runCA1.outの中身からこれが生成されるのは
/~/home/user/software folder/anaconda3/envs/masurca/CA8/Linux-amd64/bin/bogart -O /home/user/作業folder/CA/genome.ovlStore -G /home/user/作業folder/CA/genome.gkpStore -T /home/user/作業folder/CA/genome.tigStore -B 970405 -eg 0.030 -Eg 500 -em 0.045 -Em 500 -o /home/user/folder/CA/4-unitigger/genome > /home/user/作業folder/CA/4-unitigger/unitigger.err 2>&1
の工程らしいので、ここだけ切り離してコマンドにベタ貼りしてみた
やはりここで問題が生じているらしく、途中で止まった
~/anaconda3/envs/masurca/CA8/Linux-amd64/bin/bogart
usage: ~/anaconda3/envs/masurca/CA8/Linux-amd64/bin/bogart -o outputName -O ovlStore -G gkpStore -T tigStore
-O Mandatory path to an ovlStore.
-G Mandatory path to a gkpStore（アウトプットがあやしい）.
-T Mandatory path to a tigStore (can exist or not)（このファイルはなかった）.
-o prefix Mandatory name for the output files
-B b Target number of fragments per tigStore (consensus) partition
-B 970405 -eg 0.030 -Eg 500 -em 0.045 -Em 500 -o unitigger.err 2>&1
Algorithm Options
When constructing the Best Overlap Graph and Promiscuous Unitigs ('g'raph):
-eg 0.020 no more than 0.020 fraction (2.0%) error 誤差は２％まで
-Eg 2 no more than 2 errors (useful with short reads) 誤差は２つまで
When popping bubbles and splitting repeat/unique junctions ('m'erging):
-em 0.045 no more than 0.045 fraction (4.5%) error when bubble popping and repeat splitting 誤差は4.5％まで
-Em 4 no more than r errors (useful with short reads) 誤差は４つまで
20201231
考える気力のなさと修論のやばさとでサボっていました。
パラメータ設定変えたら少し進んだということはパラメータ設定をいじることには意味があるということなので、とりあえず、パラメータ設定について論文を調べてみた（遅くない・・・？）
が、そもそもこのアセンブラでのde novo アセンブリはバクテリア以外でやられてない説が浮上した。ガーン。ただ、公式サイトには特にバクテリアのみとは書かれていない。
考える気力/zero なので、ひとまず一個パラメータ設定をいじってみることにしました。
cnsMaxCoverage=7→3→10
何で増やしたのかと言うと、バクテリア論文がユニティグのカバレッジをかなり高めに設定していたからです（正しく読めていれば）。
さてやってみましょ。
