Bismarkには大変お世話になりました
N年前にメチロームの解析を始めて以来ずっとbismark一筋でしたが1サンプル半日とかかかるのにそろそろうんざりしてきました。
Methylpyがよいらしい
Methylpy:
Matthew D. Schultz, Yupeng He, John W.Whitaker, Manoj Hariharan, Eran A. Mukamel, Danny Leung, Nisha Rajagopal, Joseph R. Nery, Mark A. Urich, Huaming Chen, Shin Lin, Yiing Lin, Bing Ren, Terrence J. Sejnowski, Wei Wang, Joseph R. Ecker. Human Body Epigenome Maps Reveal Noncanonical DNA Methylation Variation. Nature. 523(7559):212-216, 2015 Jul.
Ecker labから出てるんだね。
river winさんもこれが良いと言っているという風の噂を聞いたので導入してみることにした。
Install
pipでいけるがconda派なのでcondaでやってみた
conda install -c bioconda methylpy
とくに問題なくインストール完了
ゲノムのBuild
methylpy build-reference \
--input-files hoge.fasta \
--output-prefix hoge
チュートリアルにある--bowtie2 Trueというオプションはバージョンが変わって無くなりました。
トリミング->マッピング->メチル化カウント
FASTQ=`find *1.fastq.gz`
for int in ${FASTQ[@]}; do
fastq=`echo "$int" | sed s/_R1.fastq.gz//g`
# trimming
trimmomatic PE -threads 20 \
${fastq}_R1.fastq.gz \
${fastq}_R2.fastq.gz \
${fastq}_paired_R1.fastq.gz \
${fastq}_unpaired_R1.fastq.gz \
${fastq}_paired_R2.fastq.gz \
${fastq}_unpaired_R2.fastq.gz \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:50 HEADCROP:10
ref_path=/path/to/ref_fasta/
#Mapping
mkdir ./${fastq}_output
methylpy paired-end-pipeline \
--read1-files ${fastq}_paired_R1.fastq.gz \
--read2-files ${fastq}_paired_R2.fastq.gz \
--sample ${fastq}_BS-seq \
--forward-ref ${ref_path}hoge_f \
--reverse-ref ${ref_path}hoge_r \
--ref-fasta ${ref_path}hoge.fasta \
--num-procs 20 \
--trim-reads FALSE \
--path-to-output ./${fastq}_output
done
フォルダにある全てのfastqをそれぞれtrimmomaticでトリミングし、メチル化カウントまで一気に処理
個人的な趣味でトリミングにはtrimmomaticを使用。それに伴いmethylpy paired-end-pipelineでは--trim-readsはFALSE指定にしている。
最終産物としてall-c_hoge.tsv.gzというファイルが得られる。bismarkでいうところのcx_reportみたいなやつ
数千万リードがだいたい3時間弱で終わった。速い！
Wigファイルの作成
methylpy allc-to-bigwig \
--allc-file allc_hoge.tsv.gz \
--output-file methylpy_hoge_CG \
--ref-fasta /path/to/ref_fasta/hoge.fasta \
--mc-type CGN \
--bin-size 50
methylpy allc-to-bigwig \
--allc-file allc_hoge.tsv.gz \
--output-file methylpy_hoge_CHG \
--ref-fasta /path/to/ref_fasta/hoge.fasta \
--mc-type CHG \
--bin-size 50
methylpy allc-to-bigwig \
--allc-file allc_hoge.tsv.gz \
--output-file methylpy_hoge_CHH \
--ref-fasta /path/to/ref_fasta/hoge.fasta \
--mc-type CHH \
--bin-size 50
allc-file: でき上がったall-c_hoge.tsv.gz
--ref-fasta :リファレンスfasta
--mc-type :wigにしたいコンテクストを指定。例えば植物だとCGN, CHG, CHHを三つを上のようにそれぞれ指定してwigファイルを作成する
なにせ速い。もうこっちに切り替えよう。
追記
allc-to-bigwigを走らせてもwigファイルしかできないのでなんでだろうと思っていたが
wigToBigWig が正常に走っていないのが原因だった。
libsslをうまく参照できていないのが動かない原因のよう。
/miniconda3/libに移動して以下のようにシムリンクを貼り直す。
(base) [hoge lib]$ ln -s libssl.so.1.1 libssl.so.1.0.0
(base) [hoge lib]$ ln -s libcrypto.so.1.1 libcrypto.so.1.0.0
そうすると正常に動いた。似たような報告があったのでそれを参考した。
https://qiita.com/kohei-108/items/9129ad953b0906e115a6