bamConcordanceは、PacificBiosciencesが管理しているレポジトリの１つで管理されている、リードのリファレンス配列とのマッピングの一致度からリードの経験的なクオリティ値を算出するpythonスクリプト。エラー修正ツールで修正された後のシークエンシングリードは、通常、元のシークエンシングリードにはあったquality値がなくなってfasta出力される（互換性がなくなるため）。bamConcordanceは、このようなクオリティ値が失われたシークエンシングリードのクオリティ品質を調べる時などに使用される。類似のツールにBBduk（紹介）のクオリティ再校正機能（link）があるが、BBdukは、fastqを入力として元のfastqのクオリティ品質を（マッピングベースで）リキャリブレートする目的で使われる。よって入力はfastqである必要があり、結果は元のクオリティ値に大きく依存している（pseudo fastqを作って使うとおかしな結果を出す）。bamConcordanceはqualityクオリティ情報が失われたシークエンシングリードを使って、正しいリファレンスへのマッピングから"リード単位"でのクオリティ値を推定する。
インストール
#pysamが必要
mamba install -c bioconda pysam -y
#本体
git clone https://github.com/PacificBiosciences/hg002-ccs.git
cd hg002-ccs/concordance/
> ./bamConcordance -h
usage: Measure concordance of read alignments to a reference genome [-h] [--hcregions hcregions.bed.gz] [--hcvariants hcvariants.vcf.gz] [--chrom S] [--outbam out.bam] ref.fasta in.bam out.csv
positional arguments:
ref.fasta
in.bam
out.csv Per aligment concordance statistics, CSV
options:
-h, --help show this help message and exit
--hcregions hcregions.bed.gz
High-confidence regions, BED.gz (.tbi required)
--hcvariants hcvariants.vcf.gz
High-confidence variants, VCF.gz (.tbi required)
--chrom S Limit analysis to the specified chromosome
--outbam out.bam Input BAM with concordance annotations for relevant records
CONCORDANCE
Measure the concordance of alignments to the reference genome.
Consider only non-variant positions in high-confidence regions.
実行方法
１，リファレンスゲノムにリードをマッピングしてbamを作成する。
minimap2 -t20 -ax sr --eqx reference.fasta sample_R* |samtools sort -@ 4 > reference.bam
minimap2やbwa1,2を使う場合、CIGARをマッチとミスマッチで区別するため、--eqxを指定する（--eqxをつけると、CIGARのマッチに=、ミスマッチにXが使用される。CIGARがマッチとミスマッチ両方Mだと、bamConcordanceがエラーが起こす）。
2、bamConcordanceを実行する。
bamConcordance reference.fasta reference.bam out.csv
出力
csvlookで出力を見てみる（pip install csvkit）。
csvlook out.csv |less
情報が少ないが、リードごとに１行ずつ、以下のフィールドがプリントされていると推測される。
- read - リード名
- readLengthBp - リード長
- effectiveCoverage - empty
- subreadPasses - empty 同じリードを何回測定したか）
- predictedConcordance - empty (リードの予測精度）
- alignmentType - Primary, Supplementary
- alignmentMapq - MAPQ (アラインメントの品質スコア）
- hcReadLengthBp - 高品質リード長
- concordance - 実際の一致率
- concordanceQv - 一致率の品質値（品質スコア）
- mismatchBp - ミスマッチの塩基数
- nonHpInsertionBp - 非ホモポリマー挿入塩基数
- nonHpDeletionBp - 非ホモポリマー欠失塩基数
- hpInsertionBp - ホモポリマー挿入塩基数
- hpDeletionBp - ホモポリマー欠失塩基数
引用
https://github.com/PacificBiosciences/hg002-ccs
関連
・リファレンスアセンブリにアライメントした後のリードの品質を評価する best
・Nanoporeのロングリードのアセンブリを中心としたツールキット Pomoxis
