完全長ゲノム配列を使って大規模SVsを反映したvariant callingをやってみる

動機

バクテリアやウイルスゲノムのバリアントコールといえばSnippyがポピュラー。しかしショートリードベースであり (ゲノムアセンブリを入力するとショートリードを生成する)、大規模な欠失や挿入などの構造変異 (structural variants; SV) を考慮できない。
これは、完全長の全ゲノム配列の情報をフルに活用できないことを意味する。
なんかもったいない気がする。

そこで、全ゲノムアラインメントから直接バリアントコールすることで、大規模なSVを反映したvariant calling format (VCF) ファイルの生成を試みた。
対象種は全長約300 kb、塩基配列一致率99%で数kb-数十kb単位のSVが存在する完全長ウイルスゲノム (WSSV) 。

Reference: NC_003225.3 White spot syndrome virus strain CN01, complete genome
Query: MH663976.1 Procambarus clarkii virus isolate 1601, complete genome

Requirements

• minimap2
• bcftools
• bgzip
• SVIM-asm

概要

Minimap2でゲノム間のアラインメントを作成し、bcftools mpileupでバリアントコールする (ここはSnippyでもいいかもしれない)。
得られたVCFファイルはSNPsや小さなindelを網羅するが、キロベース単位の欠失など大きなSVsはカバーできない。
そこで、SVIM-asmを用いてより大きなSVsのバリアントコールを行う。大規模SVsのコールに優れるSVIM-asmは、SNPs等の漏れが多い。
2つのバリアントコーラーを併用することでできる限りコール漏れをなくす。
最後に二つのVCFをマージして完成。

REF=CN01
QUERY=1601

# Map query sequence against the reference by Minimap2
minimap2 -k 27 -ax asm5 --cs -r2k,40k -z 500000 -t 8 ${REF}.fa ${QUERY}.fa|samtools sort  > ${QUERY}.bam
samtools index ${QUERY}.bam

# Call variants by bcftools mpileup
bcftools mpileup  -Ov -f ${REF}.fa ${QUERY}.bam | bcftools call --ploidy 1 -P 1 -mv -Ov -o ${QUERY}.mpileup.tmp.vcf

# Make some necessary modifications
cat ${QUERY}.mpileup.tmp.vcf |sed "s/GT:PL/GT/g" |sed "s@1:60,0@1/1@g" >${QUERY}.mpileup.vcf 

# Compress and index the VCF files
bgzip -f ${QUERY}.mpileup.vcf
bcftools index -f ${QUERY}.mpileup.vcf.gz

# Run SVIM-asm
conda activate svim-asm-1.0.3
svim-asm haploid ${QUERY}_svimasm --min_sv_size 1 --max_sv_size 50000 ${QUERY}.bam ${REF}.fa --sample ${QUERY}.bam
conda deactivate

# Compress and index the VCF files
bgzip -f ${QUERY}_svimasm/variants.vcf 
bcftools index -f ${QUERY}_svimasm/variants.vcf.gz

# Concatenate, sort and normalize variants
bcftools concat --ligate-force ${QUERY}.mpileup.vcf.gz  ${QUERY}_svimasm/variants.vcf.gz |bcftools sort |bcftools norm   -f ${REF}.fa -m +any -Oz -o ${QUERY}.concat.norm.vcf.gz

# Index the VCF files
bcftools index -f ${QUERY}.concat.norm.vcf.gz

# Verify saved variants by applying the variants onto the reference sequence and get the consensus sequence
cat ${REF}.fa | bcftools consensus ${QUERY}.concat.norm.vcf.gz > ${QUERY}.consensus.fa

結果

以下のような結果となった。

上から

• mpileup
• SVIM-asm
• 完成形 (mpileup + SVIM-asm)
  となる。
  下は一部拡大。
  
  mpileupとSVIM-asmが互いを補うことで、SNPsから大規模な欠失までを網羅したVCFファイルができた。
  今回の例はbcftools consensusでクエリ配列を完全再現することができた。
  しかし他の細菌・ウイルスゲノムではところどころ抜けがあり、完全再現はできていない。
  minimap2のパラメータ (-k 27 -ax asm5など) を調節することでうまくいくかもしれない。
  そもそもパンゲノム解析が当たり前になりつつある現在VCFファイルを作ること自体時代遅れなのかもしれない。