macでインフォマティクス

macでインフォマティクス

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

SRA Toolkitのfasta-dumpを高速化した fasterq-dump

2019 4/29 複数ファイルダウンロード例、8/13 ダウンロード例のコード修正、12/18 インストールエラー修正、12/21 実行例追記

2020 1/21 ダウンロード例のコード修正、4/1 リンク追加

2023/07/22 docker イメージ例追加

 

 

 

タイトルの通りのコマンド。 使い方だけ簡単に紹介します。

 

 

 

インストール

sra-toolsを導入すれば使える。

 

downloandリンク

https://www.ncbi.nlm.nih.gov/sra/docs/toolkitsoft/

Github

#bioconda
mamba create -n sratools -y
conda activate sratools
mamba install -c bioconda -y sra-tools

#docker(ayyuan氏のイメージ。おそらく非公式)
docker pull ayyuan/fasterq-dump:1.0

#後日(2019 12/18)別のマシンでテスト時にfasterq-dumpがインストールされなかったので、バイナリをanacondaサーバ(link)からダウンロードした。

#darwin
wget https://anaconda.org/bioconda/sra-tools/2.10.0/download/osx-64/sra-tools-2.10.0-pl526h6de7cb9_0.tar.bz2
#解凍
tar zxvf sra-tools-2.10.0-pl526h6de7cb9_0.tar.bz2
#パスの通ったディレクトリに移動
mv bin/fastq-dump-orig.2.10.0 /usr/local/bin/fasterq-dump

#linux 64bit
wget https://anaconda.org/bioconda/sra-tools/2.10.1/download/linux-64/sra-tools-2.10.1-pl526haddd2b5_0.tar.bz2

fasterq-dump

$ fasterq-dump 

 

Usage:

  fasterq-dump <path> [options]

 

Options:

  -o|--outfile                     output-file 

  -O|--outdir                      output-dir 

  -b|--bufsize                     size of file-buffer dflt=1MB 

  -c|--curcache                    size of cursor-cache dflt=10MB 

  -m|--mem                         memory limit for sorting dflt=100MB 

  -t|--temp                        where to put temp. files dflt=curr dir 

  -e|--threads                     how many thread dflt=6 

  -p|--progress                    show progress 

  -x|--details                     print details 

  -s|--split-spot                  split spots into reads 

  -S|--split-files                 write reads into different files 

  -3|--split-3                     writes single reads in special file 

  --concatenate-reads              writes whole spots into one file 

  -Z|--stdout                      print output to stdout 

  -f|--force                       force to overwrite existing file(s) 

  -N|--rowid-as-name               use row-id as name 

  --skip-technical                 skip technical reads 

  --include-technical              include technical reads 

  -P|--print-read-nr               print read-numbers 

  -M|--min-read-len                filter by sequence-len 

  --table                          which seq-table to use in case of pacbio 

  --strict                         terminate on invalid read 

  -B|--bases                       filter by bases 

  -h|--help                        Output brief explanation for the program. 

  -V|--version                     Display the version of the program then 

                                   quit. 

  -L|--log-level <level>           Logging level as number or enum string. One 

                                   of (fatal|sys|int|err|warn|info|debug) or 

                                   (0-6) Current/default is warn 

  -v|--verbose                     Increase the verbosity of the program 

                                   status messages. Use multiple times for more 

                                   verbosity. Negates quiet. 

  -q|--quiet                       Turn off all status messages for the 

                                   program. Negated by verbose. 

  --option-file <file>             Read more options and parameters from the 

                                   file. 

 

fasterq-dump : 2.9.1 ( 2.9.1-1 )

 

 

 

実行方法

8スレッド指定でSRR000001をダウンロードする。進捗も表示させる。

fasterq-dump SRR000001 -e 8 -p
  • -e    how many thread dflt=6
  • -p    show progress 

カレントにペアエンドfastqが出力される。

 

作業ディレクトリを/tmpに、出力ディレクトリをカレントのoutput/にして実行する。

fasterq-dump SRR000001 -O output -t /tmp -e 8 -p
  • -O   output-dir
  • -t     where to put temp. files dflt=curr dir

$ ls -alth output/

total 593920

-rw-r--r--    1 kazuma  wheel   208M  4 18 00:13 SRR000001_1.fastq

-rw-r--r--    1 kazuma  wheel    75M  4 18 00:13 SRR000001_2.fastq

 

追記

複数SRAをダウンロード(1例)

#以下の3つのペアエンドfastqをダウンロードする
#ヒアドキュメント
cat >sra_ids.txt <<EOF
SRR10712689
SRR10713826
SRR10709253
EOF

#whileで回す。12スレッド並列(*1)。pigzでgzipping(ここでは16スレッド指定)
cat sra_ids.txt | while read line; do
fasterq-dump $line -O ./ -e 12 -p
pigz -p 16 ${line}_1.fastq
pigz -p 16 ${line}_2.fastq
done

 

Githubでは、RAMディスクが利用できる場合、作業ディレクトリをRAMディスクにして高速化する事も提案されています。

引用

Using the SRA Toolkit to convert .sra files into other formats

https://www.ncbi.nlm.nih.gov/books/NBK158900/#SRA_download.what_is_the_purpose_of_the

 

2020  4/1追記

こちらを使ってみて下さい。簡単に個別のデータ、プロジェクト全体のシーケンスデータ、メタデータをダウンロードできます。

 

 

 

 

 

 

 

 

関連


 

 

*1

大規模なデータセットをダウンロードするときはテストしてベストな数値を探ることをお勧めします。 

倍数性レベルを可視化して推測する smudgeplot

2022/09/02 論文引用、Githubリンク修正

 

 性別:それは何の利点があるか?直接的な選択肢が利用可能であるとき、ほとんどの真核生物が繁殖に複雑な迂回路をとる理由は、進化生物学の中心的かつ主として未解決の問題であり続けている。無性生殖を唯一の複製形態として使用する種は系統発生の先端で起こり、成功したのはそのうちのわずか数だけである。言い換えれば、ほとんどの無性系統は最終的に絶滅の運命にあるかもしれない。しかしながら、これらの初期の進化的失敗は、無性種の進化的運命を理解することによって、性の適応的価値について学ぶことができ非常に貴重である。

 多くの研究が無性生殖の動物ゲノムを決定してきたが、それらの多くは、それらを性のある種と区別する特徴を同定することを目的としていた(論文図1)。無性動物では、雌は未受精卵からいわゆるthelytokous parthenogenesis(以下、無性生殖)を介してdaughtersを産む(ref.4)。無性はゲノム進化に多くの結果をもたらすと予測されており、それは減数分裂による配偶子の生産と受精による体細胞倍数性レベルの回復はもはや行われないためである。予測される結果には、例えば有害な変異の蓄積、ヘテロ接合性レベルの変化、転移因子(TE)ダイナミクスなどがある。いくつかの予測は、少数のハウスキーピング遺伝子を用いてゲノムデータなしでテストされている。ハイスループットシークエンシングの到来により、ゲノム規模での無性の古典的予測を評価することができ、さらに以下のような新しい予測を試験することが可能である。パリンドロームの蓄積(論文参照)、これは単一遺伝子アプローチでは研究できなかった。この研究で本著者らは、無性動物に特徴的な任意の重要な特徴を識別できるかどうかを評価するために無性動物種の公表されているゲノムを比較する(論文図1)。 24種は、4種のbdelloid rotifersを含む、これは正統な性別がなくても4000万年以上の間存続し多様化したグループである(ref.12)。Bdelloidsは、これまでのところ予想された無性の行き止まり運命を克服した。メカニズムは絶滅からそれらを保護し、そしてこれらのメカニズムがそれらのゲノムの特定の特徴において見えるかどうか。

(複数段落省略)

 倍数性レベルを推定するためにsmudgeplot v0.1.3(https://github.com/tbenavi1/smudgeplotで入手可能)を使用した。この方法は、互いに1SNPだけ異なるユニークなkmerペアをリードセットから抽出する。これらのkmerペアは、ヘテロ接合性領域に由来すると推測される。そのkmerペアのカバレッジsumはそれらのカバレッジ比率と比較される。この比較は異なるハプロタイプ構造を分離する(論文補足図1b)。最も一般的な構造は、ゲノムの全体的な倍数性を示している。 A. vagaを除くすべての種でこの倍数性推定を使用した。最も一般的な構造は、この種が二倍体であることを示唆していた。 A. vagaは、四倍体として十分に特徴付けられているが、ホモログにあまりに多様性がありすぎて、kmerベースのsmudgeplot法では四倍体を検出することができなかった。
 推定された倍数性レベルを使用して、ゲノムサイズとヘテロ接合性を拡張バージョンのGenomeScopeを使用して推定した。GenomeScopeは、kmerスペクトル分析によってゲノム全体のヘテロ接合性を推定する。入力倍数性を考慮して、近似分布を決定する。推定分布は、ゲノム内で1回、2回など発生するk-merに対応している。次にヘテロ接合性、ゲノム中のリピートの割合、ならびに1nの範囲を推定する。後者はその後ゲノムサイズの推定に使用される。多倍数体のヘテロ接合性の定義は十分に確立されていない(論文のボックス3を参照)が、GenomeScopeは多倍数体の異なる種類のヘテロ接合性遺伝子座を区別する(図3を参照)。
k-merスペクトル分析は、k-mer長の選択によって影響を受ける。より長いk-merはより高いシーケンシングカバレッジを必要とするが、より有益なk-merスペクトルをもたらす。マーブルザリガニを除くすべての種に対してデフォルトの21merサイズを選択した。ここでは、シーケンスの適用範囲が狭いため、17merを選択した。

 

Githubより

 smudgeplotは、(jellyfishやKMCからの)kmerダンプファイルからヘテロ接合のkmerペアを抽出する。 kmerペアカバレッジの合計(CovA + CovB)をそれらの相対的カバレッジ(CovA /(CovA + CovB))と比較することによってゲノム構造を解くことができる。 そのようなアプローチはまた、重複、様々な倍数性レベルなどを用いて不明瞭なゲノムを分析することを可能にする。現在進行中の作業であるが、すでに素晴らしい洞察を提供している。

すべてのハプロタイプ構造はグラフ上でユニークなsmudge(しみ)を持っており、ヒートマップの色が熱くなるしみ部分はハプロタイプ構造が他の構造と比較してゲノムでどのくらい頻繁に表されるかを示している。 GithubのTop画像は理想的なケースである。シーケンスカバレッジがすべてのしみを美しく分離するのに十分であるなら、三倍性の非常に強く明確な証拠を提供する。このツールは近い将来GenomeScopeの一部になる予定である。

 

f:id:kazumaxneo:20190417222139p:plain

Supplementary Figure 1: Genomic evidence of triploidy in M. floridensis .a | the smudgeplot shows dominance of a triploid (AAB) genome structure.

 Preprintより転載。

 

 

インストール

本体 Github

git clone https://github.com/tbenavi1/smudgeplot
cd smudgeplot
pip3 install -r requirements.txt
pip3 install PySAIS==1.0.4

#finally install this package
python3 setup.py install

#conda
mamba install -c bioconda smudgeplot -y

> smudgeplot

# smudgeplot   

usage: smudgeplot <task> [options] 

 

tasks: cutoff    Calculate meaningful values for lower/upper kmer histogram cutoff.

       hetkmers  Calculate unique kmer pairs from a Jellyfish or KMC dump file.

       plot      Generate 2d histogram; infere ploidy and plot a smudgeplot.

       map       Map kmers of individual smudges to a genome.

ERROR:root:"" is not a valid task name

smudgeplot cutoff -h

# smudgeplot cutoff -h

usage: smudgeplot cutoff [-h] infile boundary

 

Calculate meaningful values for lower/upper kmer histogram cutoff.

 

positional arguments:

  infile      Name of the input kmer histogram file (default "kmer.hist")."

  boundary    Which bounary to compute L (lower, default) or U (upper)

 

optional arguments:

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

smudgeplot hetkmers -h

# smudgeplot hetkmers -h

usage: smudgeplot hetkmers [-h] [-o O] [-k K] [-t T] [--all] [infile]

 

Calculate unique kmer pairs from a Jellyfish or KMC dump file.

 

positional arguments:

  infile      Alphabetically sorted Jellyfish or KMC dump file (stdin).

 

optional arguments:

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

  -o O        The pattern used to name the output (kmerpairs).

  -k K        The length of the kmer.

  -t T        Number of processes to use.

  --all       Get all kmer pairs one SNP away from each other (default: just

              the middle one).

smudgeplot plot -h

# smudgeplot plot -h

usage: smudgeplot plot [-h] [-o O] [-q Q] [-L L] [-n N] [-t TITLE] [-m M]

                       [-nbins NBINS] [-kmer_file KMER_FILE] [--homozygous]

                       [infile]

 

Generate 2d histogram for smudgeplot

 

positional arguments:

  infile                name of the input tsv file with covarages (default

                        "coverages_2.tsv")."

 

optional arguments:

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

  -o O                  The pattern used to name the output (smudgeplot).

  -q Q                  Remove kmer pairs with coverage over the specified

                        quantile; (default none).

  -L L                  The lower boundary used when dumping kmers (default

                        min(total_pair_cov) / 2).

  -n N                  The expected haploid coverage (default estimated from

                        data).

  -t TITLE, --title TITLE

                        name printed at the top of the smudgeplot (default

                        none).

  -m M, -method M       The algorithm for annotation of smudges (default

                        'local_aggregation')

  -nbins NBINS          The number of nbins used for smudgeplot matrix (nbins

                        x nbins) (default autodetection).

  -kmer_file KMER_FILE  Name of the input files containing kmer seuqences

                        (assuming the same order as in the coverage file)

  --homozygous          Assume no heterozygosity in the genome - plotting a

                        paralog structure; (default False).

smudgeplot map -h

# smudgeplot map -h

usage: smudgeplot map [-h] [-o O] [-s S] genomefile

 

Identify extracted kmer pairs from individual smudges to a genome.

 

positional arguments:

  genomefile  The reference genome (fasta file; can be gunzipped).

 

optional arguments:

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

  -o O        The pattern used to read kmers and output (default

              "smudgeplot").

  -s S        Smudge to be mapped (intiger, ordered by size; default is to map

              all smudges)

dockerイメージも上げておきます。

docker pull kazumax/smudgeplot

#カレントと/dataをシェアしてラン
docker run -itv $PWD:/data/ kazumax/smudgeplot
> source ~/.profile && ls ~/smudgeplot/
> smudgeplot #help

 

実行方法

1、KMCでk-merヒストグラムを作成する。

#まず作業ディレクトリを作る
kdir tmp
ls *.fastq.gz > FILES
kmc -k21 -t16 -m64 -ci1 -cs10000 @FILES kmer_counts tmp
kmc_tools transform kmer_counts histogram kmer_k21.hist -cx10000
  • -k <len>     k-mer length (k from 1 to 256; default: 25) 
  • -m <size>  max amount of RAM in GB (from 1 to 1024); default: 12
  • -cs <value>   maximal value of a counter (default: 255)

 

2、適切なk-merカバレッジの幅を決定する。smudgeplot/kmer_cov_cutoff.Rを使う。

smudgeplot/kmer_cov_cutoff.R kmer_k21.hist
10 400 #出力

#optional step Rで視覚化
> R #Rに入る
input <- read.table("kmer_k21.hist") #読み込み
plot(input[2:100,],type="l") #2列目が頻度情報

  

3、決定したk-mer頻度レンジ情報を指定し、heterozygous kmersを計算。ここでは30-500とした。

kmc_tools transform kmer_counts -ci30 -cx500 dump -s kmer_k21.dump
smudgeplot hetkmers -k 21 -o kmer_pairs < kmer_k21.dump

 

4、ploidyを推測する2次元プロットを描く。

smudgeplot plot kmer_pairs_coverages.tsv

 

引用

Genomic features of asexual animals
Kamil S. Jaron, Jens Bast, T. Rhyker Ranallo-Benavidez, Marc Robinson-Rechavi, Tanja Schwander

bioRxiv preprint first posted online Dec. 17, 2018

 

GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes

T. Rhyker Ranallo-Benavidez, Kamil S. Jaron & Michael C. Schatz 
Nature Communications volume 11, Article number: 1432 (2020) 

 

関連

 

 

Nanoporeのsquiggle plot (basecallされたリード含む) のシミュレータ DeepSimulator

2019 4/17 誤字修正

2020 1/30タイトル修正

2020 2/1リンク追加

 

 次世代シークエンシング(NGS)技術により、研究者はDNAおよびRNAをハイスループットでシーケンシングすることが可能になり、これはゲノミクス、トランスクリプトミクスおよびエピゲノミクスにおける多数の突破口を開いた(MacLean et al、2009; Metzker、2010;Shi et al、2016; Wu et al、2017)。市場で最も人気のあるNGSテクノロジにはIllumina、PacBio、Nanoporeがある。他のシークエンシング技術とは異なり、ナノポアは、コア成分がナノポアに埋め込まれた電圧バイアス膜を含むポアケミストリーであり、DNAまたはRNA分子が電圧によってポアを通過するように強制されたときに変化する。検出されたシグナルを、Nanopore用に特別に設計されたベースコールに入力すると、ヌクレオチド配列リードを得ることができる。基礎となる設計から恩恵を受け、ナノポアシークエンシングは、ロングリード(Byrne et al、2017)、point-of-care(Lu et al、2016)、およびPCRフリー(Simpsonら、2017)の利点を所有している。これらは、リピート領域のあるゲノムのde novoアセンブリまたはトランスクリプトームアセンブリ、フィールドリアルタイム分析およびダイレクトエピジェネティック検出をそれぞれ可能にする。

 ナノポアシークエンシングの急速な発展と共に、下流のデータ分析方法およびツールもまた急速に出現している。例えば、Graphmap(Sović et al、2016)(紹介)、Minimap2(Li、2017)(紹介)およびMashMap2(Jainら、2017)(MashMap紹介)は、Nanoporeデータをゲノムにマッピングするように設計された。 Canu(Koren et al、2017)(紹介)とRacon(Vaser et al、2017)(紹介)は、Nanoporeによって生成されたロングノイジーリードをまとめるために作成された。近い将来、さらに多くの方法やツールが開発されることが予想される。したがって、経験的(empirical)データ(すなわち実験的に得られたもの)またはシミュレートされたデータ(Escalona et al、2016)のいずれかを用いてこれらの新しい方法をベンチマークすることは非常に重要である。empiricalデータに対して最終的にこの方法を実行することが不可欠であるが、empiricalデータは、未知の根拠のある真実を伴って、取得するのが困難で高価なことがある。それどころか、シミュレートされたデータは低コストで容易に入手することができ、そしてその真実は完全に制御することができる。これらの機能により、シミュレートしたデータを新しい方法のベンチマークの基礎として使用できる。

 NGS技術のために20以上のシミュレータが存在するにもかかわらず(Escalona et al、2016)、Nanoporeシーケンス用に作成されたシミュレータは3つしかない。すなわちReadSim(Lee et al、2014)、SiLiCO(Baker et al、2016)(紹介) )およびNanoSim(Yang et al、2017)(紹介)である。 3つのシミュレータ(論文のセクションS1に示す)の間にはいくつかの違いがあるが、それらは入力ヌクレオチド配列と明示的プロファイル(ここでプロファイルは挿入や欠失などのパラメータのセットを指す)を利用してシミュレーションデータを生成するという同じ特性を共有する。たとえば、ReadSimは固定プロファイルを使用し、SiLiCOはユーザー提供のプロファイルを使用し、NanoSimはシミュレーション段階で使用されるプロファイルを学習するためにユーザー提供のempiricalデータを使用する。しかし、それらのシミュレータは、サンプル調製、電流信号収集、およびベースコールを含む複数の段階を含むナノポアシーケンス手順の複雑な性質を真に捉えていない(論文図1A)。さらに重要なことに、電気信号がナノポアシーケンシングの本質だが、電気信号生成ステップを模倣しようとするようなシミュレータはない。

 統計的観点からシミュレータを設計するという一般的に適応されたシナリオに従う代わりに、本著者らは異なる角度から問題に取り組み、ナノポアシーケンシングのためにより自然に設計された新しいシミュレータを提案する。シミュレータを実行するには、ユーザーは、カバレッジまたはリード数を指定して、リファレンスゲノムまたはアセンブリされたコンティグを入力するだけである。このシーケンスは最初に前処理段階を経て、入力カバレッジの要件と実際のNanoporリードのリード長の分布を満たす、いくつかの短いシーケンスが生成される。次に、それらのシーケンスは、ポアモデル成分とシグナル反復成分を含むシグナル生成モジュールを通過する。ナノポアモデルコンポーネントは、与えられたk-mer(kは通常5または6に等しく、ここでは一般性を失うことなく5-merを使用する)の予測電流信号をモデル化するために使用される。これらのシミュレートされた模擬電流信号は、強度と規模の両方で実際の信号と似ている。最後に、シミュレートされた信号は最終的にAlbacore(https://community.nanoporetech.com/protocols/albacore-offline-basecalli/v/abec_2003_v1_revad_29nov2016/linux)、すなわちオックスフォードナノポアテクノロジー(ONT)の公式ベースコーラーを通過してシミュレートされる。

 本シミュレータのコアコンポーネントが信号生成モジュールのポアモデルであることは明らかである。現在、既存のすべてのポアモデル(https://github.com/nanoporetech/kmer_models)はコンテキストに依存せず、ヌクレオチド配列上の位置に関係なく、予想される電流シグナルに対して5-merごとに固定値を割り当てる。このシミュレータをさらに磨くために、本著者らは、バイオインフォマティクスにおいて大きな可能性を示しているディープラーニングを利用して、新しいコンテキスト依存性のナノポアモデルを提案する(Dai et al、2017; Li et al、2018)。それにもかかわらず電気信号シグナルは通常のヌクレオチド配列よりも8〜10倍長いという事実のため、ディープラーニングモデルを訓練することは簡単ではない。この困難を克服するために、bi-directional long short-term memory(Bi-LSTM)(Graves and Schmidhuber、2005)とdeep canonical time warping (DCTW) (Trigeorgis et al., 2016) を組み合わせた、新規なディープラーニング戦略であるBiLSTM-extended Deep Canonical Time Warping(BDCTW)を提案する。

 上記のように、また図1Bに示されているように、このDeepSimulatorは2つの点で「深く」なっている。まず、結果をまねるだけのシミュレータではなく、処理パイプライン全体をシミュレートすることで、ナノポアシーケンスを深くまねている。次に、シーケンスを電気信号に変換するときに、ディープラーニング法を使用してコンテキスト依存のポアモデルを構築する。 Nanoporeの動作を模倣することにより、本シミュレータは完全なNanoporeシーケンスプロセスをシミュレートし、シミュレートされた電気信号と最終リードの両方を生成する。その上、公式のベースコーラーを使用して、本シミュレータはプロファイルのパラメータを学習する手順を排除するだけでなく、実際に暗黙のうちに実際のパラメータを展開する。さらに、シミュレーション手順を複数のモジュールに分割することによって、本シミュレータはより高い柔軟性を提供する。例えば、ユーザは、異なるベースコーラーを使用すること(Boža et al、2018; Teng et al、2018)を選択するか、または信号生成モジュール内のパラメータを調整して、異なる正確さで最終的なリードを得ることができる。

 

f:id:kazumaxneo:20190416193632p:plain

The Nanopore sequencing procedure.  左が実際のwetの流れで、右がDeepSimulatorの流れ。Githubより転載 

 

 

インストール

ubuntu16.04のMiniconda2.4.0.5環境でテストした。

依存

 Anaconda2 (https://www.anaconda.com/distribution/) or Minoconda2 (https://conda.io/miniconda.html). 

Github

git clone https://github.com/lykaust15/DeepSimulator.git
cd ./DeepSimulator/
./install.sh

> ./deep_simulator.sh

# ./deep_simulator.sh 

DeepSimulator v0.21 [Mar-14-2019] 

    A Deep Learning based Nanopore simulator which can simulate the process of Nanopore sequencing. 

 

USAGE:  ./deep_simulator.sh <-i input_genome> [-n simu_read_num] [-o out_root] [-c CPU_num] [-m sample_mode] [-M simulator] 

                [-C cirular_genome] [-u tune_sampling] [-e event_std] [-f filter_freq] [-s noise_std] [-P perfect] [-H home] 

Options:

 

***** required arguments *****

-i input_genome   : input genome in FASTA format. 

 

***** optional arguments *****

-n simu_read_num  : the number of reads need to be simulated. [default = 100] 

                    Set -1 to simulate the whole input sequence without cut (not suitable for genome-level). 

 

-o out_root       : Default output would the current directory. [default = './${input_name}_DeepSimu'] 

 

-c CPU_num        : Number of processors. [default = 8]

 

-m sample_mode    : choose from the following distribution for the read length. [default = 3] 

                    1: beta_distribution, 2: alpha_distribution, 3: mixed_gamma_dis. 

 

-M simulator      : choose either context-dependent(0) or context-independent(1) simulator. [default = 1] 

 

-C cirular_genome : 0 for linear genome and 1 for circular genome. [default = 0] 

 

-u tune_sampling  : 1 for tuning sampling rate to around eight and 0 for not. [default = 1] 

 

-e event_std      : set the standard deviation (std) of the random noise of the event. [default = 1.0] 

 

-f filter_freq    : set the frequency for the low-pass filter. [default = 850] 

 

-s noise_std      : set the standard deviation (std) of the random noise of the signal. [default = 1.5] 

                    '1.0' would give the base-calling accuracy around 92\%, 

                    '1.5' would give the base-calling accuracy around 90\%, 

                    '2.0' would give the base-calling accuracy around 85\%, 

 

-P perfect        : 0 for normal mode (with length repeat and random noise). [default = 0]

                    1 for perfect context-dependent pore model (without length repeat and random noise). 

                    2 for generating almost perfect reads without any randomness in signals (equal to -e 0 -f 0 -s 0). 

 

-H home           : home directory of DeepSimulator. [default = 'current directory'] 

dockerイメージも置いておきます。

docker pull kazumax/deepsimulator

#カレントと/dataをシェアしてラン
docker run -itv $PWD:/data/ kazumax/deepsimulator
> source ~/.profile && cd ~/DeepSimulator/
> ./deep_simulator.sh #help

 

 テストラン

./deep_simulator.sh -i example/artificial_human_chr22.fasta

artificial_human_chr22_DeepSimu/ が出力される。

f:id:kazumaxneo:20190416231516j:plain

シミュレートされた信号、fast5、basecallingされたfastq、最後にminimap2で評価した精度に関するファイルが出力される。

 

実行方法

fastaを指定して実行する。

./deep_simulator.sh -i example/artificial_human_chr22.fasta

 

出力ディレクトリに生成されるファイルついてはGithubで出力の時系列順に説明されています。詳細はGithubで確認してください。

実行方法

以下のパラメータ設定で50000リード発生させる。

./deep_simulator.sh -i reference.fasta -o output \
-c 16 \
-m 3 \
-M 1 \
-C 1 \
-s 2\
-n 50000
  • -i      input genome in FASTA format.
  • -o     Default output would the current directory.

  • -c     Number of processors. [default = 8]

  • -M    choose either context-dependent(0) or context-independent(1) simulator. [default = 1]

  • -C    cirular_genome : 0 for linear genome and 1 for circular genome. [default = 0]

  • -n    the number of reads need to be simulated. [default = 100]
  • -s    noise_std : set the standard deviation (std) of the random noise of the signal. [default = 1.5]
    '1.0' would give the base-calling accuracy around 92\%,
    '1.5' would give the base-calling accuracy around 90\%,
    '2.0' would give the base-calling accuracy around 85\%,  

重たいので、CPUリソースはあるだけ指定した方が良い。

 

 

オーサーが学習させたケミストリと異なるケミストリをシミュレートさせたい場合、poreモデルもトレーニングもできますが、CPUのみだとかなりの時間を要します。 どうしても行いたい場合、Tensorflow GPU版と対応ライブラリを導入し、グラフィックカードを使って学習時間を短縮させることが推奨されています。

 

追記

引用

DeepSimulator: a deep simulator for Nanopore sequencing
Yu Li, Renmin Han, Chongwei Bi, Mo Li, Sheng Wang, Xin Gao
Bioinformatics, Volume 34, Issue 17, 01 September 2018

 

論文のsupplementary(パラメータ) リンク

 

 

関連


 

 

 

 

 

リファレンスなしでnanopore Direct RNA seqのリードの向きを予測する ReorientExpress

 

 ロングリードシークエンシング技術は、あらゆる種からのトランスクリプトームの体系的な調査を可能にする。ただし、機能評価には5 'から3'への方向を正しく決定する必要がある。 complementary DNA(cDNA)ライブラリーのシーケンシングは、一般に多数のリードをもたらすが(Workmanら、2018)、Oxford Nanopore Technologies(ONT)はRNA分子の直接測定を可能にする(Garalde et al、2018)。strand specificなアダプターを使用できるが、エラー率の多さはアダプター検出を妨げる。現在の方法は、ゲノムまたはトランスクリプトームリファレンスとの比較(Wyman and Mortazavi 2018; Workman et al。2018)、またはさらなる技術の使用(Fu et al、2018)に依存しているため、迅速で費用効果の高いロングリードシーケンシング適用は制限される。ゲノムやトランスクリプトームリファレンスが利用できない種やサンプルでのトランスクリプトームのde-novoの問い合わせを容易にするために、著者らはReorientExpress(https://github.com/comprna/reorientexpress)という、cDNAライブラリーからONTのリードの向きを、アダプターなしで読み取る方法を開発した。 ReorientExpressは、ディープニューラルネットワーク(DNN)を使用して、アダプターとは無関係に、またリファレンスを使用せずにcDNAロングリードの方向を予測する。
  ReorientExpressアプローチは、シークエンシングエラーにもかかわらず、転写産物から得られたリードは依然としてRNA代謝に関与するstrand specificRNPコードに対応する複数のシーケンスモチーフを保持しなければならないという仮説に基づく(Rissland 2017 pubmed; Hentze et al. 2018 pubmed; Gerstberger et al. 2014).。このモデルは、トランスクリプトームアノテーションまたはONTダイレクトRNAシークエンシング(DRS)リードなど、既知の方向を持つ任意のRNAシーケンスセットからトレーニングできる。ReorientExpressは、入力シーケンスから、1から任意の指定された長さまでのkで、正規化されたk-merカウントの行列を計算する(補足Methods)。その後、DNNモデルは、ロングリードを正しい5から3の方向または逆補完の方向に分類する。著者らは、5つの隠れ層を有するDNNを使用し、最初の1つの層には500のノードを配置し、その後は前層の量の半分を使用した(すなわち、500、250、125、62、1)。最後の層はただ1つのノードを持ち、その入力スコアから確率を近似するためにシグモイド関数を適用する。この確率は、リードが正しい方向ではないという確実性を表す。さらに、DNNは過剰適合を減らすためにDropoutのレイヤーを使用する。 ReorientExpressの実装は他のDNNアーキテクチャを可能にする。

(最終段落まで省略)

 近縁種で訓練されたモデルを使用してcDNAナノポアリードの5 'から3'方向の向きを予測する機能により、ReorientExpressは非リードモデル生物からのトランスクリプトームの研究に重要な処理ツールとなる。 

 

ReorientExpressに関するツイート

 

インストール 

ubuntu16.04のminiconda3-4.3.21環境でテストした(docker使用、ホストOS macos10.14)。

本体 Github

pip3 install reorientexpress

reorientexpress.py -h

# reorientexpress.py -h

Using TensorFlow backend.

usage: reorientexpress.py [-h] [-train] [-test] [-predict] -data D -source

                          {annotation,experimental,mapped}

                          [-format {fasta,fastq,auto}] [-annotation A]

                          [-use_all_annotation] [-kmers K] [-reads R]

                          [-trimming T] [-verbose] [-epochs E] [-output O]

                          [-model M]

 

Builds, test and uses models for the orientation of cDNA reads.

 

optional arguments:

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

  -train                Set true to train a model.

  -test                 Set true to test a model.

  -predict              Set true to use a model to make predictions

  -data D, --d D        The path to the input data. Must be either fasta or

                        fastq. Can be compressed in gz.

  -source {annotation,experimental,mapped}, --s {annotation,experimental,mapped}

                        The source of the data. Must be either 'experimental',

                        ' annotation' or 'mapped'. Choose experimental for

                        experiments like RNA-direct, annotation for

                        transcriptomes and mapped for mapped cDNA reads.

                        Mapped reads require a paf file to know the

                        orientation.

  -format {fasta,fastq,auto}, --f {fasta,fastq,auto}

                        The format of the input data. Auto by deafult. Change

                        only if inconsistencies in the name.

  -annotation A, --a A  Path to the paf file if a mapped training set which

                        requires a paf reference is being used.

  -use_all_annotation, -aa

                        Uses all the reads, instead of only keeping

                        antisense,lincRNA,processed_transcript,

                        protein_coding, and retained_intron. Use it also if

                        the fasta has unconventional format and gives errors.

  -kmers K, --k K       The maximum length of the kmers used for training,

                        testing and using the models.

  -reads R, --r R       Number of reads to read from the dataset.

  -trimming T, --t T    Number of nucleotides to trimm at each side. 0 by

                        default.

  -verbose, --v         Whether to print detailed information about the

                        training process.

  -epochs E, --e E      Number of epochs to train the model.

  -output O, --o O      Where to store the outputs. using "--train" outputs a

                        model, while using "-predict" outputs a csv.

                        Corresponding extensions will be added.

  -model M, --m M       The model to test or to predict with.

ro

 

実行方法 

 モデルをトレーニングする。

reorientexpress -train -data path_to_data -source annotation --v -output my_model

 

実行する。

reorientexpress -predict -data path_to_data -source experimental -model path_to_model -output my_predictions

  

引用

ReorientExpress: reference-free orientation of nanopore cDNA reads with deep learning

Angel Ruiz-Reche, Joel A. Indi, Ivan de la Rubia, Eduardo Eyras

bioRxiv preprint first posted online Feb. 18, 2019

rawロングリードから直接MLSTタイピングを行う Krocus

2019 4/16 コマンド修正

 

 急速にコストが下がる中、Pacific Biosciences(PacBio)やOxford Nanopore Technologies(ONT)のようなロングリードシークエンシング技術がアウトブレイク調査に使われ始めている(Faria et al、2017; Quick et al、2015)。そして迅速な臨床診断のために(Votintseva et al、2017)。オックスフォードナノポアのロングリードシーケンサーはほんの数分でシークエンスリードを作成でき、PacBioのシーケンサーは数時間でシークエンスを作成することができる。 MLSTは、細菌を分類するために広く使用されている分類システムである。それはアウトブレイクから単離株を迅速に除外するために使用でき、そしてシーケンスタイプ(ST)を知ることはしばしば細菌の一般的な特徴が推測を可能にする。綿棒から実用的な答えまでの時間を短縮することによって、ゲノミクスは臨床決定に直接影響を及ぼし始めることができ、患者に大きなプラスの影響を与える可能性がある(Gardy&Loman、2018)。

 ロングリードシーケンシング技術によって得られる速度の向上に伴い、ベースエラー率が向上している。ロングリードシーケンシングリードに固有の高いエラー率は、リードを修正するための特殊なツールを必要とする(Koren et al、2017)が、これらの方法はシーケンシングデータを生成するために元の時間より実行に時間がかかるかなりの計算リソース要件をしばしば有する。 

 ショートリードシークエンシング技術のためのMLSTソフトウェアの完全な概説は、Page et al. (2017)から入手可能である。 Page et al. (2017)でレビューされているように、入力としてde novoアセンブリをとる方法だけがロングリードシークエンシング技術では使用できる。しかしながら、de novoアセンブリは、後処理のかなりの計算上のオーバーヘッドを有し、それはシーケンスを実行するのにかかる時間を超える可能性がある。 StringMLST(Gupta、Jordan&Rishishwar、2017)は、k-mer分析を実行することによってraw リードセットからMLSTを迅速に予測するように設計されている。 MentaLiST(Feijao et al、2018)は、同様のk-mer分析アプローチを採用しており、cgMLSTやwgMLSTなどの大規模なタイピングスキーム用に設計されている。それらは、高品質のショートリードシーケンシングデータでのみ機能するように設計されていた。本著者らの知る限り、未修正のロングリードシーケンスデータからMLSTをコールする方法は現在存在しない。

 単離株からの未修正ロングリードから直接STを直接推定することができるKrocusを紹介する。結果は、PacBioシークエンシング技術を使用して生成された700を超えるサンプルを含む細菌のロングリードシーケンシングデータの最大の公開データセット、およびONTデータのリアルデータセットとシミュレートデータセットについて調べられた。平均して、それは未補正のPacbioシーケンシングデータについて94%の感度および97%の特異性で90秒で正しい配列STを生成する。大腸菌リファレンスゲノムに基づく524の模擬ONTサンプルからなるデータセットでは、モデル化したフローセルに応じて感度は82〜94%だった。 12 のKlebsiella pneumoniaeのリアルONTデータセットでは、100%の感度が達成された。 Krocusは、未修正のロングリードから直接MLSTを高精度でコールできる唯一のツールである。これは完全にPython3で書かれており、オープンソースライセンスGNU GPL 3の下でhttps://github.com/andrewjpage/krocusから入手できる。

 

f:id:kazumaxneo:20190416003819p:plain

Figure 1: Flowchart of the Krocus method.  論文より転載

 


インストール

依存

  • Python3

To install Python3 on Ubuntu, as root run:(Githubより)

apt-get update -qq
apt-get install -y git python3 python3-setuptools python3-biopython python3-pip
pip3 install krocus

本体 Github

#anacondaを使っているならcondaで
conda install -c bioconda -y krocus

#pipを使うなら
pip3 install krocus

#or latest development version:
pip3 install git+git://github.com/andrewjpage/krocus.git

krocus -h

# krocus -h

usage: krocus [options] allele_directory input.fastq

 

multi-locus sequence typing (MLST) from uncorrected long reads

 

positional arguments:

  allele_directory      Allele directory

  input_fastq           Input FASTQ file (optionally gzipped)

 

optional arguments:

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

  --filtered_reads_file FILTERED_READS_FILE, -f FILTERED_READS_FILE

                        Filename to save matching reads to (default: None)

  --output_file OUTPUT_FILE, -o OUTPUT_FILE

                        Output file [STDOUT] (default: None)

  --max_gap MAX_GAP     Maximum gap for blocks to be contigous, measured in

                        multiples of the k-mer size (default: 4)

  --margin MARGIN       Flanking region around a block to use for mapping

                        (default: 50)

  --min_block_size MIN_BLOCK_SIZE

                        Minimum block size in bases (default: 150)

  --min_fasta_hits MIN_FASTA_HITS, -m MIN_FASTA_HITS

                        Minimum No. of kmers matching a read (default: 10)

  --min_kmers_for_onex_pass MIN_KMERS_FOR_ONEX_PASS

                        Minimum No. of kmers matching a read in 1st pass

                        (default: 10)

  --max_kmers MAX_KMERS, -r MAX_KMERS

                        Dont count kmers occuring more than this many times

                        (default: 10)

  --print_interval PRINT_INTERVAL, -p PRINT_INTERVAL

                        Print ST every this number of reads (default: 500)

  --kmer KMER, -k KMER  k-mer size (default: 11)

  --divisible_by_3, -d  Genes which are not divisible by 3 are excluded

                        (default: False)

  --target_st TARGET_ST

                        For performance testing print time to find given ST

                        (default: None)

  --verbose, -v         Turn on debugging [False]

  --version             show program's version number and exit

> krocus_database_downloader -h

# krocus_database_downloader -h

usage: krocus_database_downloader [options]

 

Download

 

optional arguments:

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

  --list_species, -l    List all available species (default: False)

  --species SPECIES, -s SPECIES

                        Species to download (default: None)

  --output_directory OUTPUT_DIRECTORY, -o OUTPUT_DIRECTORY

                        Output directory (default: mlst_files)

  --verbose, -v         Turn on debugging (default: False)

  --version             show program's version number and exit

 

データベース

利用できるデータセットを確認する。

krocus_database_downloader -l

# krocus_database_downloader -l

Achromobacter spp.

Acinetobacter baumannii#1

Acinetobacter baumannii#2

Aeromonas spp.

Anaplasma phagocytophilum

Arcobacter spp.

Aspergillus fumigatus

Bacillus cereus

Bacillus licheniformis

Bacillus subtilis

Bartonella bacilliformis

Bartonella henselae

Bordetella spp.

Borrelia spp.

Brachyspira hampsonii

Brachyspira hyodysenteriae

Brachyspira intermedia

Brachyspira pilosicoli

Brachyspira spp.

Brucella spp.

Burkholderia cepacia complex

Burkholderia pseudomallei

Campylobacter concisus/curvus

Campylobacter fetus

Campylobacter helveticus

Campylobacter hyointestinalis

Campylobacter insulaenigrae

Campylobacter jejuni

Campylobacter lanienae

Campylobacter lari

Campylobacter sputorum

Campylobacter upsaliensis

Candida albicans

Candida glabrata

Candida krusei

Candida tropicalis

Candidatus Liberibacter solanacearum

Carnobacterium maltaromaticum

Chlamydiales spp.

Citrobacter freundii

Clonorchis sinensis

Clostridium botulinum

Clostridium difficile

Clostridium septicum

Corynebacterium diphtheriae

Cronobacter spp.

Dichelobacter nodosus

Edwardsiella spp.

Enterobacter cloacae

Enterococcus faecalis

Enterococcus faecium

Escherichia coli#1

Escherichia coli#2

Flavobacterium psychrophilum

Gallibacterium anatis

Haemophilus influenzae

Haemophilus parasuis

Helicobacter cinaedi

Helicobacter pylori

Helicobacter suis

Kingella kingae

Klebsiella aerogenes

Klebsiella oxytoca

Klebsiella pneumoniae

Kudoa septempunctata

Lactobacillus salivarius

Leptospira spp.

Leptospira spp.#2

Leptospira spp.#3

Listeria monocytogenes

Macrococcus canis

Macrococcus caseolyticus

Mannheimia haemolytica

Melissococcus plutonius

Moraxella catarrhalis

Mycobacteria spp.

Mycobacterium abscessus

Mycobacterium massiliense

Mycoplasma agalactiae

Mycoplasma bovis

Mycoplasma hyopneumoniae

Mycoplasma hyorhinis

Mycoplasma iowae

Mycoplasma pneumoniae

Mycoplasma synoviae

Neisseria spp.

Orientia tsutsugamushi

Ornithobacterium rhinotracheale

Paenibacillus larvae

Pasteurella multocida#1

Pasteurella multocida#2

Pediococcus pentosaceus

Photobacterium damselae

Piscirickettsia salmonis

Porphyromonas gingivalis

Propionibacterium acnes

Pseudomonas aeruginosa

Pseudomonas fluorescens

Pseudomonas putida

Rhodococcus spp.

Riemerella anatipestifer

Salmonella enterica

Saprolegnia parasitica

Sinorhizobium spp.

Staphylococcus aureus

Staphylococcus epidermidis

Staphylococcus haemolyticus

Staphylococcus hominis

Staphylococcus lugdunensis

Staphylococcus pseudintermedius

Stenotrophomonas maltophilia

Streptococcus agalactiae

Streptococcus bovis/equinus complex (SBSEC)

Streptococcus canis

Streptococcus dysgalactiae equisimilis

Streptococcus gallolyticus

Streptococcus oralis

Streptococcus pneumoniae

Streptococcus pyogenes

Streptococcus suis

Streptococcus thermophilus

Streptococcus thermophilus#2

Streptococcus uberis

Streptococcus zooepidemicus

Streptomyces spp

Taylorella spp.

Tenacibaculum spp.

Treponema pallidum

Trichomonas vaginalis

Ureaplasma spp.

Vibrio cholerae

Vibrio cholerae#2

Vibrio parahaemolyticus

Vibrio spp.

Vibrio tapetis

Vibrio vulnificus

Wolbachia

Xylella fastidiosa

Yersinia pseudotuberculosis

Yersinia ruckeri

Yersinia spp.

例えばSalmonella entericaのデータセットをダウンロードする。

krocus_database_downloader --species "Salmonella enterica" --output_directory Salmonella_enterica

>ls -alh Salmonella_enterica/

# ls -alh Salmonella_enterica/

total 3.2M

drwxr-xr-x 10 root root  320 Apr 15 15:49 .

drwxr-xr-x  6 root root  192 Apr 15 15:48 ..

-rw-r--r--  1 root root 417K Apr 15 15:48 aroC.tfa

-rw-r--r--  1 root root 415K Apr 15 15:48 dnaN.tfa

-rw-r--r--  1 root root 339K Apr 15 15:48 hemD.tfa

-rw-r--r--  1 root root 582K Apr 15 15:49 hisD.tfa

-rw-r--r--  1 root root 145K Apr 15 15:48 profile.txt

-rw-r--r--  1 root root 348K Apr 15 15:49 purE.tfa

-rw-r--r--  1 root root 429K Apr 15 15:49 sucA.tfa

-rw-r--r--  1 root root 489K Apr 15 15:49 thrA.tfa

root@1fd6

 

実行方法

ダウンロードしたデータベースディレクトリと入力のロングリードfastqを指定し、Krocusを実行する。ここではSalmonella_entericaのディレクトリを指定している。

krocus Salmonella_enterica input.fastq -k 11 -o output
  • -k    k-mer size (default: 11)
  • -o   Output file [STDOUT] (default: None)
  • --filtered_reads_file    Filename to save matching reads to (default: None)

出力の詳細についてはGithubで確認して下さい。

引用

Rapid multi-locus sequence typing direct from uncorrected long reads using Krocus
Andrew J. Page​, Jacqueline A. Keane

PeerJ. 2018 Jul 31;6:e5233

 

 

 

コンソールでbamのカバレッジを素早く確認できる bamcov

 

bamcovは、Florian Breitwieserさん(Github)が公開されている、bamのカバレッジを計算してコンソール上で表示するユーティリティ。

 

インストール

macps10.14でテストした。

Github

git clone --recurse-submodules https://github.com/fbreitwieser/bamcov
cd bamcov
make
make test

> make test

$ make test

./bamcov -H test.bam | column -ts$'\t'

NW_002479957.1  1  514      16    335     65.1751  2.05642   35    6.38

NW_002477246.1  1  2012514  3966  234753  11.6647  0.147286  34.5  43.4

./bamcov -m test.bam

NW_002479957.1 (514bp)

>  90.00% │    ██████████████████    ██████████████          │ Number of reads: 16

>  80.00% │    ██████████████████    ██████████████           

>  70.00% │    ██████████████████    ██████████████          │ Covered bases:   335bp

>  60.00% │    ██████████████████    ██████████████          │ Percent covered: 65.18%

>  50.00% │    ███████████████████  ███████████████          │ Mean coverage:   2.06x

>  40.00% │    ███████████████████  ███████████████          │ Mean baseQ:      35

>  30.00% │    ███████████████████  ███████████████          │ Mean mapQ:       6.38

>  20.00% │   ████████████████████  ███████████████           

>  10.00% │   ████████████████████  ███████████████          │ Histo bin width: 10bp

>   0.00% │   ████████████████████  ███████████████          │ Histo max bin:   100%

          1        100       200       300       400        514   

 

NW_002477246.1 (2.01Mbp)

>  14.45% │                                     █ ▄         │ Number of reads: 3966

>  12.84% │     ▄█          ▄▄      ▄▄    █ ▄  ███      ▄│     (20 filtered)

>  11.24% │▄   ███ █ ▄▄ ▄  ███▄ ██   ▄██  █ ███  ███ ▄██   ██│ Covered bases:   234.8Kbp

>   9.63% │█  ▄███ █ ██▄█  ████ ██ █ ████▄█ ████▄███ ████ ███│ Percent covered: 11.66%

>   8.03% │██ ████ █▄████  ███████▄█▄██████ ████████▄████████│ Mean coverage:   0.147x

>   6.42% │███████ ██████ ███████████████████████████████████│ Mean baseQ:      34.5

>   4.82% │██████████████ ███████████████████████████████████│ Mean mapQ:       43.4

>   3.21% │██████████████▄███████████████████████████████████│ 

>   1.61% │██████████████████████████████████████████████████│ Histo bin width: 40.2Kbp

>   0.00% │██████████████████████████████████████████████████│ Histo max bin:   16.055%

          1       402.5K    805.0K    1.21M     1.61M      2.01M  

 

#help

> ./bamcov -h

$ ./bamcov -h

Usage: bamcov [options] in1.bam [in2.bam [...]]

 

Input options:

  -l, --min-read-len INT  ignore reads shorter than INT bp [0]

  -q, --min-MQ INT        base quality threshold [0]

  -Q, --min-BQ INT        mapping quality threshold [0]

  --rf <int|str>          required flags: skip reads with mask bits unset []

  --ff <int|str>          filter flags: skip reads with mask bits set 

                                      [UNMAP,SECONDARY,QCFAIL,DUP]

Output options:

  -m, --histogram         show histogram instead of tabular output

  -o, --output FILE       write output to FILE [stdout]

  -H, --no-header         don't print a header in tabular mode

  -U, --full-utf          full UTF8 mode for histogram for finer resolution

  -w, --n-bins INT        number of bins in histogram. If 0, use terminal width [50]

  -r, --region REG        show specified region. Format: chr:start-end. 

  -h, --help              help (this page)

 

See manpage for additional details. Tabular format columns::

  rname       Reference name / chromosome

  startpos    Start position

  endpos      End position (or sequence length)

  numreads    Number reads aligned to the region (after filtering)

  covbases    Number of covered bases with depth >= 1

  coverage    Proportion of covered bases [0..1]

  meandepth   Mean depth of coverage

  meanbaseq   Mean baseQ in covered region

  meanmapq    Mean mapQ of selected reads

 

 

実行方法

カバレッジをタブ区切りで表示

bamcov input.bam

-Hをつけるとヘッダー部分(コメント行)は出力しない。

 

カバレッジヒストグラム表示

bamcov -m input.bam

chromosome/plasmidごとに、 chr全体のカバレッジが出力される。

chr (4.5Mbp)

>  90.00% │██████████████████████████████████████████████████│ Number of reads: 9166

>  80.00% │██████████████████████████████████████████████████│     (43 filtered)

>  70.00% │██████████████████████████████████████████████████│ Covered bases:   4.5Mbp

>  60.00% │██████████████████████████████████████████████████│ Percent covered: 99.97%

>  50.00% │██████████████████████████████████████████████████│ Mean coverage:   21.2x

>  40.00% │██████████████████████████████████████████████████│ Mean baseQ:      255

>  30.00% │██████████████████████████████████████████████████│ Mean mapQ:       59.9

>  20.00% │██████████████████████████████████████████████████│ 

>  10.00% │██████████████████████████████████████████████████│ Histo bin width: 71.4Kbp

>   0.00% │██████████████████████████████████████████████████│ Histo max bin:   100%

          1       714.2K    1.43M     2.14M     2.86M      3.57M  

 

chr2の10000-20000のカバレッジだけ表示する。 

bamcov -m -r chr2:10000-20000 input.bam

$ bamcov -m -r chr2:10000-20000 sorted.bam

chr1 (3.4Mbp)

>  90.00% │██████████████████████████████████████████████████│ Number of reads: 35

>  80.00% │██████████████████████████████████████████████████│ 

>  70.00% │██████████████████████████████████████████████████│ Covered bases:   10.0Kbp

>  60.00% │██████████████████████████████████████████████████│ Percent covered: 100%

>  50.00% │██████████████████████████████████████████████████│ Mean coverage:   13.3x

>  40.00% │██████████████████████████████████████████████████│ Mean baseQ:      255

>  30.00% │██████████████████████████████████████████████████│ Mean mapQ:       60

>  20.00% │██████████████████████████████████████████████████│ 

>  10.00% │██████████████████████████████████████████████████│ Histo bin width: 200bp

>   0.00% │██████████████████████████████████████████████████│ Histo max bin:   100%

        10.0K     12.0K     14.0K     16.0K     18.0K      20.0K  

(base)

 

-wで幅を調節できる。-w1で最小幅になり、"-w0"で最大の幅になる。

bamcov -mw2 input.bam

chr (4.5Mbp)

>  89.99% │██│ Number of reads: 9166

>  79.99% │██│     (43 filtered)

>  70.00% │██│ Covered bases:   4.5Mbp

>  60.00% │██│ Percent covered: 99.97%

>  50.00% │██│ Mean coverage:   21.2x

>  40.00% │██│ Mean baseQ:      255

>  30.00% │██│ Mean mapQ:       59.9

>  20.00% │██│ 

>  10.00% │██│ Histo bin width: 1.79Mbp

>   0.00% │██│ Histo max bin:   99.993%

          1         3.57M  

 

plasmid (43.3Kbp)

>  89.83% │█▄│ Number of reads: 242

>  79.85% │██│     (2 filtered)

>  69.87% │██│ Covered bases:   40.1Kbp

>  59.89% │██│ Percent covered: 98.28%

>  49.90% │██│ Mean coverage:   17.8x

>  39.92% │██│ Mean baseQ:      255

>  29.94% │██│ Mean mapQ:       59.9

>  19.96% │██│ 

>   9.98% │██│ Histo bin width: 51.7Kbp

>   0.00% │██│ Histo max bin:   99.808%

          1         103.3K 

 

引用

GitHub - fbreitwieser/bamcov: Quickly calculate and visualize sequence coverage in alignment files

 

関連


 

 

ラージゲノムにもスケールする高速なドラフトゲノム配列polishingツール ntEdit

2019 5/17 論文引用、タイトル修正

2020 10/9 コマンド修正

2021 9/15 インストール手順追加

2022/06/05 condaインストール追記

 

 この10年間で、次世代シーケンシングテクノロジはスループットを大幅に向上させた。例えば、今日では、20 Gbpの針葉樹ゲノムの50倍のカバレッジシーケンシングもIllumina HiSeq-Xマシンなら8レーンフローセル1回で達成できる。しかし、この膨大なデータはバイオインフォマティクスパイプラインにボトルネックを生み出している。典型的には、ショートリードデータでは、未解決の対立遺伝子mixtureを表す偽一倍体ドラフトゲノムは、二倍体配列アセンブリから生じる。使用される方法によっては、これらのアセンブリにかなりの誤差が含まれる可能性がある。ハイスループットシークエンシングプラットフォームのリードのエラー修正のための多くのツールが存在するが[ref.1]、ゲノムアセンブリpolishingツールはほとんど利用可能なものがない。
 アセンブリpolishingの主な用途には、GATK [ref.2]、Pilon [ref.3]、Racon [ref.4]がありる。 PilonとGATKは、ゲノム改良のための確立された包括的なツールであり、短いギャップを埋め、局所的な間違いを修正し、そして変異の塩基を同定し報告する能力を含む。比較すると、Raconはもともと高速ナノポアリード訂正ツールとして設計された、より最近のユーティリティである。後者は、イルミナのデータ、Pacific Biosciences(PacBio)やOxford Nanopore(Nanopore)のシーケンスリードからアセンブリされたものなどの1分子シークエンシング(SMS)ゲノムドラフトを使用して、polishingをかけている[ref. 5]。最新のpolishing精度を考慮すると、Pilonは、微生物および小さな真核生物ゲノム(<100 Mbp)のpolishingに日常的に使用されている堅牢なゲノムアセンブリ改善ツールである。これは人のアセンブリにも適用されている[ref.6]が、残念なことに時間的に2次的に拡大縮小される。
 前述のツールはすべてリードのアライメントを採用している。このパラダイムは、実行時間を犠牲にしても、精査の下で塩基にコンテキストを与える。これらのスケーラビリティの限界に対処するために、本著者らはntEditを開発した。これは、非常に大きなゲノム(> 3Gbp)アセンブリのホモ接合エラーを修正するために長さk(kmer)のワードを使用するユーティリティである。 ntEditは、評価と修正に簡潔なブルームフィルタのデータ構造を採用している。それを他のツールの基本polishing能力、すなわち塩基置換とindelsを修正する能力と比較することによって、ntEditがどのように匹敵する結果を生み出し、そしてヒトの3Gbpゲノムと、トウヒのlarge 20Gbpゲノムに直線的に比例するかを示す。

 

 

インストール

ubuntu16.04でテストした(docker使用。ホストOS macos10.14)。

ビルド依存

c++コンパイラ、zlib、make、autoconf、automakeなど

依存

  1. ntHits (https://github.com/bcgsc/nthits)
  2. BloomFilter utilities (provided in ./lib)
  3. kseq (provided in ./lib)
#ntHitsのビルド
git clone https://github.com/bcgsc/ntHits.git
cd ntHits/
./autogen.sh
./configure
make
make install

nthits --help

# nthits --help

Usage: nthits [OPTION]... FILES...

Reports the most frequent k-mers in FILES(>=1).

Accepatble file formats: fastq, fasta, gz, bz, zip.

 

 Options:

 

  -t, --threads=N use N parallel threads [16]

  -k, --kmer=N the length of kmer [64]

  -c, --cutoff=N the maximum coverage of kmer in output

  -p, --pref=STRING the prefix for output file name [repeat]

  --outbloom output the most frequent k-mers in a Bloom filter

  --solid output the solid k-mers (non-errornous k-mers)

  --help display this help and exit

  --version output version information and exit

 

本体 Github


#nthits
git clone https://github.com/bcgsc/ntHits.git
cd ntHits/
./autogen.sh
./configure
make 
sudo make install 

#ntEdit
git clone https://github.com/bcgsc/ntEdit.git
cd ntEdit
make ntEdit

#conda (未テスト)
mamba create -n ntedit -y
conda activate ntedit
mamba install -c bioconda ntedit nthits -y

> ./ntedit --help

# ./ntedit --help

ntEdit v1.2.0

 

Scalable genome sequence polishing.

 

 Options:

-t, number of threads [default=1]

-f, Draft genome assembly (FASTA, Multi-FASTA, and/or gzipped compatible), REQUIRED

-r, Bloom filter file (generated from ntHits), REQUIRED

-b, output file prefix, OPTIONAL

-k, kmer size, REQUIRED

-z, minimum contig length [default=100]

-i, maximum number of insertion bases to try, range 0-5, [default=4]

-d, maximum number of deletions bases to try, range 0-5, [default=5]

-x, k/x ratio for the number of kmers that should be missing, [default=5.000]

-y, k/y ratio for the number of editted kmers that should be present, [default=9.000]

-c, cap for the number of base insertions that can be made at one position, [default=k*1.5]

-m, mode of editing, range 0-2, [default=0]

0: best substitution, or first good indel

1: best substitution, or best indel

2: best edit overall (suggestion that you reduce i and d for performance)

-v, verbose mode (-v 1 = yes, default = 0, no)

 

--help, display this message and exit 

--version, output version information and exit

 

 

 

実行方法

1、ショートリードを指定する。ここではペアエンドのシーケンシングリードを指定している。カバレッジが30x以上のデータを使うことが推奨されているが (-c2 (>=20X)) 、--c(--cutoff=N)フラグの変更でlow coverageのデータにも対応する(Github参照)。論文ではk=25とk=20が使われている。

nthits --outbloom -p solidBF -k 25 -t 24 pair_R1.fq.gz pair_R2.fq.gz
  • -c    the maximum coverage of kmer in output
  • -k    the length of kmer [64]
  • -t     use N parallel threads [16]
  • -p    the prefix for output file name [repeat]

solidBF_k25.bfが出力される。 

 

2、ドラフトアセンブリ配列を指定してpolishingを実行する。ここでは20スレッド指定。

ntedit -f draft_assembly.fa -r solidBF_k25.bf -b output -t 20
  • -t     number of threads [default=1]
  • -f     Draft genome assembly (FASTA, Multi-FASTA, and/or gzipped compatible), 
  • -k    kmer size, REQUIRED
  • -b    output file prefix, OPTIONAL

ランが終わるとエラーの位置を示したoutput_changes.tsvと、polishingされたoutput_edited.faが出力される。

 

 

ジョブは非常に早く終わる。small data-setとして使用したロングリードのドラフトアセンブリ4Mbpのpolishingは、およそ2秒ほどで終わった(*1)。

引用

ntEdit: scalable genome assembly polishing
Ren ́e L Warren, Lauren Coombe, Hamid Mohamadi, Jessica Zhang, Barry Jaquish, Nathalie Isabel, Steven JM Jones, Jean Bousquet, Joerg Bohlmann, Inanc ̧ Birol

bioRxiv preprint first posted online Mar. 26, 2019

 

ntEdit: scalable genome sequence polishing
René L Warren Lauren Coombe Hamid Mohamadi Jessica Zhang Barry Jaquish Nathalie Isabel Steven J M Jones Jean Bousquet Joerg Bohlmann Inanç Birol
Bioinformatics, btz400, Published: 16 May 2019

 

*1

mac pro (2012) X5690 dualの全CPUリソースを使用。ショートリードはペアエンドfastq1GBx2。

 

関連