テロメア単位の完全なゲノムアセンブリは、解析の向上や新しいバリアントの発見を期待できるが、多くの重要なゲノムリソースは古いリファレンスゲノムと関連したままである。そのため、リファレンスゲノム間のゲノムフイーチャーやリードアラインメントをトランスレートする必要がある。本論文では、levioSAM2と呼ばれる新しい手法を紹介する。この手法は、リファレンスの変更を考慮し、全ゲノムマップを用いたアセンブリ間のリフトオーバーを高速かつ正確に実行する。複数のリファレンスを使用できることに加え、高品質のリファレンス（T2T-CHM13など）にリードをアライメントし、古いリファレンス（GRCh38など）にリフトオーバーすると、古いリファレンスでの結果のバリアントコールの精度が実際に向上することを実証した。T2T-CHM13の品質向上を活用することにより、levioSAM2は、実際のイルミナデータセットを用いたGRCベースのマッピングと比較して、スモールバリアントコールエラーを11.4～39.5％減少させるる。LevioSAM2は、ロングリードベースの構造バリアントコールも改善し、PacBio HiFiデータセットに対して3.8～11.8%のエラーを減少させた。特に、GRCリファレンスが低品質である複雑な医療関連遺伝子のセットで性能が向上している。このソフトウェアは、MITライセンスのもと、https://github.com/milkschen/leviosam2 で入手できる。
T2T references give drastic assembly improvements, but older ones are useful with rich annotations. LevioSAM2 is an alignment method aiming for the best of both worlds, exploiting the best assembly while reporting w/r/t an older assembly (like GRC) 1/ https://t.co/zNhwU3PiqS
— Nae-Chyun Chen (@naechyun_chen) April 28, 2022
Small- and structural-variant calling for long reads improves as well. We observed large-scale mapping artifacts resolved by levioSAM2, e.g., parts of genes KMT2C and MAP2K3. Top: GRCh37; bottom: levioSAM2. 5/ pic.twitter.com/gaPTQ2YxSW
— Nae-Chyun Chen (@naechyun_chen) April 28, 2022
workflow（小さなデータを使ったcase studyもあり）
https://github.com/milkschen/leviosam2/blob/main/workflow/README.md
インストール
ubuntu18でcondaで環境を作ってテストした。
#conda(link)
mamba create -n leviosam2 -y
conda activate leviosam2
mamba install -c conda-forge -c bioconda leviosam2 -y
#docker
docker pull naechyun/leviosam2:latest
#singularity
singularity pull docker://naechyun/leviosam2:latest
> leviosam2 -h
Program: leviosam2 (lifting over alignments)
Version: 0.2.1
Usage: leviosam2 <command> [options]
Commands: index Index a lift-over map.
lift Lift alignments in the SAM/BAM formats.
bed Lift intervals in the BED format.
collate Collate lifted paired-end alignments.
reconcile Reconcile alignments.
Options: -h Print detailed usage.
-V Verbose level [0].
Index a lift-over map using either a chain file.
Usage: leviosam2 index -c <chain> -p <out_prefix> -F <fai>
Inputs: -c string Index a lift-over map from a chain file.
-F string Path to the FAI (FASTA index) file of the dest reference.
-p string The prefix of the output file.
Perform efficient lift-over using leviosam2.
Usage: leviosam2 lift [options] -C <clft>
Inputs: -C string Path to an indexed ChainMap.
Options: -a string Path to the SAM/BAM file to be lifted.
Leave empty or set to "-" to read from stdin.
-t INT Number of threads used. [1]
-T INT Chunk size for each thread. [256]
Each thread queries <-T> reads, lifts, and writes.
Setting a higher <-T> uses slightly more memory but might benefit thread scaling.
-m add MD and NM to output alignment records (requires -f option)
-f string Fasta reference that corresponds to input SAM/BAM (for use w/ -m option)
-x string Alignment preset []
-G INT Number of allowed CIGAR changes for one alingment. [0]
Commit/defer rule options:
-S string<:int/float> Key-value pair of a split rule. We allow appending multiple `-S` options.
Options: mapq:<int>, aln_score:<int>, isize:<int>, hdist:<int>, clipped_frac:<float>, lifted. [none]
* mapq INT Min MAPQ to commit (pre-liftover). [30]
* aln_score INT Min AS:i (alignment score) to commit (pre-liftover). [100]
* isize INT Max TLEN/isize to commit (post-liftover). [1000]
* hdist INT Max NM:i (Hamming dist.) to commit (post-liftover). `-m` and `-f` must be set. [5]
* clipped_frac FLOAT Min fraction of clipped to commit (post-liftover). [0.95]
Example: `-S mapq:20 -S aln_score:20` commits MQ>=20 and AS>=20 alignments.
-r string Path to a BED file (source coordinates). Reads overlap with the regions are always committed. [none]
-D string Path to a BED file (dest coordinates). Reads overlap with the regions are always deferred. [none]
実行方法
１、leviosam2 indexの実行。LevioSAM2は、リフトオーバークエリーのためにchainファイル（UCSCのhg38 Chain Files）からインデックスを作成する。
leviosam2 index -c source_to_target.chain -p target -F target.fai
- -c string Index a lift-over map from a chain file.
- -F string Path to the FAI (FASTA index) file of the dest reference.
- -p string The prefix of the output file.
target .clft拡張子のインデックスファイルが出力される。
２、leviosam2 liftの実行。インデックスファイルと入力のbamファイルを指定する。bamファイルはcoordinate sortされている必要がある。
leviosam2 lift -C source_to_target.clft -a source.bam -p lifted_from_source -O bam
- -C string Path to an indexed ChainMap.
- -a string Path to the SAM/BAM file to be lifted.
- -t Number of threads used. [1]
lifted_from_source.bamが出力される。
チュートリアルの実行結果
引用
Improved sequence mapping using a complete reference genome and lift-over
Nae-Chyun Chen, Luis F Paulin, Fritz J Sedlazeck, Sergey Koren, Adam M Phillippy, Ben Langmead
bioRxiv, Posted April 28, 2022
LevioSAM: fast lift-over of variant-aware reference alignments
Taher Mun, Nae-Chyun Chen, Ben Langmead Author Notes
Bioinformatics, Volume 37, Issue 22, 15 November 2021, Pages 4243–4245
関連
