アラインメントのPAFファイルを扱うユーティリティ paftools

2019,11,15 paftools.js call のコマンドエラー修正

2019 11/26 vatiantコールコマンドのエラー修正

2021 7/21 コマンド例追加

Minimap2には、PAFフォーマットのアライメントを処理する（java）スクリプトpaftools.jsが付属している。paftoolsを使うことで、 アセンブリをリファレンスゲノムにアラインメントしてバリアントをコールしたり、PAF/SAMからBEDなどのフォーマットに変換したりできる。

インストール 

Github

#minimap2インストールで導入できる。（link）
mamba install -c bioconda -y minimap2

#またはソースからビルドする（Github参照）

> paftools.js

# paftools.js

Usage: paftools.js <command> [arguments]

Commands:

  view       convert PAF to BLAST-like (for eyeballing) or MAF

  splice2bed convert spliced alignment in PAF/SAM to BED12

  sam2paf    convert SAM to PAF

  delta2paf  convert MUMmer's delta to PAF

  gff2bed    convert GTF/GFF3 to BED12

  stat       collect basic mapping information in PAF/SAM

  asmstat    collect basic assembly information

  asmgene    evaluate gene completeness (EXPERIMENTAL)

  liftover   simplistic liftOver

  call       call variants from asm-to-ref alignment with the cs tag

  bedcov     compute the number of bases covered

  version    print paftools.js version

  mapeval    evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ

  mason2fq   convert mason2-simulated SAM to FASTQ

  pbsim2fq   convert PBSIM-simulated MAF to FASTQ

  junceval   evaluate splice junction consistency with known annotations

  ov-eval    evaluate read overlap sensitivity using read-to-ref mapping

テストラン

git clone https://github.com/lh3/minimap2.git
cd minimap2/test/

1、paftools.js view  - アラインメント表示

Usage: paftools.js view [options] <in.paf>

Options:

  -f STR    output format: aln (BLAST-like), maf or lastz-cigar [aln]

  -l INT    line length in BLAST-like output [80]

minimap2 --cs MT-human.fa MT-orang.fa | paftools.js view -

#実際にはlessなどに渡したり保存（"> alignment"）して見る
minimap2 --cs MT-human.fa MT-orang.fa | paftools.js view - |less

• --cs   output the cs tag; STR is 'short' (if absent) or 'long' [none]

2、paftools.js stat  - アラインメントstatistics

Usage: paftools.js stat [-l gapOutLen] <in.sam>|<in.paf>

minimap2 -c MT-human.fa MT-orang.fa | paftools.js stat -

• -c    output CIGAR in PAF

出力

Number of mapped sequences: 1

Number of primary alignments: 1

Number of secondary alignments: 0

Number of primary alignments with >65535 CIGAR operations: 0

Number of bases in mapped sequences: 16499

Number of mapped bases: 16025

Number of insertions in [0,50): 51

Number of insertions in [50,100): 0

Number of insertions in [100,300): 0

Number of insertions in [300,400): 0

Number of insertions in [400,1000): 0

Number of insertions in [1000,inf): 0

Number of deletions in [0,50): 50

Number of deletions in [50,100): 0

Number of deletions in [100,300): 0

Number of deletions in [300,400): 0

Number of deletions in [400,1000): 0

Number of deletions in [1000,inf): 0

3、paftools.js call  - バリアントコール

Usage: sort -k6,6 -k8,8n <with-cs.paf> | paftools.js call [options] -

Options:

  -l INT    min alignment length to compute coverage [10000]

  -L INT    min alignment length to call variants [50000]

  -q INT    min mapping quality [5]

  -g INT    short/long gap threshold (for statistics only) [50]

  -f FILE   reference sequences (enabling VCF output) [null]

  -s NAME   sample name in VCF header [sample]

#sortコマンドに渡しcoordinate sortしてからバリアントコールする
minimap2 -c --cs test/MT-human.fa test/MT-orang.fa \
 | sort -k6,6 -k8,8n \
 | paftools.js call -L 15000 - 

#haploid assemblyのバリアントコール。de novoアセンブリ配列とリファレンスを比較し、結果を保存。
minimap2 -cx asm5 -t8 --cs ref.fa asm.fa \
 | sort -k6,6 -k8,8n - \
 | paftools.js call - > variant

#vcf output (ここではmap-ontプリセット)（100bp以上の違いのみ）
minimap2 -cx map-ont -t8 --cs ref.fa targets.fa | sort -k6,6 -k8,8n - | paftools.js call -f ref.fa -L 100 -l 100 - > variant.vcf

• -L     min alignment length to call variants [50000] 
• -x asm5|asm10|asm20    asm-to-ref mapping, for ~0.1/1/5% sequence divergence
• -f FILE   reference sequences (enabling VCF output) [null]

出力最後尾

V MT_human 16528 16529 1 60 t c MT_orang 15983 15984 +

V MT_human 16531 16531 1 60 - g MT_orang 15986 15987 +

V MT_human 16533 16534 1 60 a c MT_orang 15989 15990 +

15993 reference bases covered by exactly one contig

2150 substitutions; ts/tv = 3.498

36 1bp deletions

34 1bp insertions

9 2bp deletions

5 2bp insertions

5 [3,50) deletions

12 [3,50) insertions

0 [50,1000) deletions

0 [50,1000) insertions

0 >=1000 deletions

0 >=1000 insertions

4、paftools.js liftover  - リフトオーバー

Usage: paftools.js liftover [options] <aln.paf> <query.bed>

Options:

  -q INT    min mapping quality [5]

  -l INT    min alignment length [50000]

  -d FLOAT  max sequence divergence (>=1 to disable) [1]

クエリのbedファイルを指定する。

minimap2 -c MT-human.fa MT-orang.fa \
 | paftools.js liftover -l 10000 - <(echo -e "MT_orang\t2000\t5000")

• -l     min alignment length [50000] 

5、paftools.js junceval  - アノテーションとスプライスジャンクションを比較

Usage: paftools.js junceval [options] <gene.gtf> <aln.sam>

Options:

  -l INT    tolerance of junction positions (0 for exact) [0]

  -p        print overlapping introns

  -e        print erroreous overlapping introns

  -c        only consider alignments to /^(chr)?([0-9]+|X|Y)$/

gtfとマッピングのsam（pafではない）を指定する。

paftools.js junceval -e anno.gtf splice.sam > out.txt

6、paftools.js splice2bed  - GTF/GFF3をBED12フォーマットに変換

Usage: paftools.js splice2bed [options] <in.paf>|<in.sam>

Options:

  -m      keep multiple mappings (SAM flag 0x100)

 gtf/gffを指定する。

paftools.js splice2bed anno.gtf > anno.bed

７、paftools.js sam2paf   - SAMをBED12フォーマットに変換

Usage: paftools.js sam2paf [options] <in.sam>

Options:

  -p      convert primary or supplementary alignments only

  -L      output the cs tag in the long form

paftools.js sam2paf input.sam > output.bed

引用

Minimap2: pairwise alignment for nucleotide sequences.

Li H

Bioinformatics. 2018 May 10.

関連