高スループットシークエンシング（HTS）は、費用対効果が高く時間効率の良いサンプルの完全な遺伝情報を得る能力を持ち、ゲノム研究に大きく革命をもたらした。多くの臨床応用において、作用可能な領域のパネルのみが調査対象である（Bellos et al、2014; Samorodnitsky et al、2015）（Diagnostic Panels. e.g., Takara ）。これらの分析では、研究者は、合成されたオリゴヌクレオチド（プローブ）のプールを用いて、ハイブリダイゼーションを用いて興味のあるゲノム断片を選択的に捕捉するターゲットキャプチャシーケンシングプロトコルを使用することが多い（Gnirke et al、2009）。これにより、ゲノムシーケンシング全体と比較して、より低いコストでより迅速な結果を得ることができ、臨床検査でのスケーラブルなアプローチが可能になる。
計算シミュレーションは、HTSデータ分析ツールの開発とベンチマーキングに不可欠である（Escalona et al、2016）。In silicoでのシミュレーションデータは、実際のデータよりも安価に生産することができる。それらは制御された条件下で生成され、完全に特徴づけることができる。さらに、シミュレーションは、研究者がシーケンシングプロトコルの性能を評価し、実験を実施する前に設計を最適化するのに役立つ。whole genome sequencing（Escalona et al、2016）とtargeted exome sequencing （Kim et al、2013）には数多くのシミュレータが利用可能だが、私たち（著者ら）の知る限りでは、現在のところキャプチャプロセスのダイナミクスをin silicoでシミュレートできるツールは存在しない。このようなツールは、計算解析パイプラインだけでなく、キャプチャデザインの効率を評価するのにも役立つと考えている。
この論文では、目標とするキャプチャシーケンシングデータをシミュレートする必要性を満たすためのソフトウェアパッケージであるCapSimを紹介する。プローブセットが与えられると、CapSimは、in silicoにおけるプローブハイブリダイゼーションのダイナミクスをシミュレートし、シーケンシングされるフラグメントセットを生成する。既存のHTSシミュレーションツールとは異なり、CapSimは、フラグメンテーション、フラグメントキャプチャ、シーケンシングなど、シーケンス処理のあらゆる段階をエミュレートする。 CapSimはJavaで書かれており、あらゆるコンピューティングプラットフォーム上でネイティブに実行できるため、バイオインフォマティクスコミュニティに簡単にアクセスできる（どのような環境のバイオインフォマティシャンでも利用できる）。
シーケンシングプロセスはDNAの断片化から始まる。CapSimは所定のサイズ分布を有するゲノム配列からフラグメントを繰り返しサンプリングしてこのプロセスをシミュレートする。いくつかのデータセットからの断片サイズ分布に最もよく適合することがわかったlog-logistic分布を用いて断片長をモデル化する（論文の補足資料参照）。
ターゲットキャプチャシーケンシングの次のステップでは、DNA断片にハイブリダイゼーションにした標的プローブとそのDNA断片はビーズからpull downされる。ビーズは標的プローブとハイブリダイズしたDNA断片と特異的に結合している。CapSimはプローブをDNA断片にマッピングすることでこのプロセスをシミュレートする。計算上効率的であるために、CapSimは最初にプローブをゲノム配列にマップし、DNA断片がゲノムからサンプリングされると、greedy algorithmを使用してDNA断片に結合できるプローブの最大数を決定する。我々（著者ら）は、キャプチャされる確率がバインドしたプローブの数に比例し、断片長に反比例するキャプチャプロセスの確率論的性質をモデル化した。ハイブリダイゼーションのダイナミクスのこのシミュレーションは、Kim et al. (2013)よりもより現実的なキャプチャシーケンシングデータを生成することが示されている。 Kim et al. (2013)は、キャプチャシークエンシングをシミュレーションするための唯一のツールである（論文の結果参照）。
次いで、キャプチャされたフラグメントはin silico sequencingされる。 Illuminaシーケンシングのために、CapSimはクラスター形成のDNA断片分布を導入ししている。 それはLog-Logistic分布を使用して、キャプチャステップからシミュレートされたDNA断片からサンプリングを行う（論文の補足資料を参照）。これらの選択されたDNA断片は、DNA断片の両末端からエラーありで配列が読み取られペアエンドシークエンシングシミュレーションに使用される。 PacBioシークエンシングのために、CapSimは、所与の分布からのポリメラーゼが読む長さをシミュレートする（補足資料参照）。DNA断片トのシークエンシングをシミュレートする際、CapSimはポリメラーゼが読める長さに達するまで、PacBioエラープロファイルを使い、2つの鎖のコピーを交互に行う。
BioinfomaticsのSEQUENCE ANALYSISに発表された論文です。大半のデータはsupplementaryにあります。
capsimに関するツイート。
Japsa（Just Another JAva Package for Sequence Analysis）マニュアル
http://japsa.readthedocs.io/en/latest/install.html
ダウンロード
mac os 10.13、java1.6環境でテストした。
リリースよりjapsaをダウンロードする。
https://github.com/mdcao/japsa/releases
tar zxvf JapsaRelease.tar.gz
cd JapsaRelease
./install.sh
コンソールに表示される指示に従いインストールする。最後の質問のHDF libraryがなければ"none"と打つ。ここではbin/にパスを通した。
cd bin/
echo 'export PATH=$PATH:`pwd`' >> ~/.bash_profile && source ~/.bash_profile
> jsa.sim.capsim --help
$ jsa.sim.capsim --help
Simulate capture sequencing
Usage: jsa.sim.capsim [options]
Options:
--reference=s Name of genome to be
(REQUIRED)
--probe=s File containing probes mapped to the reference in bam format
(default='null')
--logFile=s Log file
(default='-')
--ID=s A unique ID for the data set
(default='')
--miseq=s Name of read file if miseq is simulated
(default='null')
--pacbio=s Name of read file if pacbio is simulated
(default='null')
--fmedian=i Median of fragment size at shearing
(default='2000')
--fshape=d Shape parameter of the fragment size distribution
(default='6.0')
--num=i Number of fragments
(default='1000000')
--pblen=i PacBio: Average (polymerase) read length
(default='30000')
--illen=i Illumina: read length
(default='300')
--seed=i Random seed, 0 for a random seed
(default='0')
--help Display this usage and exit
(default='false')
——
ラン
genepanelのfasta配列が見つからなかったので、biomartで生成する。ここでは、指定した領域から配列を出力している（gene idを持っているなら、biomartでidから配列を生成する方が王道。IDを別種のID、または遺伝子名に変換したいならDavidを使う）。
１、配列の準備
Gene panelの領域は、Perkin Elmer: Amplicon Panel Kitsの領域を使う。
http://www.biooscientific.com/Illumina-Amplicon-Panels
ここではColorectal Cancer -1 kitのターゲット領域を選んだ。
http://www.biooscientific.com/NEXTflex-Colorectal-Cancer-1-Amplicon-Panel
Ensembl - BIomartにアクセスする。
http://asia.ensembl.org/biomart/martview/f5938f5b76d2d943d8c4a176709e376a
Ensembl Genes 93を選択、Human genome GRCh38を選択。
左のタブからFilterを選択。様々な条件で配列などを抽出できる。代表的なデータベースのIDからも配列は抽出できるが、ここではbedファイルしかないので、ポジション情報から配列を抽出する。Multiple regionにbedの配列をペーストする。bedファイルのタブはあらかじめ":"に置換してある。
貼ったところ。
左のメニューからAttributeタブを選択。SequencesとcDNA sequencesを選択。
SequencesとcDNA sequencesを選んだところ。
左上のResultsをクリックして、配列を生成する。Compressed fileをダウンロードする。
ターゲット領域のMulti-FASTAファイルが準備できた。テストランでは、これをprobes.faとして使う。 capsimのランに移る。
２、capsim実行。
probes.faをリファレンスにmappingする。GRCh38はEnsembl humanからダウンロードした。
bowtie2-build ref.fa ref
bowtie2 --local --very-sensitive-local --mp 8 --rdg 10,8 --rfg 10,8 -k 10000 -f -x ref -U probes.fa -S probes.sam
samtools sort -O BAM probes.sam > probes.bam && samtools index probes.bam
capsimを実行。
#pacbio
jsa.sim.capsim --reference ref.fa --probe probes.bam --ID someid --fmedian 5000 --pacbio output --pblen 3000 --num 20000000
#Miseq
jsa.sim.capsim --reference ref.fa --probe probes.bam --ID someid --fmedian 500 --miseq output --illen 300 --num 20000000
- --ID A unique ID for the data set (default='')
- --fmedian Median of fragment size at shearing (default='2000')
- --pacbio Name of read file if pacbio is simulated (default='null')
- --pblen PacBio: Average (polymerase) read length (default='30000')
- --num Number of fragments (default='1000000')
- --miseq Name of read file if miseq is simulated (default='null')
- --illen Illumina: read length (default='300')
マッピング結果。
figv -g ref.fa miseq.bam,Colorectal_Cancer_-1.bed
引用
Simulating the dynamics of targeted capture sequencing with CapSim
Cao MD, Ganesamoorthy D, Zhou C, Coin LJM
Bioinformatics. 2018 Mar 1;34(5):873-874.
参考
TruSeq Amplicon Cancer Panel
