macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

De Novo Variantsを正確に発見するためのマッピング不要のフレームワーク Kevlar

 

 遺伝性変異は複雑な遺伝性疾患における主要な寄与因子であると推測されている。多くの遺伝性疾患の遺伝率は比較的高いと推定されている。例えば、自閉症スペクトラム障害ASD)の遺伝率は0.6を超え、統合失調症の遺伝率は0.5を超える。この遺伝性のごく一部のみが既知の遺伝性変異によって説明されており、"missing heritability"と名付けられた現象が存在する(ref.3 link)。一つの仮説はde novo突然変異である(特にindelと構造変異)。しかし、特にSV発見の複雑さから、これらの障害に対するそれらの寄与の説明は不完全になっている。一般的な遺伝性変異、特にde novoバリアントの発見は、依然として強い研究関心のあるトピックである。複雑な病因における遺伝性変異の役割を明らかにすることに加えて、多くのサンプルまたはコホートにわたるde novoバリアントの発見およびカタログ化の改善は、ヒトゲノムにおける重要な未解決問題にさらなる光を投げかける。
 単一家系の全ゲノムシークエンシング(遺伝的疾患の孤立した症例を提示)は、生殖系列のde novo変異の発見から証明された成功例である。次に、ミスマッチ、ギャップ、カバレッジの急激なシフト、ペアリードのインサートサイズやペアリードの向きの不一致など、リードのアライメントに見られるアーティファクトに基づいて、各個人のバリアントが予測される。これは、継承されたバリエーション、真のデノボバリエーション、および偽のバリアントコールと識別するために検査する必要がある。
 リファレンスに基づいた変異発見法は複雑な遺伝性疾患の研究において価値があることが証明されているが、本著者らはそれらの限界のいくつかに注目する。リードアラインメントアルゴリズムの一貫した改善にもかかわらず、各リードについて正しいマッピングを見つけることは、シーケンシングエラー、 repetitive DNA含有量、およびリファレンスにおけるミスアセンブリによって依然として複雑である。また、変異ブレイクポイントにまたがるか、または新規配列を含むためにリファレンスゲノムにマッピングされないリードは、マッピングベースの変異予測器では完全に無視される。最後に、リードアラインメントの分析によって決定されたほとんどのバリアントコールは、関心のある個人(子供または発端者)に固有のものではなく、代わりにファミリーとリファレンスゲノムドナーとの間の先祖の相違を反映している。ヒト生殖細胞系突然変異率の推定値は、一世代あたりおよそ80の新規突然変異という予想を与え、そして何百万もの遺伝的または偽のバリアントから真のde novo変異イベントを区別することは大きな課題である。

 より一般的には、正確かつ包括的なde novoバリアントの発見は、いくつかの計算上および生物学的要因により複雑であり、依然としてとらえどころのない目標である。proband(発端者)にバリアントが存在することだけでなく、両親には存在しないことについてどのアルゴリズムも確信している必要がある。そして、一塩基変異(SNV)が最も一般的な変異のタイプであるが、よりまれなより大きな変異は全体的により多くのヌクレオチドに影響を及ぼし、そして遺伝性障害ではより大きな効果を有すると仮定される。 indelとSVの予測の本質的な複雑さのために、これらのより大きなde novoバリアントの正確な発見は特にチャレンジングになる。リファレンスマッピングのコンテキストでは、indelを確実にコールするには、すべてのギャップを一貫してアライメントした状態で、各リードをindelにまたがって正確にマッピングする必要がある。これは短いindelsに対してのみ可能であり、エラーやポジションギャップを起こしがちである。したがって、長さ10 bpを超えるindelsの予測は非常に困難であり、高い偽陽性率および高い偽陰性率を伴うことが証明されている。さらに、SVの予測は、リードデプス変化またはリードペアのシグネチャのような間接的なシグネチャによってのみ可能である。これらのシグネチャは非常にノイズが多くなる可能性があり、偽陽性および偽陽性の予測率が高くなる。結果として、それらの発生率を含むデノボSVのいくつかの基本的な特性は未知のままである。逆位を伴うduplicationなど、より複雑なタイプのde novo SVを予測する方法も存在しないことに注意することが重要である。
 新規のバリアント予測に関する課題の多くは、リファレンスゲノムを介して間接的にではなく、直接関連する個体間のシーケンシング内容を比較するアプローチによって軽減することができる。このようなアプローチでは、リードの調整は不要である。また、offターゲットの共有/継承されたバリエーションには影響されない。マッピングフリーのアプローチが必要とするのは、アライメントアーチファクトに関して定義されていないバリエーションのシグネチャである。このマッピングフリーのコンテキストでは、k-mer(固定長kのシーケンス)の分析を介して、シーケンスの内容と量に関してバリアントを定義する。最近、ガン腫瘍における体細胞変異の正確な発見および全エキソームシーケンシングデータの完全な発見のために同様のアイデアが利用されている(ref.18, 19)。
 直感的には、de novo生殖細胞系変異は、親ゲノムと比較してコーラーに新しい配列内容をもたらすはずである。最も単純な場合でさえ、十分に大きなkの値を考えると、突然変異をまたぐ k-merの全部ではないにしても、ほとんどの場合、単一のヌクレオチド置換はユニークであるべきである。ちなみに、これはindelやさまざまなタイプの構造的なバリエーションなど、他のクラスのバリエーションにも当てはまる。そして、probandゲノムは十分深くサンプリングされており、よってこれらの新規なk-merはシーケンシングデータ中に豊富に存在すると期待される。したがって、単一のモデルを使用して、SNVおよびより大きなバリアント(indel、構造変化)の両方を同時に検出することが可能であるはずである。
 この直感に基づいて、本著者らは、新規のバリアント発見問題のマッピングフリーな定式化に基づく新しい方法であるKevlarを開発した。Kevlarは、両親からのリードには事実上欠けているが、proband/子供のデータのみに有意に高い存在量を持つ「興味深い」k-merを特定するために、k-merの存在量を調べる。これらのk-merは、デノボバリアントの潜在的存在の指標である。次に、ペアリード間で共有される新規のk-merに基づいてリードをグループ分けする。このそれぞれは推定バリアントを反映している。その後、Kevlarは標準的なアルゴリズムを使用して、リードの各セットをコンティグにアセンブルし、アセンブルされたコンティグをリファレンスゲノムにアライメントさせて予備的なバリアントコールを行う。最後に、Kevlarはバリアント予測をスコアリングし、誤ってコールされた遺伝性変異と真の de novo変異を区別するために新規な確率モデルを使用する。
 論文では、シミュレートされたデータとリアルデータを使いこの新しい方法の有用性を実証する。 KevlarがSNVおよびshort indel発見のためのクラス最高のツールと同様の予測精度を達成すると同時に、より大きなイベントを高い精度で予測することを示す。ヌクレオチドレベルの精度でブレイクポイントを定義し、大規模なSVイベントを正確に予測するケブラーの能力を実証する。
 Kevlarはオープンソースソフトウェアプロジェクトとして入手可能で、Python APIコマンドラインインターフェース、または標準のSnakemakeワークフローを介して呼び出すことができる。

 

f:id:kazumaxneo:20190217224907p:plain

The Kevlar pipeline. Preprintより転載。

 

Documentation

https://kevlar.readthedocs.io/en/latest/

Quick start guide

Quick start — kevlar 0.7rc1+1.ge55ba25 documentation

Tutorial

https://github.com/kevlar-dev/kevlar/tree/master/kevlar/workflows/mark-I

 

インストール

依存

The kevlar package requires Python 3 and has several dependencies that are not in the standard Python libraries.

  • the networkx package
  • the pysam module
  • the pandas library
  • the scipy library
  • the intervaltree library
  • the khmer package

Also, kevlar requires the bwa command to be callable from your $PATH environmental variable.

 本体 Github

mamba create -n Kevlar python=3.6 -y
conda activate Kevlar

pip3 install pysam networkx pandas scipy intervaltree git+https://github.com/dib-lab/khmer.git
pip3 install biokevlar

#snakemakeも必要
pip3 install snakemake

 

依存も含めて正しくインストールされているかテスト

pip3 install pytest
pytest --pyargs kevlar.tests

=========================================================================================================== test session starts ============================================================================================================

platform linux -- Python 3.6.7, pytest-4.4.0, py-1.8.0, pluggy-0.9.0

rootdir: /

collected 373 items                                                                                                                                                                                                                        

 

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_alac.py ...........................                                                                                                                                          [  7%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_assemble.py ...............                                                                                                                                                  [ 11%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_augment.py .....                                                                                                                                                             [ 12%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_call.py ...........................                                                                                                                                          [ 19%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_cigar.py .........                                                                                                                                                           [ 22%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_cli.py ....................                                                                                                                                                  [ 27%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_count.py ...........................                                                                                                                                         [ 34%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_dist.py ..........                                                                                                                                                           [ 37%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_evaluate.py .                                                                                                                                                                [ 37%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_filter.py ......                                                                                                                                                             [ 39%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_gentrio.py ............................                                                                                                                                      [ 46%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_localize.py ............................                                                                                                                                     [ 54%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_mutablestring.py .....                                                                                                                                                       [ 55%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_mutate.py .........                                                                                                                                                          [ 58%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_novel.py ................                                                                                                                                                    [ 62%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_partition.py .......                                                                                                                                                         [ 64%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_progress.py ..                                                                                                                                                               [ 64%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_readgraph.py ..                                                                                                                                                              [ 65%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_readpair.py .........                                                                                                                                                        [ 67%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_reference.py .....                                                                                                                                                           [ 69%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_seqio.py ..............                                                                                                                                                      [ 72%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_simlike.py ....................................                                                                                                                              [ 82%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_sketch.py .................                                                                                                                                                  [ 87%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_split.py ..                                                                                                                                                                  [ 87%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_timer.py .                                                                                                                                                                   [ 87%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_unband.py ...                                                                                                                                                                [ 88%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_varfilter.py ....                                                                                                                                                            [ 89%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_varmap.py ..........................                                                                                                                                         [ 96%]

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_vcf.py ............                                                                                                                                                          [100%]

 

============================================================================================================= warnings summary =============================================================================================================

usr/local/lib/python3.6/dist-packages/kevlar/tests/test_dist.py::test_tsv

  /usr/local/lib/python3.6/dist-packages/kevlar/tests/test_dist.py:120: FutureWarning: read_table is deprecated, use read_csv instead.

    data = pandas.read_table(distfile.name, sep='\t')

 

-- Docs: https://docs.pytest.org/en/latest/warnings.html

================================================================================================= 373 passed, 1 warnings in 70.72 seconds ==================================================================================================

ヘルプ

> kevlar count -h

# kevlar count -h

usage: kevlar count [-h] [-k K] [-c C] [-M MEM] [--max-fpr FPR] [--mask MSK]

                    [--count-masked] [--num-bands N] [--band I] [-t T]

                    counttable seqfile [seqfile ...]

 

Compute k-mer abundances for the provided sample. Supports k-mer banding:

see http://kevlar.readthedocs.io/en/latest/banding.html for more details.

 

positional arguments:

  counttable            name of the file to which the output (a k-mer count

                        table) will be written; the suffix ".counttable" will

                        be applied if the provided file name does not end in

                        ".ct" or ".counttable"

  seqfile               input files in Fastq/Fasta format

 

optional arguments:

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

  -k K, --ksize K       k-mer size; default is 31

  -c C, --counter-size C

                        number of bits to allocate for counting each k-mer;

                        options are 1 (max count: 1), 4 (max count: 15), and 8

                        (max count: 255); default is 8

  -M MEM, --memory MEM  memory to allocate for the count table

  --max-fpr FPR         terminate if the estimated false positive rate for any

                        sample is higher than "FPR"; default is 0.2

  --mask MSK            counttable or nodetable of k-mers to ignore when

                        counting k-mers

  --count-masked        by default, when a mask is provided k-mers in the mask

                        are ignored; this setting inverts the behavior so that

                        only k-mers in the mask are counted

  --num-bands N         number of bands into which to divide the hashed k-mer

                        space

  --band I              a number between 1 and N (inclusive) indicating the

                        band to be processed

  -t T, --threads T     number of threads to use for file processing; default

                        is 1

 

Example::

 

    kevlar count --memory 500M case1.ct case1-reads.fastq

 

Example::

 

    kevlar count --ksize 25 --memory 12G --max-fpr 0.01 --threads 8 \

        proband.counttable \

        proband-R1.fq.gz proband-R2.fq.gz proband-unpaired.fq.gz

> kevlar novel -h

# kevlar novel -h

usage: kevlar novel --case F [F ...] [--case-counts F [F ...]]

                    [--control F [F ...]] [--control-counts F [F ...]] [-x X]

                    [-y Y] [-M MEM] [--max-fpr FPR] [--num-bands N] [--band I]

                    [-o FILE] [--save-case-counts CT [CT ...]]

                    [--save-ctrl-counts CT [CT ...]] [-h] [-k K]

                    [--abund-screen INT] [-t T] [--skip-until ID]

 

Identify "interesting" (potentially novel) k-mers and output the

corresponding reads. Here we define "interesting" k-mers as those which are

high abundance in each case sample and effectively absent (below some

specified abundance threshold) in each control sample.

 

Case/control config:

  Specify input files, as well as thresholds for selecting "interesting"

  k-mers. A single pass is made over input files for control samples (to

  compute k-mer abundances), while two passes are made over input files for

  case samples (to compute k-mer abundances, and then to identify

  "interesting" k-mers). The k-mer abundance computing steps can be skipped

  if pre-computed k-mer abunandances are provided using the "--case-counts"

  and/or "--control-counts" settings. If "--control-counts" is declared, then

  all "--control" flags are ignored. If "--case-counts" is declared,

  FASTA/FASTQ files must still be provided with "--case" for selecting

  "interesting" k-mers and reads.

 

  --case F [F ...]      one or more FASTA/FASTQ files containing reads from a

                        case sample; can be declared multiple times

                        corresponding to multiple case samples, see examples

                        below

  --case-counts F [F ...]

                        counttable file(s) corresponding to each case sample;

                        if not provided, k-mer abundances will be computed

                        from FASTA/FASTQ input; only one counttable per

                        sample, see examples below

  --control F [F ...]   one or more FASTA/FASTQ files containing reads from a

                        control sample; can be declared multiple times

                        corresponding to multiple control samples, see

                        examples below

  --control-counts F [F ...]

                        counttable file(s) corresponding to each control

                        sample; if not provided, k-mer abundances will be

                        computed from FASTA/FASTQ input; only one counttable

                        per sample, see examples below

  -x X, --ctrl-max X    k-mers with abund > X in any control sample are

                        uninteresting; default is X=1

  -y Y, --case-min Y    k-mers with abund < Y in any case sample are

                        uninteresting; default is Y=6

  -M MEM, --memory MEM  total memory allocated to k-mer abundance for each

                        sample; default is 1M; ignored when pre-computed k-mer

                        abundances are supplied via counttable

  --max-fpr FPR         terminate if the expected false positive rate for any

                        sample is higher than the specified FPR; default is

                        0.2

 

K-mer banding:

  If memory is a limiting factor, it is possible to get a linear decrease in

  memory consumption by running `kevlar novel` in "banded" mode. Splitting

  the hashed k-mer space into N bands and only considering k-mers from one

  band at a time reduces the memory consumption to approximately 1/N of the

  total memory required. This implements a scatter/gather approach in which

  `kevlar novel` is run N times, after the results are combined using

  `kevlar filter`.

 

  --num-bands N         number of bands into which to divide the hashed k-mer

                        space

  --band I              a number between 1 and N (inclusive) indicating the

                        band to be processed

 

Output settings:

  -o FILE, --out FILE   file to which interesting reads will be written;

                        default is terminal (stdout)

  --save-case-counts CT [CT ...]

                        save the computed k-mer counts for each case sample to

                        the specified count table file(s)

  --save-ctrl-counts CT [CT ...]

                        save the computed k-mer counts for each control sample

                        to the specified count table file(s)

 

Miscellaneous settings:

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

  -k K, --ksize K       k-mer size; default is 31

  --abund-screen INT    discard reads with any k-mers whose abundance is < INT

  -t T, --threads T     number of threads to use for file processing; default

                        is 1

  --skip-until ID       when re-running `kevlar novel`, skip all reads in the

                        case input until read with name `ID` is observed

 

Example::

 

    kevlar novel --out novel-reads.augfastq --case proband-reads.fq.gz \

        --control father-reads-r1.fq.gz father-reads-r2.fq.gz \

        --control mother-reads.fq.gz

 

Example::

 

    kevlar novel --out novel-reads.augfastq.gz \

        --control-counts father.counttable mother.counttable \

        --case-counts proband.counttable --case proband-reads.fastq \

        --ctrl-max 0 --case-min 10 --ksize 27

 

Example::

 

    kevlar novel --out output.augfastq \

        --case proband1.fq --case proband2.fq \

        --control control1a.fq control1b.fq \

        --control control2a.fq control2b.fq \

        --save-case-counts p1.ct p2.ct --save-ctrl-counts c1.ct c2.ct

> kevlar filter -h

# kevlar filter -h

usage: kevlar [-h] [-v] [-l F] [--tee] cmd ...

kevlar: error: argument cmd: invalid choice: 'filter\udcc2' (choose from 'count', 'novel', 'filter', 'augment', 'assemble', 'mutate', 'gentrio', 'partition', 'localize', 'call', 'alac', 'varfilter', 'simlike', 'split', 'dist', 'unband')

> kevlar partition -h

# kevlar partition -h

usage: kevlar partition [-h] [-s] [--min-abund X] [--max-abund Y] [--no-dedup]

                        [--gml FILE] [--split OUTPREFIX] [-o FILE]

                        infile

 

Construct a graph to group reads by shared interesting k-mers. Nodes in the

graph represent reads, and edges between a pair of nodes indicate that the two

corresponding reads have one or more interesting k-mers in common. Connected

components in the undirected graph correspond to distinct variants (or

variant-related breakpoints).

 

positional arguments:

  infile               input reads in augmented Fast[q|a] format

 

optional arguments:

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

  -s, --strict         require perfect identity between overlapping reads for

                       inclusion in the same partition; by default, only a

                       shared interesting k-mer is required

  --min-abund X        ignore k-mers with abundance lower than X; default is 2

  --max-abund Y        ignore k-mers with abundance higher than Y; default is

                       200

  --no-dedup           skip step to remove duplicates

  --gml FILE           write read graph to .gml file

  --split OUTPREFIX    write each partition to a separate output file, each

                       with a filename like "OUTPREFIX.cc#.augfastq.gz"

  -o FILE, --out FILE  output file; default is terminal (stdout)

> kevlar call -h

# kevlar call -h

usage: kevlar call [-A A] [-B B] [-O O] [-E E] [--gen-mask FILE]

                   [--mask-mem MEM] [--mask-max-fpr FPR] [-h] [-d]

                   [--no-homopoly-filter] [--refr FILE] [-o FILE] [-k K]

                   queryseq targetseq

 

Align variant-related reads to the reference genome and call the variant from

the alignment.

 

Required inputs:

  queryseq              contigs assembled by "kevlar assemble"

  targetseq             region of reference genome identified by "kevlar

                        localize"

 

Alignment scoring:

  -A A, --match A       match score; default is 1

  -B B, --mismatch B    mismatch penalty; default is 2

  -O O, --open O        gap open penalty; default is 5

  -E E, --extend E      gap extension penalty; default is 0

 

Mask generation settings:

  --gen-mask FILE       generate a nodetable containing all k-mers that span

                        any variant call

  --mask-mem MEM        memory to allocate for the node table

  --mask-max-fpr FPR    terminate if the estimated false positive rate is

                        higher than "FPR"; default is 0.01

 

Miscellaneous settings:

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

  -d, --debug           show debugging output

  --no-homopoly-filter  by default, short indels adjacent to homopolymers are

                        filtered out; use this flag to disable that filter

  --refr FILE           reference genome indexed for BWA search; if provided,

                        mates of interesting reads will be used to diambiguate

                        multi-mapping contigs

  -o FILE, --out FILE   output file; default is terminal (stdout)

  -k K, --ksize K       k-mer size; default is 31

> kevlar alac -h

# kevlar alac -h

usage: kevlar alac [-p ID] [--max-reads N] [-z Z] [-d D] [-x X]

                   [--include REGEX] [--exclude REGEX] [-A A] [-B B] [-O O]

                   [-E E] [--gen-mask FILE] [--mask-mem MEM]

                   [--mask-max-fpr FPR] [-h] [-o FILE] [-i I] [-k K] [-t T]

                   infile refr

 

Assemble partitioned reads, localize to the reference genome, align the

assembled contig to the reference, and call variants.

 

Required inputs:

  infile                partitioned reads in augmented Fastq format

  refr                  reference genome in Fasta format (indexed for bwa

                        search)

 

Read assembly:

  -p ID, --part-id ID   only process partition "ID" in the input

  --max-reads N         do not attempt to assemble any partitions with more

                        than N reads (default: 10000)

 

Target extraction:

  -z Z, --seed-size Z   seed size; default is 51

  -d D, --delta D       retrieve the genomic interval from the reference by

                        extending beyond the span of all k-mer starting

                        positions by D bp

  -x X, --max-diff X    split and report multiple reference targets if the

                        distance between two seed matches is > X; by default,

                        X is set dynamically for each partition and is equal

                        to 3 times the length of the longest contig in the

                        partition; each resulting bin specifies a reference

                        target sequence to which assembled contigs will

                        subsequently be aligned

  --include REGEX       discard alignments to any chromosomes whose sequence

                        IDs do not match the given pattern

  --exclude REGEX       discard alignments to any chromosomes whose sequence

                        IDs match the given pattern

 

Alignment scoring:

  -A A, --match A       match score; default is 1

  -B B, --mismatch B    mismatch penalty; default is 2

  -O O, --open O        gap open penalty; default is 5

  -E E, --extend E      gap extension penalty; default is 0

 

Mask generation settings:

  --gen-mask FILE       generate a nodetable containing all k-mers that span

                        any variant call

  --mask-mem MEM        memory to allocate for the node table

  --mask-max-fpr FPR    terminate if the estimated false positive rate is

                        higher than "FPR"; default is 0.01

 

Miscellaneous settings:

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

  -o FILE, --out FILE   output file; default is terminal (stdout)

  -i I, --min-ikmers I  do not report calls that a supported by fewer than `I`

                        interesting k-mers

  -k K, --ksize K       k-mer size; default is 31

  -t T, --threads T     process T partitions at a time using T threads

> kevlar unband -h

# kevlar unband -h

usage: kevlar unband [-h] [-n N] [-o FILE] infile [infile ...]

 

When kevlar is run in k-mer banding mode, the same read will typically

appear in multiple output files, annotated with a different set of

potentially novel k-mers in each case. This command will consolidate these

duplicated records across files into a single non-redundant set of reads

with the complete set of novel k-mer annotations.

 

positional arguments:

  infile               input files in augmented Fasta/Fastq format

 

optional arguments:

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

  -n N, --n-batches N  number of batches into which records will be split;

                       default is 16; N temporary files are created and each

                       record from the input is written to a temporary file

                       based on its read name; then each batch is loaded into

                       memory and duplicate records are resolved

  -o FILE, --out FILE  output file; default is terminal (stdout)

> kevlar augment -h

# kevlar augment -h

usage: kevlar augment [-h] [-o FILE] augseqs seqs

 

Internally, kevlar annotates sequences with "interesting k-mers" and uses

"augmented" Fastq and Fasta formats. Processing sequences with third-part

tools usually requires discarding these annotations. This command is used to

augment/reaugment a set of sequences using annotations from an already

augmented sequence file.

 

positional arguments:

  augseqs              augmented sequence file

  seqs                 sequences to annotate

 

optional arguments:

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

  -o FILE, --out FILE  output file; default is terminal (stdout)

> kevlar mutate -h

# kevlar mutate -h

usage: kevlar mutate [-h] [-o FILE] mutations genome

 

Apply the specified mutations to the genome provided.

 

positional arguments:

  mutations            mutations file

  genome               genome to mutate

 

optional arguments:

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

  -o FILE, --out FILE  output file; default is terminal (stdout)

 

 

テストラン

1、両親、子のfastq.gz、リファレンスゲノムをダウンロードする。

curl -L https://s3-us-west-1.amazonaws.com/noble-trios/neon-mother-reads.fq.gz -o mother.fq.gz
curl -L https://s3-us-west-1.amazonaws.com/noble-trios/neon-father-reads.fq.gz -o father.fq.gz
curl -L https://s3-us-west-1.amazonaws.com/noble-trios/neon-proband-reads.fq.gz -o proband.fq.gz
curl -L https://s3-us-west-1.amazonaws.com/noble-trios/neon-refr.fa.gz -o refr.fa.gz

f:id:kazumaxneo:20190409175141j:plain

 

2、リファレンスのindexing

bwa index refr.fa.gz

 

 3、snakemakeのconfigの用意

curl -L https://s3-us-west-1.amazonaws.com/noble-trios/helium-config.json \
| sed "s:/home/user/Desktop:$(pwd):g" > helium-config.json

 

 4、snakemakeのパイプラインの実行

#test run (dry run)
git clone https://github.com/kevlar-dev/kevlar.git
snakemake --configfile kevlar/kevlar/workflows/mark-I/config.json --snakefile kevlar/kevlar/workflows/mark-I/Snakefile --cores 64 --directory dir/ --printshellcmds --dryrun calls

#configファイルを書き換えてから実行する(ncbi/UniVec.faも必要)
snakemake \
--snakefile kevlar/workflows/mark-I/Snakefile \
--configfile helium-config.json --cores 4 --directory workdir -p calls

 

チュートリアルには使用する上での注意点もいくつか書かれています。確認してから使うようにして下さい。

kevlar/kevlar/workflows/mark-I at master · kevlar-dev/kevlar · GitHub

 

引用

Kevlar: a mapping-free framework for accurate discovery of de novo variants

Daniel S. Standage, C. Titus Brown, Fereydoun Hormozdi

bioRxiv preprint first posted online Feb. 13, 2019

 

Kevlar: A Mapping-Free Framework for Accurate Discovery of De Novo Variants
Daniel S. Standage, C. Titus Brown, Fereydoun Hormozdiar

iScience. 2019 Aug 30; 18: 28–36.