macでインフォマティクス

macでインフォマティクス

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

倍数性レベルを可視化して推測する 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) 

 

関連