Nanopore データ、Megalodon から phasing まで
1. ONT、basecalling & mapping
＜環境 GPU ノード＞
Intel Xeon Gold 6326 CPU 2.90GHz x 2, 256GB memory
NVIDIA GA100 [A100 PCIe 40GB] x 4
Rocky Linux release 8.5
megalodon v2.5.0 (https://github.com/nanoporetech/megalodon)
Remora v2.0.0
guppy v6.4.2
＜データ＞
R9.4.1 フローセル、PromethION 出力の fast5 データ。
$ myref=GRCh38.mmi
$ nproc=16
$ gpu_id="0 1"
$ megalodon ${fast5_pass}
--guppy-config dna_r9.4.1_450bps_sup_prom.cfg
--remora-modified-bases dna_r9.4.1_e8 sup 0.0.0 5mc CG 0
--outputs basecalls mappings mod_mappings mods
--reference ${myref}
--output-directory ${output_dir}
--processes ${nproc}
--devices ${gpu_id}
--guppy-server-path path_to_guppy_basecall_serverguppy-config はguppyインストールパス/data/にある。
2023.3 現在、dna_r9.4.1_450bpsが最新。
r9.4.1はフローセルバージョン。
fast/hac/sup がある（速度⇄ 精度？）。
prom は promethION。
--remora-modified-bases の選択肢は
$ remora model list_pretrained
で確認できる。
--reference は minimap2 で使用する、マッピングするゲノム。
ゲノムのfastaファイルもしくはminimap2のインデックスファイルを指定できるが、minimap2のインデックスファイルだとメモリ使用量が改善されるらしい。こちらの環境ではR9.4.1 フローセル1ランでおよそ２４時間程度かかっていたものが、minimap2 インデックスへの変更でGPU使用効率が向上し、1/3程度に短縮された。
--process は指定した数値よりCPU使用量が多くなるので調整する必要あり。こちらの環境で私はマシンスペックの半分程度にしている。
--devices はGPUのID。GridEngine で割り振られる場合はその番号に準ずる。
2. SNP calling（リファレンスと対象ゲノム間のSNP）
Clair3でSNP callを行う。
＜環境＞
CentOS Linux 7.6
clair3 v0.1.12 (https://github.com/HKU-BAL/Clair3)
whatshap v1.4 (https://github.com/whatshap/whatshap)
anaconda で clair3 + whatshapをインストール。
$ conda create -n clair3 -c bioconda clair3 -y当時こちらの環境では
AttributeError: module 'numpy' has no attribute 'int'.
というエラーが出たため、エラーが出ない組み合わせを試行錯誤してインストールした。
$ conda create -n clair3-01r12-nmp1235 -c conda-forge -c bioconda clair3 python=3.9.0 numpy=1.23.5 -y$ conda activate clair3
$ bam=mapping.bam
$ ref=GRCh38.fa
$ MODEL_NAME="r941_prom_sup_g5014"
$ THREADS=24
$ run_clair3.sh
--bam_fn=${bam}
--ref_fn=${myref}
--threads=${THREADS}
--platform="ont"
--model_path="path_to_model_dir/${MODEL_NAME}"
この時点でphasingデータは得られているが、いくつかの理由で、今回得られたSNPから、リファレンスゲノムと 対象ゲノムとのSNPから対象ゲノムのリファレンスを作成して、マッピングからやり直した。
3. 対象ゲノムリファレンスの作成
bcftools v1.15.1
bcftoolsインストール
conda create -n bcftools -c bioconda bcftoolsclair3の出力 merge_output.vcf.gz から、REFとALTの長さが１塩基かつ、 10列目の最初のマークが、
A) 1/1 （EF/ALT が ALT1/ALT1）になっているものはREFをALTに、
B) 1/2 （EF/ALT が ALT1/ALT2）になっているものはALTがコンマで2つ書かれているので両方とも長さが１のものについて1つ目をREFに、
変更して出力。出力されものは、リファレンスゲノムと対象ゲノムの違うSNVをREFとして記したVCFとなる。
bgzipで圧縮、tabix も行う。
# merge_output.vcf.gzの一部
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE
chr1 10061 . T C 5.26 PASS F GT:GQ:DP:AF 0/1:5:29:0.4828
chr1 10108 . C CT 0.00 LowQual F GT:GQ:DP:AF 1/1:0:29:0.3103
chr1 10150 . CT C 0.22 LowQual F GT:GQ:DP:AF 0/1:0:29:0.2759
chr1 10177 . A AC 2.34 PASS F GT:GQ:DP:AF 0/1:2:29:0.2414
chr1 10222 . T C 5.74 PASS F GT:GQ:DP:AF 0/1:5:29:0.3793$ gzip -dc ./merge_output.vcf.gz | perl -F'\t' -anle 'if ($=~/^#/) {print $;} elsif ($F[6] eq "PASS" and length($F[3])==1) {if ($F[9] =~/1/1/ and length($F[4])==1) {print $_;} elsif ($F[9]=~/1/2/) {@a=split/,/,$F[4];if (length($a[0])==1 and length($a[1])==1) {$F[9]=~s@1/2@0/1@;print join("\t",@F[0..3],$a[0],@F[5..9]);}}}' > GRCh38_line1.vcf
$ bgzip GRCh38_line1.vcf
$ tabix -p vcf GRCh38_line1.vcf.gz$ conda activate bcftools
$ myref=GRCh38.fa
$ out=GRCh38_line1.fa
$ vcf=GRCh38_line1.vcf.gz
$ bcftools consensus -f ${myref} ${vcf} > ${out}作成した対象ゲノムリファレンスをもとにもう一度megalodonを行う。
# minimap2の index 作成。
$ minimap2 -x map-ont -d GRCh38_line1.mmi GRCh38_line1.fa
$ myref=GRCh38_line1.mmi
$ nproc=16
$ gpu_id="0 1"
$ megalodon ${fast5_pass}
--guppy-config dna_r9.4.1_450bps_sup_prom.cfg
--remora-modified-bases dna_r9.4.1_e8 sup 0.0.0 5mc CG 0
--outputs basecalls mappings mod_mappings mods
--reference ${myref}
--output-directory ${output_dir}
--processes ${nproc}
--devices ${gpu_id}
--guppy-server-path path_to_guppy_basecall_server4. SNP calling （アリル間のSNP）
clair3でSNP callを行う。
$ conda activate clair3
$ bam=mapping.bam
$ ref=GRCh38_line1.fa
$ MODEL_NAME="r941_prom_sup_g5014"
$ THREADS=24
$ run_clair3.sh
--bam_fn=${bam}
--ref_fn=${myref}
--threads=${THREADS}
--platform="ont"
--model_path="path_to_model_dir/${MODEL_NAME}"whatshapでphasingされたVCFファイルが tmp/phase_output/phase_vcf に入っているので、これを合体し、bgzip圧縮。
$ cd clair3_output_dir/tmp/phase_output/phase_vcf
$ gzip -dc phased_chr1.vcf.gz phased_chr2.vcf.gz phased_chr3.vcf.gz phased_chr4.vcf.gz phased_chr5.vcf.gz phased_chr6.vcf.gz phased_chr7.vcf.gz phased_chr8.vcf.gz phased_chr9.vcf.gz phased_chr10.vcf.gz phased_chr11.vcf.gz phased_chr12.vcf.gz phased_chr13.vcf.gz phased_chr14.vcf.gz phased_chr15.vcf.gz phased_chr16.vcf.gz phased_chr17.vcf.gz phased_chr18.vcf.gz phased_chr19.vcf.gz phased_chr20.vcf.gz phased_chr21.vcf.gz phased_chr22.vcf.gz phased_chrX.vcf.gz > ${dir}/GRCh38_line1_phased.vcf
$ conda activate tabix
$ bgzip /${dir}/GRCh38_line1_phased.vcf
$ tabix -p vcf /${dir}/GRCh38_line1_phased.vcf.gz4. ONTのBAMをアリル分け
whstshap v1.14
$ conda create -n whatshap14 -c bioconda whatshap=1.4まずhaplotag情報を作成する。
$ samtools sort -o mappings_sort.bam mappings.bam
$ samtools index mappings_sort.bam
$ conda activate whatshap14
$ in=mappings_sort.bam
$ out=mappings_haplotag.bam
$ vcf=GRCh38_line1_phased.vcf.gz
$ myref=GRCh38_line1.fa
$ haplist=GRCh38_line1_phased_haplotypes.tsv
$ cpu=24
$ whatshap haplotag
-o ${out}
--reference ${myref}
--output-threads=${cpu}
--output-haplotag-list=${haplist}
--ignore-read-groups
${vcf} ${in}
# phasing blocksの確認。
$ vcf=GRCh38_line1_phased.vcf
$ whatshap stats --gtf=${vcf}.gtf ${vcf} > ${vcf}_stats.txthaplotag情報をもとに、BAMを分ける。
$ haplist=GRCh38_line1_phased_haplotypes.tsv
$ whatshap split --output-h1 mappings_H1.bam --output-h2 mappings_H2.bam mappings.bam ${haplist}
$ whatshap split --output-h1 mod_mappings_H1.bam --output-h2 mod_mappings_H2.bam mod_mappings.bam ${haplist}WGBS、EM-seq、RNA-seqなどをphasingするには、SNP calling （アリル間のSNP）で作成したVCFファイルを用いる。
SNPsplit (https://www.bioinformatics.babraham.ac.uk/projects/SNPsplit/)を用いるには、SNP座標をすべてマスクしてマッピングする必要があるので、bedtools maskfastaを用いてマスクしたリファレンスfastaを作成し、これですべて再マッピングを行う。
$ conda activate bedtools
$ in=GRCh38_line1.fa
$ out=GRCh38_line1_mask.fa
$ vcf=GRCh38_line1_phased.vcf.gz
$ bedtools maskfasta -fi ${in} -fo ${out} -bed ${vcf}6. Short read sequencing データをSNP情報をもとに分ける。
SNPsplit v0.3.2
$ conda create -n snpsplit -c bioconda snpsplit -yVCFファイルをSNPsplit用にREFとALTを入れ替える。
10列目の 0|1 or 1|0 をもとに、列4、5の REF,ALT を入れ替える。REF 側を H1 (genome1)、ALT側をH2(genome2)にする。
9列目の最後の:PSがphasingできたフラグ、10列目の最後の数字がphasing blockのID （phasing blockの最初のSNPの座標）なので、この番号が異なる場所はgenome1/genome2 が逆転するかどうかは判定できない。
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE
chr1 10517 . C T 13.06 PASS P GT:GQ:DP:AF:PS 0|1:13:23:0.3913:10517
chr1 11134 . A G 18.09 PASS P GT:GQ:DP:AF:PS 1|0:18:23:0.3913:10517
chr1 11409 . A G 14.99 PASS P GT:GQ:DP:AF:PS 1|0:14:24:0.3333:10517
chr1 11457 . C G 19.01 PASS P GT:GQ:DP:AF:PS 1|0:19:25:0.4:10517$ gzip -dc GRCh38_line1_phased.vcf.gz | perl -F'\t' -anle '@t=split(/:/,$F[9]);@F[3,4] = @F[4,3] if $t[0] eq "1|0";print join("\t","$F[0]_" . $i++,@F[0,1],1,"$F[3]/$F[4]","$F[0]:$t[4]") if $t[4] =~/^\d+$/;' > GRCh38_line1_phased_snpsplit.txtこれをもとにSNPsplitでBAMをgenome1とgenome2に分ける。
--bisulfiteオプションは、WGBS／EMseqで有効。Transcriptomeなどでは外す。SNP位置をマスクしていないと分けられない。
$ conda activate SNPsplit
$ SNPsplit --snp_file GRCh38_line1_phased_snpsplit.txt --bisulfite sample1.bam2023/03/14