T2Tゲノムアセンブリの評価ツール Genome Continuity Inspector (GCI)

　最近のロングリードシーケンス技術の進歩により、高品質なゲノムアセンブリの作成が大幅に容易になった。テロメアtoテロメアなギャップレス（T2T）アセンブリは、ゲノムアセンブリの新たなゴールデンスタンダードとなっている。最近、T2Tレベルのリファレンスゲノムを作ろうとする試みがいくつかなされている。しかし、ゲノムアセンブリがT2Tレベルであることを証明する普遍的な基準はまだ見つかっていない。従来のゲノムアセンブリー評価指標（N50やその派生指標）では、T2Tに近いアセンブリと真のT2Tアセンブリを、グローバルにもローカルにも連続的に区別することはできなかった。また、これらの評価指標は生のリードに依存しないため、人為的な操作によって簡単に膨れ上がってしまう。したがって、完全ゲノム時代において、真の完全性を反映する一塩基分解能でのギャップレス評価ツールが早急に必要である。

ゲノムアセンブリの連続性を1塩基分解能で評価するGenome Continuity Inspector (GCI)と呼ばれるツールを紹介する。これはゲノムアセンブリがT2Tレベルにどれだけ近いかを評価する。GCIは、複数のプラットフォームからのロングリードをアセンブリにマッピングするために、複数のアライナーを利用する。GCIは、信頼性の高いリードアライメントのキュレーションされたマッピングカバレッジを組み込むことで、潜在的なアセンブリの問題を特定する。また、全ゲノムまたは染色体スケールでのアセンブリ全体の連続性を定量化するためにGCIスコアを報告する。オープンソースのGCIコードは、MITライセンスの下、Github (https://github.com/yeeus/GCI)で自由に利用できる。

レポジトリより

Genome Continuity Inspector (GCI) は、T2Tゲノムなどの高品質ゲノムを対象とした、塩基分解能でのアセンブリ評価ツールである。ロングリード（PacBio HiFiおよび/またはOxford Nanoporeロングリード）をゲノムアセンブリにマッピングすることによって生成されたアラインメントを厳密にフィルタリングした後、GCIはアセンブリの潜在的な問題を報告し、さらにアセンブリの連続性を定量化するスコアを報告する。

インストール

condaで環境を作ってテストした。

依存

canu (for trio-binning)
minimap2 (for mapping)
winnowmap (for mapping)
samtools (for sam/bam processing)
paftools.js (for converting sam to paf)

As for GCI, it requires:

python3.x (tested in python3.10)
pysam (stable version)
biopython (stable version)
numpy (stable version)
matplotlib (stable version)
bamsnap (for plotting in utility filter_bam.py)

Github

mamba create -n CGI python=3.11 -y
conda activate CGI
mamba install -c bioconda winnowmap -y
mamba install -c bioconda samtools -y
mamba install -c bioconda canu -y
mamba install -c bioconda -y minimap2 -y
pip install pysam numpy biopython matplotlib bamsnap 

#本体
git clone https://github.com/yeeus/GCI.git
cd GCI/

> python GCI.py --help

usage: GCI.py [-r FILE] [--hifi [...]] [--nano [...]] [--chrs] [-R FILE] [-ts INT] [-dp FLOAT] [-d PATH] [-o STR] [-mq INT] [--mq-cutoff INT] [-ip FLOAT] [-op FLOAT] [-cp FLOAT] [-fl INT] [-p] [-dmin FLOAT] [-dmax FLOAT] [-ws INT] [-it STR] [-f]

[-h] [-v]

A program for assessing the T2T genome

Input/Output:

-r FILE, --reference FILE

The reference file

--hifi [ ...] PacBio HiFi reads alignment files (at least one bam file)

--nano [ ...] Oxford Nanopore long reads alignment files (at least one bam file)

--chrs A list of chromosomes separated by comma

-R FILE, --regions FILE

Bed file containing regions

Be cautious! If both specify `--chrs` and `--regions`, chromosomes in regions bed file should be included in the chromosomes list

-ts INT, --threshold INT

The threshold of depth to be reported as issues [0]

-dp FLOAT, --dist-percent FLOAT

The distance between the candidate gap intervals for combining in chromosome units [0.005]

-d PATH The directory of output files [.]

-o STR, --output STR Prefix of output files [GCI]

Filter Options:

-mq INT, --map-qual INT

Minium mapping quality for alignments [30]

--mq-cutoff INT The cutoff of mapping quality for keeping the alignment [50]

(only used when inputting more than one alignment files)

-ip FLOAT, --iden-percent FLOAT

Minimum identity (num_match_res/len_aln) of alignments [0.9]

-op FLOAT, --ovlp-percent FLOAT

Minimum overlapping percentage of the same read alignment if inputting more than one alignment files [0.9]

-cp FLOAT, --clip-percent FLOAT

Maximum clipped percentage of the alignment [0.1]

-fl INT, --flank-len INT

The flanking length of the clipped bases [15]

Plot Options:

-p, --plot Visualize the finally filtered whole genome (and regions if providing the option `-R`) depth [False]

-dmin FLOAT, --depth-min FLOAT

Minimum depth in folds of mean coverage for plotting [0.1]

-dmax FLOAT, --depth-max FLOAT

Maximum depth in folds of mean coverage for plotting [4.0]

-ws INT, --window-size INT

The window size when plotting [50000]

-it STR, --image-type STR

The format of the output images: png or pdf [png]

Other Options:

-f, --force Force rewriting of existing files [False]

-h, --help Show this help message and exit

-v, --version Show program's version number and exit

Examples:

python GCI.py -r ref.fa --hifi hifi.bam hifi.paf ... --nano nano.bam nano.paf ...

実行方法

レポジトリでは、canuのTrio binning assemblyを想定して作業を進めている。canuのTrio binning assemblyとは、2倍体ゲノムでハプロタイプを区別したアセンブリを行いたい時に2倍体の複雑な対立遺伝子のバリエーションにアセンブルが妨げられてしまうことに対処する方法。ハプロタイプゲノム配列を組み立てる前に対立遺伝子の変異を解消することで、ハプロタイプの組み立てを単純化する（引用）。具体的には、Trio binningは、両親それぞれのショートリードシークエンシングリードを用いて、まず子（F1）のロングシークエンシングリードをハプロタイプ固有のセットに分割する。その後、各ハプロタイプを独立にアセンブルし、完全な2倍体の再構成を行う。

GCIは、このTrio binning assemblyのデータが利用できること、あるいは両親のロングリードも利用出来る事（step1はskipできる）を想定して進めている。

１，Canuの”canu -haplotype”*1 をランする。両親のリードとbinning対象のF1のロングリードを指定する。レポジトリではゲノムサイズ”3g”となっているのでヒトゲノムアセンブリを想定していると思われる。分散コンピューティング環境ではないならuseGridはfalseとする。

canu -haplotype \
-p $prefix -d $dictionary \
genomeSize=3g \
maxThreads=20 \
-haplotypePat $pat \
-haplotypeMat $mat \
-pacbio-hifi $hifi \ ## binning ONT reads with -nanopore $ont
useGrid=false

-haplotype generate haplotype-specific reads

canu -haplotypeのラン後、レポジトリでは、父母に正しく分類できなかったリードもcatで結合して１つにすることを提案している。

２，binningしたロングリードを評価したいアセンブリへマッピング

評価したいF1のアセンブリにbinningしたロングリードをマッピングする。winnowmapとminimap2でそれぞれマッピングし、winnowmapからはソートしたbam、minimap2からはソートされたpafを得る。

#1 まずminimap2(paf出力がdefaultだがsamで書き出す)
minimap2 -t 20 -ax map-hifi $mat_asm $mat_hifi > ${mat.minimap2.hifi.sam}

#2 本プログラムはターゲット名でファイルをソートしないので、ソートしたpafファイルを作る
paftools.js sam2paf ${mat.minimap2.hifi.sam} | sort -k6,6V -k8,8n > ${mat.minimap2.hifi.paf}

#3 続いてwinnowmap 
meryl count k=15 output $mat_merylDB $mat_asm
meryl print greater-than distinct=0.9998 $mat_merylDB > $mat_repetitive_k15.txt
winnowmap -W $mat_repetitive_k15.txt -ax map-pb $mat_asm $mat_hifi > ${mat.winnowmap.hifi.sam}   ## mapping ONT reads with -ax map-ont

#4 bamとbam.baiも必要
samtools sort -@ 12 ${mat.winnowmap.hifi.sam} -o ${mat.winnowmap.hifi.bam}
samtools index ${mat.minimap2.hifi.bam}

winnowmapはリピートの多いゲノム配列へのロングリードの正確なend-to-endアラインメントを特徴としている（時間がやや長めにかかる）。

3、GCI.pyを使ってマッピングファイルのフィルタリングとゲノム連続性インデックスの取得。レポジトリではONTとpacbio両方指定している。すなわち、pacbioのbamとpaf、ONTのbamとpaf。

python GCI.py -r ref.fa --hifi hifi.bam hifi.paf --nano ont.bam ont.paf -d mat -o mat -p -it pdf

-r The reference file
--hifi PacBio HiFi reads alignment files (at least one bam file)
--nano Oxford Nanopore long reads alignment files (at least one bam file)
-d The directory of output files [.]
-o Prefix of output files [GCI]
-p Visualize the finally filtered whole genome (and regions if providing the option `-R`) depth [False]
-it The format of the output images: png or pdf [png]

出力例

outdir/

images/

mat.tig000000003.pdf

カバレッジプロット。水平の赤い破線は全ゲノム平均マッピングデプスを表す。緑はマッピングデプスを表す。複数のアライナーの結果が指定された場合、緑や青で色分けして表示される。水色とピンクの影はそれぞれ、高信頼性リードのマッピングサポートが低い（0.1*平均デプス未満）領域とサポートがない（深度ゼロ）領域を示す。

mat.gci

論文より

このツールは、厳格なフィルター基準をパスしたアラインメントを必要とする。マッピングされていないアラインメント、セカンダリーアラインメント、補足的アラインメントはすべて破棄される。さらに低品質なアラインメントを除去するために、mapping quality (<30 in default)、mapping identity (<90% in default)、clipped proportion (>10% in default)、が採用されている。
アライナー間のマッピングアルゴリズムによって生じるアラインメントバイアスの可能性を考慮し、少なくとも2つの一般的な配列アライナー（例：minimap2 (Li, 2018)、Winnowmap2 (Jain et al., 2022)、VerityMap (Bzikadze et al., 2022)）を同一データ上で実行することが推奨される。2つのアライナーによって生成されたアラインメントのうち、マッピングの品質要件に合格し、マッピング座標が一貫しているもの（デフォルトではオーバーラップ≧90%）が保持される。マッピングの精度と感度はアライナーによって大きく異なる。例えば、minimap2はWinnowmap2よりはるかに高速に実行されるが、マッピング品質が低く（通常はゼロに近い）、高度に反復した配列のアラインメントでは性能が劣る(Jain et al., 2022)。繰り返し領域のアラインメントをレスキューするため、1つのアライナーが高いマッピング品質（デフォルトでは50以上）を生成した場合、そのリードは保持される。
一連の厳密なアラインメントフィルトレーションの後、各塩基のマッピングカバレッジがカウントされる。GCIワークフローでは、まずクリッピングアラインメントを除外し、デプス推定を改善するために、各フランキングアラインメント（図では1塩基となっている）をトリミングする。
キュレーションされたマッピングデプスを使い、GCIは、デプスがゼロまたは極端に低い領域を潜在的なアセンブリの問題として報告する。物理的に近い問題（例えば距離が染色体長の0.5%未満）はマージされる。その後、元の染色体または配列は、リードアライメントがサポートされていない遺伝子座でキュレートされたコンティグに分割される。curatedアセンブリに対してcurated contig N50とcurated contig numberが計算される。最後に、コンティグN50の不連続性を考慮し、GCIはキュレーションされたアセンブリのコンティグN50値とコンティグ数の両方を統合し、GCIスコア（ゼロから100までのスケール）を用いて、真のギャップレスT2Tアセンブリに対するアセンブリの連続性のギャップを定量化する（論文図1）。コンティグN50値が飽和している場合でも、コンティグ数を用いてアセンブリ間の連続性の違いを定量することができる。真のT2Tアセンブリの場合、問題やギャップは検出されないため、コンティグN50は理論的な最大N50（染色体N50）に等しく、コンティグ数は染色体数に等しいため、このようなアセンブリのGCIスコアは100となる。