研究内容を高めて（深めて）いくために大事なこと

随時追加します。Tweetはクリックするとスレッドが読めます。

Here are my 12 guidelines for data exploration and analysis with the right attitude for discovery:
1. You never really finish analyzing a dataset. You just decide to stop and move on at some point, leaving some things undiscovered. 🧵

— Itai Yanai 💔 (@ItaiYanai) 2023年1月10日

How to become a SUCCESSFUL academic: a guide 1/n

— Maarten van Smeden (@MaartenvSmeden) 2021年1月16日

This is how I advise my #PhD students to write research manuscripts (in case someone finds it helpful).

General points:

1. Research questions addressed by your manuscript are key and should guide you.

2. Don’t view your manuscript as an article. See it as a STORY.

3. Pick the…

— Andrew Akbashev (@Andrew_Akbashev) 2023年8月29日

In the LLM-science discussion, I see a common misconception that science is a thing you do and that writing about it is separate and can be automated. I’ve written over 300 scientific papers and can assure you that science writing can’t be separated from science doing. Why? 1/18

— Michael Black (@Michael_J_Black) 2022年12月3日

1. Inspired by @BiophysicalFrog insightful comments about the lack of training in science in how to be a PI I wanted to share some thoughts about something that helped me a lot over the years in my own learning journey. This is a thread about PI Superpowers!

— Tanentzapf Lab (@TanentzapfLab) 2019年3月13日

• 技術者が研究職をするうえで考えること

biochem-fan.hatenablog.com

Docker daemonが ルートディレクトリの容量を圧迫するのを回避する

Dockerのイメージを大量にビルドなどしていると、気づくとルート(/)の容量がいっぱいになっていることがあります。

$ df -h
ファイルシス          サイズ  使用  残り 使用% マウント位置
udev                    126G     0  126G    0% /dev
tmpfs                    26G  2.5M   26G    1% /run
/dev/sda3               275G  273G     0  100% /    ←100%になっている

フォルダ容量を調べてみると、 /var/lib/docker/overlay2 が大きく容量を取っていることがわかります。

/var/lib/docker# du -sh *
88K     buildkit
120K    containers
266M    image
64K     network
211G    overlay2
16K     plugins
4.0K    runtimes
4.0K    swarm
4.0K    tmp
4.0K    trust
28K     volumes

既存のイメージを削除していく方法はいくつか見つけました。

dockerで/var/lib/docker/overlay2 が肥大化した時の対処 - 1 Minute Tech Tips

/var/lib/docker/overlay2のストレージ肥大化対応 - /home/jiro4989

ここでは /var/lib/docker/overlay2 の容量を確保し、いくらdocker buildしても大丈夫な環境を構築したいと思います。

物理パーテションには空きはないので、ルートパーテションを確保するにはOSの再インストールから行わないといけません。それはやりたくない。 そこで、容量に余裕がある他のディレクトリ（ここではホームディレクトリ (/home/rnakato)とする）にディレクトリを移動し、そこにリンクを貼るやり方を試してみます。

# systemctl stop docker
# mv /var/lib/docker/overlay2  /home/rnakato/var_lib_docker_overlay2/
# ln -s /home/rnakato/var_lib_docker_overlay2/  /var/lib/docker/overlay2
# systemctl start docker

しかしこの方法だとうまくいきません。docker buildを使うと以下のようなエラーが出ます。 どうやらシンボリックリンクディレクトリだとdocker は動作しないようです。

failed to create shim: OCI runtime create failed: runc create failed: invalid rootfs: not an absolute path, or a symlink: unknown

Docker 'is not an absolute path or is a symlink' error - Stack Overflow

そこで、ホームディレクトリに新しく作ったディレクトリを mount コマンドでbindしてみることにします。

# systemctl stop docker
# mv /var/lib/docker/overlay2  /home/rnakato/var_lib_docker_overlay2/
# mkdir  /var/lib/docker/overlay2
# mount --bind  /home/rnakato/var_lib_docker_overlay2/  /var/lib/docker/overlay2
# systemctl start docker

Dockerの/var/lib/dockerを移動する - Qiita

この方法だとうまくいきました。

$ df -h
ファイルシス          サイズ  使用  残り 使用% マウント位置
udev                    126G     0  126G    0% /dev
tmpfs                    26G  2.6M   26G    1% /run
/dev/sda3               275G   63G  199G   24% /

NGS解析のための共有サーバ環境構築を考える（2022年度版）（2022/9/12追記）

5年前にこんな記事を書きました。 この時はPython2系と3系の共存（というか、Anaconda環境構築のベストプラクティス）に苦労していました。

rnakato.hatenablog.jp

今はほとんどPython3系のみで事足りますが、依然としてdependenciesの問題は解決されていません。 特に、bioconda から多数のツールをインストールしようとするとconflictが多発します。conda-forge と比べて、登録されたツールのdependencyがきちんと管理されていないような印象です。 また、Pythonの公式ライブラリも後方互換をあまり重視していない感じであり、Python3.7で書かれたツールをPython3.9にインストールしても動かないというようなこともあります。

私のようなサーバ管理者からすると、不特定多数の一般ユーザが共有できる解析環境をサーバ上に構築・管理したいものです。 うちはローカルサーバですが、クラウドでも基本的に変わらないでしょう。 色々な管理ツールを使ってツール群を管理していますが、Pythonは特に維持が大変です。

Rは管理者と一般ユーザでインストールされるディレクトリが分かれるのでわかりやすいですね。 Pythonの場合も共有ディレクトリにanacondaやpyenvをインストールしてパスを通すことで環境を共有することは可能です。 管理者ユーザでインストールすれば、一般ユーザに環境を書き換えられることもありません。 ただしRと異なり、このやり方だと一般ユーザが自分でツールをインストールすることができません。 毎回管理者にpipやcondaのインストールをお願いするということになります。これはお互いにとって負担です（管理者ユーザが十分数いれば別ですが）。

（2022/9/12追記）anacondaでも、管理者ユーザでインストールした環境にユーザがpipでインストールした場合は~/.local/binにインストールされるようです。そういえば pip --userってオプションがありましたね。

代替案として、各ユーザが自分のホームディレクトリにanacondaをインストールし、自力で管理するという方式が考えられます。 この場合、「自分でanacondaや各種ツールをインストール・管理できるユーザである」ことがサーバ利用の必要条件になります。 しかし実際にはそういうユーザばかりではありません。とにかく教えられた・決まった作業しかできないユーザもいます。 そのようなユーザは解析を行う資格なし、と厳しいポリシーで運用するというのは一案ですが（管理者は楽ですが）、実際にはそうも言っていられません。 ラボに入った学生は卒業させなければいけませんし、解析を習いに来たウェット研究者に環境構築から行わせていたらそれだけに日数が取られ、それは本意ではありません。

また、正しく運用したとしても上述のconflictは解決しません。 condaとpipを共存させてはいけないというのはよく言われることですが、condaのみに限ってもconflictは不可避です。 特にbiocondaはかなりの鬼門です。今の私のローカル環境にはもうdeepToolsはインストールできません。

まだ問題があります。「pipとcondaは共存できない」時にどちらかひとつに限るとしたならばおそらくcondaになるのですが、condaでインストールできるツールのバージョンはしばしば古いものになります（aptと似たような問題）。最新バージョンをインストールするためにpipを使いたいという場合があり得ます。 また、pipと違ってcondaはdependencyとしてRやらperlやらも一緒にインストールしてしまう「行儀が悪い」インストーラなので、既に構築した既存の解析環境を上書きし、破壊する危険性があります。様々なユーザが様々な用途で様々なツールをインストールしていくことになると、完全に制御することはほとんど不可能です。

これらの問題回避策としてよく利用されるのが、condaまたはvenvで仮想環境を構築するという方法です（「conda 仮想環境 」でググるとたくさん出てきます）。 たとえばreadのアダプタ配列をトリミングする cutadaptのインストールページ を見ると、仮想環境にcutadaptを入れる方法が例示されています。 現状では、この仮想環境構築がベストソリューションのように感じます。試してみましたが、仮想環境もちゃんとユーザ間で共有できるようです。 この方式では、base環境には必要最低限の（たとえばAnaconda同梱の）パッケージのみにとどめ、各解析ツールは全て個別の仮想環境に切り分ける、という運用が考えられます。

ただ、ここは個人的な意見になりますが、「仮想環境の切り替え」は一般ユーザにはややハードルが高いように感じます（私が仮想環境の切り替えになじんでいないだけかもしれませんが）。 cutadapt, deeptools, etc.. の個別のツールを使うたびに仮想環境切り替えするのは個人的にはわずらわしく感じます。シェルスクリプトのコーディングも書きづらいような。各仮想環境のアップデートも面倒なのでは…？

また、複数のサーバを管理している場合、各サーバでの環境統一が難しいという問題もあります（これはクラウドサーバにすれば解決する問題かもしれません）。解析環境を凍結保存したり、再構築することも難しいかもしれません。

ここで出てくるのが、切り分けしたい環境をDocker/Singularityのイメージとして保存するという方法です。 やろうとしていることはcondaの仮想環境構築と大差ないのですが、singularityイメージとして目に見えるかたちで保存することで、管理がしやすくなります。複数サーバ間での同期も簡単ですし、バージョン番号をファイル名につけておけば過去の環境をキープしつつの更新も可能です。

さらに、この方式の場合はPythonに限らず、あらゆる言語のツールをパッケージとして保存することが可能になり、NGSツール群を統一的に管理することができるようになります。 また、ツールを使うためのスクリプトを同梱しておくことでユーザビリティを高めることもできます。 例えば、ChromHMMのjarファイルと一緒に java -mx40000M -jar <path>/ChromHMM.jar "$@" と記載したシェルスクリプトを作成しパスを通しておくことで、一般ユーザはjavaファイルであることを意識せずにChromHMMが使えるようになりますし、将来的に起こりうるメモリエラーも未然に防ぐことができそうです（condaパッケージもおそらく同じようなことを目指しているのでしょう）。

ただこの場合でも、「どのくらい細かく環境を切り分けるか」という問題は残ります。 個々のツールを別々のイメージにする「１ツール１イメージ」方式が一つの案ですが、一方で、「様々な解析が可能であるイメージ」も需要が高いです。 様々なツールを内部で使うパイプラインのようなケースです。 ひとつのアイデアは、NextFlowやSnakemakeなどのワークフロー言語で全体のワークフローを記述し、「１ツール１イメージ」のイメージ群を上手に使いこなしてパイプラインを完遂する方法です。ただ、このやり方は基本的には決まったワークフローで処理することが主な目的であり、自分の研究のためにあれこれ工夫したり拡張したりという目的には向きません。

「国際プロジェクトで標準となった解析フローを実行するためのDocker image」や「様々なツールをインストール済のシングルセル解析プラットフォーム」のような、オールインワン方式のイメージを作りたい場合、ではその内部でconflictをどう回避するのか？という最初の問題に戻ります。 Dockerコンテナを起動して内部で解析する場合、通常の解析環境と同じように、仮想環境の切り替えで対応可能です（ただし「仮想環境敷居が高い問題」は回避されません）。単一のツールを仮想環境でインストールする場合、その仮想環境へのパスを直接通してしまうことでツールを使えるようになります。スクリプトを用意することで、ユーザの負担なくさまざまなツールを使えるようにすることはある程度可能でしょう。Singularityでコマンド実行する場合は仮想環境の切り替えはおそらく難しいので、そのような方式になります。いずれにしても、複数の仮想環境を包含するようなイメージを維持可能かどうかは程度問題であり、個別のケースに応じて工夫を凝らすしかないと言えそうです（例外として、Jupyter notebookを使う場合は多数の仮想環境を使い分けることが可能になります）。 さらなる代替案として、「多数のsingularityイメージを内部に保存したsingularityイメージ」というのを思いつきましたが…そもそも内部に多数の個別環境を内包するイメージはコンテナ化の趣旨からは外れたものであるかもしれません。

結局、ある程度ユーザを自力で啓蒙しつつ、ある程度ローカル環境で対応しつつ、ある程度個別環境を用意する、という道を模索するしかなさそうです。 ５年後にはこのあたりの微妙な感じを全て打ち払うベストソリューションが生まれるとよいのですが。

chromapを試す（その２）

高速にマッピングできるchromap*1 の続きです。今回はHi-Cデータに対してマッピングし、BWAと比較します。

前回の記事はこちら

rnakato.hatenablog.jp

Hi-Cファイルの取得

今回はHaarhuisらのCell論文*2 のデータを使いました。ヒトHap1細胞でのHi-Cデータです。

fastq-dump --split-files --gzip SRR5266584

BWAでマッピング

BWAを使ってfastqデータをマッピングします。マッピングのパラメータは 4D Nucleome project *3 を参考にしています。用いるCPU数は32です。

$ time bwa mem -t 32 -SP5M bwa-indexes/UCSC-hg38 \
   SRR5266584_1.fastq.gz SRR5266584_2.fastq.gz \
   > mapped.bwa.sort.bam
... (中略)
[M::mem_pestat] analyzing insert size distribution for orientation RR...
[M::mem_pestat] (25, 50, 75) percentile: (1778, 3622, 6305)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 15359)
[M::mem_pestat] mean and std.dev: (4152.72, 2755.39)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 19886)
[M::mem_process_seqs] Processed 604710 reads in 302.831 CPU sec, 9.743 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 32 -SP5M bwa-indexes/UCSC-hg38 SRR5266584_1.fastq.gz SRR5266584_2.fastq.gz
[main] Real time: 1874.166 sec; CPU: 56598.868 sec

________________________________________________________
Executed in   31.24 mins   fish           external
   usr time  928.66 mins  595.00 micros  928.66 mins
   sys time   14.65 mins  231.00 micros   14.65 mins

32コアで31分程度で終わりました。実際、マッピングのみであればHi-Cであってもそこまで遅くはありません。

実際のHi-C解析では、以下のようにpairtools を用いてその後のマップリードのソートや冗長リードの除去などに更に時間が必要になります。コマンドの詳細は割愛しますが、ここでは72分かかっています。31分 + 72分 = 103分を使っています。

$ time pairtools parse --min-mapq 30 \
   --walks-policy 5unique --max-inter-align-gap 30 \
   --nproc-in 8 --nproc-out 8 --chroms-path genometable.hg38.txt mapped.bwa.bam \
 | pairtools sort --nproc 8 \
 | pairtools dedup --nproc-in 8 --nproc-out 8 --mark-dups \
 > mapped.bwa.pairsam

real    72m54.795s
user    79m7.940s
sys     17m2.803s

chromapでマッピング

chromapでマッピングしてみます。Hi-Cデータのため、--preset hicオプションをつけています。これにより、出力形式が4DNで用いられているpairs形式 になります。 また、--remove-pcr-duplicatesオプションで冗長リードの除去ステップを追加しています。

なお、chromapではマッピングの時にもゲノムデータ (hg38.fa)を指定する必要があります。

$ time chromap --preset hic --remove-pcr-duplicates \
   -t 32 -x chromap-indexes/UCSC-hg38 -r hg38.fa \
   -1 SRR5266584_1.fastq.gz -2 SRR5266584_2.fastq.gz \
   -o mapped.chromap.pairs
Preset parameters for Hi-C are used.
Start to map reads.
Parameters: error threshold: 4, min-num-seeds: 2, max-seed-frequency: 500,1000, max-num-best-mappings: 1, max-insert-size: 1000, MAPQ-threshold: 1, min-read-length: 30, bc-error-threshold: 1, bc-probability-threshold: 0.90
Number of threads: 32
Analyze bulk data.
Won't try to remove adapters on 3'.
Will remove PCR duplicates after mapping.
Will remove PCR duplicates at bulk level.
Won't allocate multi-mappings after mapping.
Only output unique mappings after mapping.
Only output mappings of which barcodes are in whitelist.
Allow split alignment.
Output mappings in pairs format.
Reference file: hg38.fa
Index file: chromap-indexes/UCSC-hg38
1th read 1 file: SRR5266584_1.fastq.gz
1th read 2 file: SRR5266584_2.fastq.gz
Output file: mapped.chromap.pairs
Loaded all sequences successfully in 28.22s, number of sequences: 25, number of bases: 3088286401.
Kmer size: 17, window size: 7.
Lookup table size: 392946789, occurrence table size: 445474231.
Loaded index successfully in 112.23s.
Mapped 500000 read pairs in 5.55s.
（中略）
Mapped 500000 read pairs in 1.56s.
Mapped 100349 read pairs in 0.56s.
Mapped all reads in 304.95s.
Number of reads: 120200698.
Number of mapped reads: 108983560.
Number of uniquely mapped reads: 100363412.
Number of reads have multi-mappings: 8620148.
Number of candidates: 4829748453.
Number of mappings: 108983560.
Number of uni-mappings: 100363412.
Number of multi-mappings: 8620148.
Sorted, deduped and outputed mappings in 57.58s.
# uni-mappings: 53267577, # multi-mappings: 868789, total: 54136366.
Number of output mappings (passed filters): 47034561
Total time: 508.85s.

________________________________________________________
Executed in  511.87 secs   fish           external
   usr time  100.27 mins  475.00 micros  100.27 mins
   sys time    1.93 mins  234.00 micros    1.93 mins

なんと、32コアで8分あまりでマッピング終了しました。BWAでは冗長リードの除去が別途pairtools で必要になることを考えると、これは相当速いです。

ただしchromapは制限酵素サイトを考慮するオプションを持っていないので、フラグメントのフィルタリングに別途pairtoolsなどを利用する必要があります*4。手元で試したところでは、chromapの出力をpairtoolsに入れるとエラーになるため、更なる調査が必要そうです。Micro-Cの場合は制限酵素サイトを考慮する必要がないので、この出力ファイルをそのまま利用することが可能です。

まとめ

前回のChIP-seqサンプルだとあまり速くなったと感じませんでしたが、Hi-Cサンプルだとその高速性がはっきりしました。

ここで使ったサンプルはどちらかといえばread depthの少ないサンプルですので、よりdeepなHi-Cデータをマップする時にはchromapの高速性は更に真価を発揮しそうです。 制限酵素サイトへの対応がまたれるところです。

*1:H. Zhang et al., Fast alignment and preprocessing of chromatin profiles with Chromap, Nature Communications, 2021, DOI: 10.1038/s41467-021-26865-w

*2:J. Haarhuis et al., The Cohesin Release Factor WAPL Restricts Chromatin Loop Extension, Cell, 2017, DOI: 10.1016/j.cell.2017.04.013

*3:https://www.4dnucleome.org/

*4:https://github.com/haowenz/chromap/issues/90

Linux ワンライナー覚書

スクリプトなどを用いずにLinuxのCUIコマンドのみで一行で行う処理（またはその処理を行うコマンド）のことをワンライナーと呼びます。 ここではゲノム解析に使うワンライナーを順次追加していきます。

一般的なファイル操作

 $ ls workdir/            # workdirに含まれるファイルのチェック
 $ ls workdir/  | less  # ファイル数が多すぎる時はパイプしてlessで表示
 $ ls workdir/ | wc -l  # workdirに含まれるファイルの数を表示
 $ du -sh workdir/     # workdirディレクトリの容量チェック
 $ du -sh workdir/*   # workdirに含まれるファイルの容量チェック
 $ ls workdir/*.bed | wc -l  # workdirに含まれるBEDファイルの数チェック
 $ grep -x -f genelist1.txt genelist2.txt   # 2つのファイルの共通行を表示
 $ sed -i 's/aaa/bbb/g' sample.txt  # sample.txt 内の文字列aaaをbbbに置換
 $ sed -i 's/aaa/bbb/g' *.txt  # カレントディレクトリ内の全ての*.txt の文字列aaaをbbbに置換

BED, refFlatファイルなどのゲノムファイルの操作

BED

 # genome.faに含まれる染色体の表示
 $ grep \> genome.fa
 # BEDファイルをsorted BEDに変換
 $ sort -k1,1 -k2,2n peak.bed > peak.sort.bed
 # BED6ファイルをBED3に変換
 $ cut -f1-3 peak.bed6 > peak.bed3
 # peak.bedのピーク数をカウント（ヘッダが含まれる場合は-1）
 $ wc -l peak.bed 
 # peak.bedに含まれるピークを染色体別にカウント
 $ cut -f1 peak.bed | sort | uniq -c 
 # peak.bedからchr1のピークのみを取り出す
 $ cat peak.bed | awk 'if($1 == “chr1"){ print }' > peak.chr1.bed 
 # BEDファイル(タブ区切り)を２列目、１列目、３列目の順に表示
 $ cat peak.bed | awk -v 'OFS=\t' '{print $2, $1, $3}' 
 # カレントディレクトリ内にある全てのBEDファイルの染色体名 chrmt を chrM に置換
 $ sed -i 's/chrmt/chrM/gI' *.bed

refFlat

 # refFlatファイルをBED3形式に変換
 $ cat gene.refFlat | awk 'BEGIN{OFS = "\t"}{ print($3, $5, $6)}' > gene.bed 
 # refFlatファイルをBED6形式に変換
 $ cat gene.refFlat | awk 'BEGIN{OFS = "\t"}{ print($3, $5, $6, $1, ".", $4)}' > gene.bed
 # refFlatファイルからForward strandの遺伝子のみを取り出す
 $ cat refFlat.txt | awk 'if($4 ==+){ print }' > refFlat.forward.txt
 # refFlatファイルからReverse strandかつchr2の遺伝子のみを取り出す (複数条件指定)
 $ cat refFlat.txt | awk '{if($3 == “chr2” && $4 ==“-”){ print }}'  > refFlat.reverse.chr2.txt
 # refFlatファイルに含まれる遺伝子のうち最大のexon数を表示
 $ cat refFlat.txt | grep -v exon | awk 'BEGIN{ a=0 }{ if(a<$9) a=$9} END{ print a}'
 # refFlatファイルに含まれる遺伝子の平均exon数を表示
 $ cat refFlat.txt | grep -v exon | awk 'BEGIN{ a=0 }{ a+=$9 } END{ print a/NR }'
 # 平均exon数を整数で表示したい場合
 $ cat refFlat.txt | grep -v exon | awk 'BEGIN{ a=0 }{ a+=$9 } END{ printf "%d\n", a/NR }'

中級者以上向け

 # ファイル名にハイフンを含むBEDファイルのハイフンをアンダーバーに一括変換
 $ for file in `ls *-*.bed`; do mv $file ${file/-/_}; done
 # カレントディレクトリのsraファイルを全てfastq-dump
 $ find . -name '*.sra' -exec fastq-dump {} \; 
 # xargsを使って上と同じ処理
 $ ls *.sra | xargs fastq-dump
 # カレントディレクトリ以下にある全ての.txtファイルに対してaaaからbbbへの置換を実行（元のファイルは*.bakとして保存される）
 $ find ./ -name '*.txt' | xargs sed -i.bak "s/aaa/bbb/g"  
 

解析ツールを利用したワンライナー

# ヒトゲノム (hg38)の染色体長を表すBEDファイルを作成
$ fetchChromSizes hg38 | awk -vOFS="\t" '($1!~/_/){ print $1, "0", $2 }' \
   > chrlen.hg38.bed
# merged siteが各ピークファイルといくつ重複しているかを調べる
$ intersectBed -wa -a peak1.bed -b peak2.bed -wb -filenames \
  | cut -f1,2,3,4 | sort | uniq \
  | cut -f1,2,3 | uniq -c \
  | awk '{sub(/^[ ]+/, "")}' \
  | awk '{print $1}' | sort | uniq -c

chromapを試す（その１）

chromapは高速・高精度にリードをゲノムにマップするツールです*1。 既存ツールのbwaやbowtie2よりも遥かに高速だとのことなので、ここではChIP-seqデータを用いて使用法及び速度・精度面を確かめてみたいと思います。

chromapについて

chromap は minimap2*2の後継です。NGS解析では通常リードがどこにどのくらいマップされたかのカウントだけを見るため、１塩基単位での詳細なアラインメントは必要なく、その部分を省くlight-weight alignmentという方法を用いることで、bwaやbowtie2よりも数倍高速にマッピングを達成しています（ざっくりした説明）。 minimap2は50bp以上のseed領域の一致を必要とするため、scATAC-seqのようにバーコード配列を含むリードのマッピングには不向きであるという弱点がありました。chromapはその弱点を克服し、アラインメントだけでなく冗長リードの除去やアダプター配列の同定までを高速かつ包括的に行います。

ChIP-seqデータはシングルエンドかつアダプター配列も含まれていないのでchromapの真価は十分に発揮されないのですが、簡単のため（あと私が一番慣れているため）フィージビリティスタディとしてここでは利用します。

インストール

condaでインストールするのが簡単です。

conda install -c bioconda chromap

indexの作成

ゲノム配列を入力にしてindexを作成します（ここではhg38.fa）。-i がindex作成のためのオプションですね。 ゲノムファイルの作成方法はこちらを参照してください。

chromap -i -r hg38.fa -o chromap-UCSC-hg38

ChIP-seqファイルの取得

なんでもよいのですが、ここではPol2 ChIP-seq (K562)のデータをSRAよりダウンロードしました。SRR227359.fastq.gzという名前で保存されます。

fastq-dump --gzip SRR227359

マッピング

chromapでhg38にマッピングしてみましょう。わかりやすさのため、打ち込むコマンドの前には$マークをつけています。 chromapではマッピングの時にもゲノムデータを指定する必要があります。用いるCPU数は24です。

$ time chromap --preset chip -t 24 \
$ -x chromap-UCSC-hg38 -r hg38.fa \
$ -1 SRR227359.fastq.gz -o SRR227359.chromap.bed
Preset parameters for ChIP-seq are used.
Start to map reads.
Parameters: error threshold: 8, min-num-seeds: 2, max-seed-frequency: 500,1000, max-num-best-mappings: 1, max-insert-size: 2000, MAPQ-threshold: 30, min-read-length: 30, bc-error-threshold: 1, bc-probability-threshold: 0.90
Number of threads: 24
Analyze bulk data.
Won't try to remove adapters on 3'.
Will remove PCR duplicates after mapping.
Will remove PCR duplicates at bulk level.
Won't allocate multi-mappings after mapping.
Only output unique mappings after mapping.
Only output mappings of which barcodes are in whitelist.
Output mappings in BED/BEDPE format.
Reference file: hg38.fa 
Index file: chromap-UCSC-hg38
1th read 1 file: SRR227359.fastq.gz
Output file: SRR227359.chromap.bed
Loaded all sequences successfully in 50.49s, number of sequences: 25, number of bases: 3088286401.
Kmer size: 17, window size: 7.
Lookup table size: 392946789, occurrence table size: 445474231.
Loaded index successfully in 174.56s.
Loaded 500000 reads in 1.01s.
（中略）
Loaded 314483 reads in 0.53s.
Mapped in 0.56s.
No more reads.
Mapped in 0.18s.
Mapped all reads in 30.17s.
Number of reads: 21814483.
Number of mapped reads: 12508635.
Number of uniquely mapped reads: 11368866.
Number of reads have multi-mappings: 1139769.
Number of candidates: 227477151.
Number of mappings: 12508635.
Number of uni-mappings: 11368866.
Number of multi-mappings: 1139769.
Sorted, deduped and outputed mappings in 8.99s.
# uni-mappings: 10833419, # multi-mappings: 1020796, total: 11854215.
Number of output mappings (passed filters): 9211499
Total time: 86.82s.

real    1m32.913s
user    4m12.370s
sys     0m27.555s

• --preset chip オプションでChIP-seqに合わせたパラメータセットになります。 このオプションでは出力はBED形式になります。
• 途中のログを見ると、冗長なリードは除去され、multi readを除いた unique readsのみが出力されているようです。
• マッピング終了後にスタッツが表示されていますね（"Number of ..."の行）。入力リード数は21814483、マップされたリードは12508635、unique readsは11368866となっています。
• Number of output mappings (passed filters): 9211499 となっていますが、調べた感じではpassed filtersはクオリティ（q-value）が30を超えたリードのことを指すようです。最終的に出力されるリード数がこれになります。

chromapでは1m32sという実行時間になりました。１回しか計測していないのでざっくりですが、だいぶ速いと言えそうです。

比較のため、bowtie2とBWAでもマッピングしてみましょう。index作成は省略します。

$ time bowtie2 -x bowtie2-UCSC-hg38 -p 24 SRR227359.fastq.gz > SRR227359.bowtie2.sam
21814483 reads; of these:
  21814483 (100.00%) were unpaired; of these:
    8003629 (36.69%) aligned 0 times
    9880708 (45.29%) aligned exactly 1 time
    3930146 (18.02%) aligned >1 times
63.31% overall alignment rate

real    1m50.454s
user    36m56.341s
sys     5m40.437s

$ time bwa mem -t24 bwa-UCSC-hg38 SRR227359.fastq.gz > SRR227359.bwa.sam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 4705884 sequences (240000084 bp)...
[M::process] read 4705884 sequences (240000084 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1092.230 CPU sec, 47.042 real sec
[M::process] read 4705884 sequences (240000084 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1236.997 CPU sec, 51.901 real sec
[M::process] read 4705884 sequences (240000084 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1199.590 CPU sec, 50.164 real sec
[M::process] read 2990947 sequences (152538297 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1202.052 CPU sec, 50.408 real sec
[M::mem_process_seqs] Processed 2990947 reads in 708.905 CPU sec, 29.680 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t24 /work/Database/bwa-indexes/UCSC-hg38 SRR227359.fastq.gz
[main] Real time: 243.498 sec; CPU: 5451.335 sec

real    4m3.885s
user    88m52.498s
sys     1m59.045s

Bowtie2は1m50s, BWAは4m3sになりました。Bowtie2も速いですね。BWAだけちょっと遅いですが、もともとこんなだったでしょうか？

ChIP-seqデータを入力にした場合、既存法と比べてそこまで劇的に高速化しているという感じではありません。 しかし内部的に冗長リードの除去までやってくれていることを考えると（通常はSAMtoolsやPicardを使って除去する必要があるため）、chromapは高速かつ簡便であると言えそうです。

出力をSAMにしてみる

chromap ではChIP-seqデータに対してBED形式の出力が推奨されています。--SAMオプションをつけることでSAM形式で出力可能ですが、低速になるため推奨しないとのことです（マニュアル参照）。

どのくらい遅くなるのかやってみましょう。

$ time chromap --preset chip -t 24 \
$ -x chromap-UCSC-hg38 -r hg38.fa \
$ -1 SRR227359.fastq.gz --SAM -o SRR227359.chromap.sam
Preset parameters for ChIP-seq are used.
Start to map reads.
（中略）
Mapped all reads in 62.61s.
Number of reads: 21814483.
Number of mapped reads: 12508635.
Number of uniquely mapped reads: 11368866.
Number of reads have multi-mappings: 1139769.
Number of candidates: 227477151.
Number of mappings: 12508635.
Number of uni-mappings: 11368866.
Number of multi-mappings: 1139769.
Sorted, deduped and outputed mappings in 31.94s.
# uni-mappings: 10859912, # multi-mappings: 1023343, total: 11883255.
Number of output mappings (passed filters): 9236028
Total time: 109.61s.

real    2m4.987s
user    6m24.888s
sys     0m48.419s

2m4sなので、ちょっと遅くなりましたが、十分許容範囲内かと思います。BEDで出力した時と微妙にスタッツが変化しているのが謎です。

chromapは結果を標準出力に出力することができないため、BAMやCRAMファイルにするためには一旦SAM形式で出力する必要があります。 また、SAMで遅くなる理由は主にI/O速度であると思われるので、リード数が多いほど書き込みのオーバーヘッド時間が増加すると想像されます。

その他のオプション: 冗長リードの除去をオフにしたい場合

--preset chip オプションに--remove-pcr-duplicatesオプションが含まれていますので、冗長リードのフィルタをオフにしたい場合は --preset chip オプションを使わないようにします。

$ time chromap -l 2000 --low-mem --BED -t 24 \
$ -x chromap-UCSC-hg38 -r hg38.fa \
$ -1 SRR227359.fastq.gz --SAM -o SRR227359.chromap.include_dup.sam
Start to map reads.
Parameters: error threshold: 8, min-num-seeds: 2, max-seed-frequency: 500,1000, max-num-best-mappings: 1, max-insert-size: 2000, MAPQ-threshold: 30, min-read-length: 30, bc-error-threshold: 1, bc-probability-threshold: 0.90
Number of threads: 24
Analyze bulk data.
Won't try to remove adapters on 3'.
Won't remove PCR duplicates after mapping.
Will remove PCR duplicates at bulk level.
Won't allocate multi-mappings after mapping.
Only output unique mappings after mapping.
Only output mappings of which barcodes are in whitelist.
（中略）
Mapped all reads in 74.75s.
Number of reads: 21814483.
Number of mapped reads: 12508635.
Number of uniquely mapped reads: 11368866.
Number of reads have multi-mappings: 1139769.
Number of candidates: 227477151.
Number of mappings: 12508635.
Number of uni-mappings: 11368866.
Number of multi-mappings: 1139769.
Sorted and outputed mappings in 42.73s.
# uni-mappings: 11368866, # multi-mappings: 1139769, total: 12508635.
Number of output mappings (passed filters): 9685899
Total time: 144.08s.

real    2m29.591s
user    8m41.938s
sys     1m0.581s

その他のオプション: クオリティによるフィルタリングをオフにしたい場合

chromapはリードのクオリティで自動的にフィルタが入るようになっていますが、（そういう需要はあまりないと思いますが）フィルタをオフにしたい場合は-q 0をつければOKです。

$ time chromap -l 2000 --low-mem -q 0 --BED -t 24 \
$ -x chromap-UCSC-hg38 -r hg38.fa \
$ -1 SRR227359.fastq.gz --SAM -o SRR227359.chromap.include_dup.q0.sam

Mapped all reads in 78.28s.
Number of reads: 21814483.
Number of mapped reads: 12508635.
Number of uniquely mapped reads: 11368866.
Number of reads have multi-mappings: 1139769.
Number of candidates: 227477151.
Number of mappings: 12508635.
Number of uni-mappings: 11368866.
Number of multi-mappings: 1139769.
Sorted and outputed mappings in 47.17s.
# uni-mappings: 11368866, # multi-mappings: 1139769, total: 12508635.
Number of output mappings (passed filters): 12508635
Total time: 188.37s.

real    3m11.691s
user    8m59.719s
sys     1m8.940s

マップデータの可視化（100kbp bin）

chromapの論文では、ピーク領域についてはbowtie2とほぼ同一の結果が得られることが示されています。 ここでは得られたマップデータのリード分布が同一かどうかをおおまかに確かめます。 染色体レベルでの大まかなリード分布を調べるため、100kbp binでサンプル間比較をしてみましょう。

データをbigWigに変換するためにDROMPAplus*3に含まれるparse2wig+を使います。 bowtie2とBWAの出力ファイルには冗長リードが含まれていますが、parse2wig+は内部で冗長リードを除去するため、得られるBigWigはどのサンプルでも冗長リードは除去されています。 parse2wig+の使い方詳細はこちらを参照してください。

gt=genome_table.hg38.txt
bin=100000 # 100kbp
parse2wig+ -i SRR227359.chromap.bed -o SRR227359.chromap.bed --gt $gt --binsize $bin -n GR -f TAGALIGN
parse2wig+ -i SRR227359.chromap.sam -o SRR227359.chromap.sam --gt $gt --binsize $bin -n GR
parse2wig+ -i SRR227359.bowtie2.sam -o SRR227359.bowtie2.sam --gt $gt --binsize $bin -n GR
parse2wig+ -i SRR227359.chromap.include_dup.sam -o SRR227359.chromap.include_dup.sam --gt $gt --binsize $bin -n GR
parse2wig+ -i SRR227359.bwa.sam -o SRR227359.bwa.sam --gt $gt --binsize $bin -n GR

1行目ではBED形式のマップデータを読み込むために-f TAGALIGNオプションをつけています。 マップデータのクオリティスタッツにはほとんど違いはありませんでしたが、chromapが最もマップ率が低い、すなわち厳しい閾値を用いている感じでした（結果は割愛）。 なおBowtie2, BWAの出力ファイルにはunmappedも含めた全リードが含まれているのに対し、chromapの出力ファイルにはmapped readしか含まれていないようです。

DROMPA+のGVコマンドを使って可視化します。結果はGV.pdfに保存されます。

gt=genome_table.hg38.txt
ideogram=DROMPAplus/data/ideogram/hg38.tsv
drompa+ GV --gt $gt \
 -i parse2wigdir+/SRR227359.chromap.bed.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,chromap.bed \
 -i parse2wigdir+/SRR227359.chromap.sam.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,chromap.sam \
 -i parse2wigdir+/SRR227359.chromap.include_dup.sam.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,chromap.include_dup.sam \
 -i parse2wigdir+/SRR227359.bwa.sam.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,bwa.sam \
 -o GV --showctag 1 --showitag 1 --ideogram $ideogram

結果は以下のようになりました。

chromap.GV.chr1

chromap.GV.chrX

緑のラインがChIP, 青のラインがInput、赤で示されているラインがChIP/Input ratioです。 ここではbowtie2のマップデータをInputにし、左に示されている４つのサンプルをChIPとしてratioを計算していますので、ratioはbowtie2のマップデータに対する増減を示しています。

結果を見ると、染色体腕部ではCNVによるリード数のばらつきが認められますが、ratioではきちんとキャンセルされており、ツール間でほぼ同数のリード数が得られていることがわかります。 一方セントロメア周辺部では、chromapのマップリード数が極端に少ないという結果になりました。 chromapとparse2wig+で冗長リードの定義が異なるせいかと思い、chromapで冗長リードを除去しないオプションで得たマップファイルでも試してみましたが（３段目）、やはりセントロメア周辺のリード数が少なくなっています。 セントロメア周辺はhighly repetitive regionsなので、chromapのキャッシュを用いたアラインメントはリピート領域に弱いか、あるいはchromapが内部的にリピート領域に対して厳しいスコアリングをしていると思われます。通常のChIP-seq解析ではリピート領域は対象にしませんのでどのツールを使っても変わらないと思いますが、H3K9me3のようにセントロメア周辺に特に濃縮するようなサンプルの場合は検討が必要そうです。なお、BWAは逆にリピート領域により多量に張り付くという結果になりました。

（2022/3/23追記） ひょっとしてクオリティフィルタでリピート領域が除外されているのかと思い、-q 0オプションを追加して得られたマップデータについても同様に可視化してみましたが、若干リード量が回復したものの同じにはならず、やはりchromapではリピート領域のリード数は少ないままでした。やはりlight-weight alignmentはリピートに弱いのかもしれません。

まとめ

chromapについて、ChIP-seqデータを用いて簡単に試してみました。出力形式をSAMにしてもそこまで遅くはならないので、SAM出力からBAM/CRAMに変換するようにしておくと、他のツールとの互換性が高まると思います。

次回はHi-Cデータに対してchromapを試してみたいと思います。

rnakato.hatenablog.jp

*1:H. Zhang et al., Fast alignment and preprocessing of chromatin profiles with Chromap, Nature Communications, 2021, DOI 10.1038/s41467-021-26865-w

*2:H. Li, Minimap2: pairwise alignment for nucleotide sequences, Bioinformatics, 2018, DOI 10.1093/bioinformatics/bty191

*3:R. Nakato, T. Sakata, Methods for ChIP-seq analysis: A practical workflow and advanced applications, Methods, 2020, DOI 10.1016/j.ymeth.2020.03.005

Docker/Singularityでも環境の統一は難しいかもしれないという話

Docker, Singularityが何かということを以下の記事で過去に紹介しました。

rnakato.hatenablog.jp

rnakato.hatenablog.jp

その後実験医学増刊号でもDockerを用いた私のシングルセル解析プラットフォームを紹介させていただきました。ここでも「Dockerは環境が統一できるので有用だよ！」と宣伝しているのですが、実際には、不特定多数の人に使ってもらう環境だと思った以上に環境統一は難しいのでは…？という気がしてきました。以下、その理由です。

・ 皆新しく自分の環境にパッケージをインストールする

「できるから」ということで、当然ではあるのですが。Dockerは実行したコンテナ上で新たにパッケージをインストール可能ですし、SingularityにもSandboxという機能があり、起動したSingularityイメージ上でpipなどでツールをインストールすると自動的にホームディレクトリ内にキャッシュを作成し、ユーザオリジナルの環境をビルドするようです。自分でパッケージをインストール可能であることはコンテナの利点の一つではあるのですが、そのためユーザごとに異なる環境になり得るので冪等性は失われます。かつ、ユーザは意外と自分が何かインストールしたという事実を覚えていません…。「こんなエラーが出るのですが」と質問をいただき、環境が再現しないのでなぜだろうか…とよくよく聞いていると何かしらパッケージをインストールしていてバージョンがコンフリクトしている、のようなこともありました。今後エラー質問メールを受けるたびにそれを確認しないといけないのか？と思うとなかなか大変です。

・ マウントしたローカルディレクトリのツール・ライブラリを読み込む可能性がある

ローカルにあるデータを読み込むためにディレクトリをマウントすることがありますが、マウントされたディレクトリ内にバイナリツールが含まれており、そこにパスが通っていたりすると、イメージ内のツールではなくローカルのツールが呼び出される可能性があります。例えばRのライブラリは一般ユーザ環境ではホームディレクトリに ~/R ディレクトリが作成され、そこに格納されますが、ホームディレクトリ上でSingularityを起動するとホームディレクトリが自動的にマウントされますから、パスによってはローカルのRライブラリが呼び出される可能性もあるわけです。

（Singularityはどこで起動してもホームディレクトリをマウントしているっぽいですね。そうなると.jupyterなどの設定はローカルPCのものが反映されることになりますね。）

また、ツールによっては内部で別のツールを呼び出すような設計になっているものがあります。例えばscATAC-seqのためのSignacはピークコールのために内部でMACS2を呼び出しますが、以下のように macs2.path オプションによって起動するmacs2バイナリを指定可能です。ここで例えばローカルのMACS2を呼び出す設定にしているとすると、それはもう同一環境とは呼べません。これもまた、外部からは見えにくいバグの要因になり得ます。

peaks <- CallPeaks(
  object = pbmc,
  group.by = "predicted.id",
  macs2.path = "/usr/local/bin/macs2"
)

・ そもそもインストールしたツールにバグ入ってたりしがち

これはイメージを作成する側の問題でユーザ依存ではないのですが、上で紹介したシングルセル解析プラットフォームに含まれているツール群は個人ベースで開発しているものが多い＋大量の依存関係を持つツールが多いので、アップデートの際にライブラリバージョン依存のバグが入ることがしばしば起きます。私のDockerfileではバージョンを指定せずにツールをインストールしているので（怠慢？すみません）、Dockerfileが同一でも、イメージ作成時によって含まれるツール群のバージョンが異なり、そのためエラーが起きる可能性があります。実際にそういうことは何度も起きています…なので同じDockerfileでもうかつにビルドをやり直さない方がいいのだなという教訓は得ました。

結論

「大量のツール群をインストールした環境」のイメージ、というのがそもそもコンセプト的に矛盾しているような気もしてきました。virtual environmentを切り分けるように、用途ごとに別々のイメージ作るのが本来の使い方だという気もするのですが、自分用ではなく不特定多数の人に公開する目的では切り分けしづらいという感覚もあります。Dockerfileをきれいにしてくれる校正業者とかいたらいいんですけど f^_^;