2020 4/23 論文追記
以前、ロングリードのアセンブリ前処理ツール yacrdを紹介した。
今回はアセンブリ過程でフィルタリングして出力を調節するfpaを紹介する。
以下のフィルタリングが行える (Githubより)。
- internal match
- containment
- dovetails
- self matching
- read name match against regex
- length of overlap
- length of read in overlap
Pierre MarijonさんのBlog (More details and some experiment are present)
https://blog.pierre.marijon.fr/binary-mapping-format/
インストール
ubuntu16.0.4のminicona3.4.0.5環境でテストした。
ビルド依存
- Rust >= 1.32
- libgz
- libbzip2
- liblzma
その他 (option)
- minimap2
- miniasm
本体 Github
#bioconda (link) (*1)
conda install -c bioconda -y fpa
> fpa -h
# fpa -h
fpa 0.5 Krabby
Pierre Marijon <pierre.marijon@inria.fr>
fpa take long read mapping information and filter them
USAGE:
fpa [OPTIONS] [SUBCOMMAND]
FLAGS:
-h, --help Prints help information
-V, --version Prints version information
OPTIONS:
-z, --compression-out <compression-out>
Output compression format, the input compression format is chosen by default [possible values: gzip, bzip2,
lzma, no]
-F, --format <format> Force the format used [possible values: paf, m4]
-i, --input <input> Path to input file, use '-' for stdin [default: -]
--internal-threshold <internal-match-threshold>
A match is internal match if overhang length > match length * internal threshold this option set internal
match [default: 0.8]
-o, --output <output> Path to output file, use '-' for stdout [default: -]
SUBCOMMANDS:
drop fpa drop mapping match this constraints
gfa fpa generate a overlap graph in gfa1 format with mapping passing filter
help Prints this message or the help of the given subcommand(s)
index fpa generate a index of mapping passing filter
keep fpa keep only mapping match this constraints
rename fpa rename reads with name you chose or with incremental counter
> fpa keep -h
# fpa keep -h
fpa-keep
fpa keep only mapping match this constraints
USAGE:
fpa keep [FLAGS] [OPTIONS] [SUBCOMMAND]
FLAGS:
-c, --containment Keep only containment mapping
-d, --dovetail Keep only dovetail mapping
-h, --help Prints help information
-i, --internalmatch Keep only internal mapping
-m, --same-name Keep only mapping where reads have same name
-V, --version Prints version information
OPTIONS:
-l, --length-lower <length_lower> Keep only mapping with length lower than value
-L, --length-upper <length_upper> Keep only mapping with length upper than value
-n, --name-match <name_match> Keep only mapping where one reads match with regex
-s, --sequence-length-lower <sequence_length_lower>
Keep only mapping where one reads have length lower than value
-S, --sequence-length-upper <sequence_length_upper>
Keep only mapping where one reads have length upper than value
> fpa drop -h
# fpa drop -h
fpa-drop
fpa drop mapping match this constraints
USAGE:
fpa drop [FLAGS] [OPTIONS] [SUBCOMMAND]
FLAGS:
-c, --containment Drop containment mapping
-d, --dovetail Drop dovetail mapping
-h, --help Prints help information
-i, --internalmatch Drop internal mapping
-m, --same-name Drop mapping where reads have same name
-V, --version Prints version information
OPTIONS:
-l, --length-lower <length_lower> Drop mapping with length lower than value
-L, --length-upper <length_upper> Drop mapping with length upper than value
-n, --name-match <name_match> Drop mapping where one reads match with regex
-s, --sequence-length-lower <sequence_length_lower> Drop mapping where one reads have length lower than value
-S, --sequence-length-upper <sequence_length_upper> Drop mapping where one reads have length upper than value
> fpa gfa -h
# fpa gfa -h
fpa-gfa
fpa generate a overlap graph in gfa1 format with mapping passing filter
USAGE:
fpa gfa [FLAGS] [OPTIONS] [SUBCOMMAND]
FLAGS:
-c, --containment Keep containment overlap
-h, --help Prints help information
-i, --internalmatch Keep internal match overlap
-V, --version Prints version information
OPTIONS:
-o, --output <output> Write mapping passing filter in gfa1 graph format in path passed as parameter
> fpa index -h
# fpa index -h
fpa-index
fpa generate a index of mapping passing filter
USAGE:
fpa index [OPTIONS] [SUBCOMMAND]
FLAGS:
-h, --help Prints help information
-V, --version Prints version information
OPTIONS:
-f, --filename <filename> Write index of mapping passing filter in path passed as parameter
-t, --type <type> Type of index, only reference read when it's query, target or both of them [default:
both] [possible values: query, target, both]
> fpa rename -h
# fpa rename -h
fpa-rename
fpa rename reads with name you chose or with incremental counter
USAGE:
fpa rename [OPTIONS] [SUBCOMMAND]
FLAGS:
-h, --help Prints help information
-V, --version Prints version information
OPTIONS:
-i, --input <input> Rename reads with value in path passed as parameter
-o, --output <output> Write rename table in path passed as parameter
実行方法
1、 keep - 制約に合致するマッピングを出力
minimap/miniasmのロングリードアセンブリの途中にfpaのフィルタリングを加えるように設計されている。そのため、fpaはpafファイルを入力とし、pafファイルを出力する。miniasmによるRaw assemblyの第一ステップはminimapまたはminimap2によるリード同士のオーバーラップ（all-vs-all read self-mappings）なので、この第一ステップの出力（minimapの出力はdefaultではPAF）を受け取り、fpaでフィルタリングしてminiasmでアセンブリ（コンカテネートしてunitigをGFA出力（エラーレートはinputのリードと同じ））を行う。
例えばONTのrawアセンブリを以下の流れで実行するとする。
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | gzip -1 > reads.paf.gz
miniasm -f nanopore.fq.gz reads.paf.gz > raw_assebly.gfa
#この流れにfpaを組み込む。"fpa drop -l <NUM>"で指定以下のリードを捨てて出力。
#1000bp以上のリードだけpaf出力。
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa drop -l 1000 | gzip - > filtered.paf.gz
#それからminiasmを実行する
#fpa dropなら逆の条件でフィルタリングできる。fpa keep -l <NUM>で指定以下のリードのみキープして出力
#10000bp以下のリードだけpaf出力。
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa keep -l 10000 | gzip - > filtered.paf.gz
#つまりkeepとdropは逆の挙動をする。1000bp以上のリードだけキープして出力するなら"keep -L <NUM>"か"drop -l <NUM>"。keepを使うなら
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa keep -L 1000 | gzip - > filtered.paf.gz
サイズセレクト
1000bp以上10000-bp以下のリードのみにフィルタリング。
#keepを使うなら
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa keep -L 1000 -l 10000| gzip - > filtered.paf.gz
#dropを使うなら
minimap2 -t 18 -x ava-ont nanopore.fq.gz nanopore.fq.gz | fpa drop -l 1000 -L 10000| gzip - > filtered.paf.gz
dovetail（蟻継ぎ）のオーバーラップを持つリードのみにフィルタリング。
minimap2 long_read.fasta long_read.fasta | fpa keep -d | gzip - > filtered.paf.gz
- -d Keep only dovetail mapping
インターナルのオーバーラップを持つリードのみにフィルタリング。
minimap2 long_read.fasta long_read.fasta | fpa keep -i | gzip - > filtered.paf.gz
- -i Keep only internal mapping
汚染配列を除去。
minimap2 long_read.fasta long_read.fasta | fpa keep -c | gzip - > filtered.paf.gz
- -c Keep only containment mapping
汚染配列除去とサイズセレクト。
minimap2 long_read.fasta long_read.fasta | fpa keep -c keep -L 1000 | gzip - > filtered.paf.gz
- -c Keep only containment mapping
他の例はGithubを参考にしてください。
fpa gfaコマンドも組み込めば、アセンブリグラフのGFA（version1）で出力できる。bandage等を使えばグラフを可視化できる。
minimap2 long_read.fasta long_read.fasta | fpa keep -d gfa -o dovetail.gfa > dovetail.paf
あくまでアセンブリ前提に設計されているので、GFAの出力（-o）以外にリダイレクトでPAFも出力する。fpa gfa -oで出力したGFAに、配列情報は含まれないので注意する。
fpaはロングリードのアセンブリを前提に設計されているので、アセンブリ前にロングリードを前処理したいならyacrdを使う。
引用
yacrd and fpa: upstream tools for long-read genome assembly
Pierre Marijon, Rayan Chikhi, Jean-Stéphane Varré
bioRxiv preprint first posted online Jun. 18, 2019
yacrd and fpa: upstream tools for long-read genome assembly
Pierre Marijon, Rayan Chikhi, Jean-Stéphane Varré
Bioinformatics, 2020
*1
2019 7/4現在、macosだと0.4までしか入らない。v0.4だとkeepなどの一部コマンドが存在しない。 keepなどを使いたければソースからビルドする。
