macでインフォマティクス

macでインフォマティクス

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

Y染色体由来リードをエンリッチする RecoverY

 ハプロイド哺乳動物Y染色体配列は、大規模な次世代配列決定(NGS)プロジェクトではいくつかの理由により適切に組み立てられないことが多い。 Yは女性には存在せず、男性に1コピーのみ存在する。したがって、所望のシーケンスデプスを得るためには、2倍シーケンスを行う必要がある。さらに、常染色体と共有するリピートファミリーの存在、およびX染色体と高い相同性を有する領域の存在は、Yの非固有領域の同定およびアセンブリを複雑にする(Tomaszkiewicz et al、2017で総説)。

 Y染色体アセンブリするいくつかのターゲット技術がある。Single-haplotype iterative mapping and sequencing (SHIMS、例: pubmed)は、ヒトのアセンブリ(Skaletsky et al、2003年)、チンパンジー(Hughes et al、2003年)、アカゲザル(Hughes et al 、2012)およびマウスY染色体(Sohら、2014)を生成するために使用されたBACベースの技術である。 SHIMSは依然として非常に正確な技術だが、ほとんどのプロジェクトではコストと時間がかかりすぎる。あるいは、男性のディープNGSでY染色体配列を含むアセンブリを生成することができる。この場合の課題は、どの種のコンティグがY由来であるかを特定することである。同種のメスのアセンブリに合致しないコンティグ(利用可能な場合)は、Y由来とフラグを付けることができる。アセンブリ前にそのようなアライメントを使用する。しかしながら、そのようなアプローチは、偽常染色体領域のようなXと相同性を有するY配列、またはDAZ遺伝子領域のような常染色体(Saxenaら、2000)を同定することが不十分である。

 男性(XY)対女性(XX)の異なるコピー数のために、オス特異的配列はまた、それらにアライメントするリード数(すなわち、そのデプス)に基づいて同定され得る。種のオスとメスの両方を配列決定し、オスのアセンブリのデプスを比較することは、Yコンティグを同定するのに役立ち得る。このアプローチに基づく方法には、Y染色体ゲノムスキャン(Carvalho and Clark、2013)および染色体商法(Hall et al。、2013)が含まれる。このようなアプローチでは、依然として、男性の高度な深度での配列決定および女性のさらなる配列決定が必要である。さらに、それらは依然として高コピーの転移因子およびXまたは常染色体と共有される他の高コピー数領域を誤って分類する可能性がある。

 配列決定のコストを低減するために、試料中に存在するY染色体物質の量をエンリッチさせるための実験技術が開発されている。染色体特異的濃縮のこのプロセスは、フローソーティング(Dolezel et al、2012)またはレーザーマイクロダイセクション(Zhou and Hu、2007)の実験技術によって達成することができる。特に、フローソーティングは、そのサイズおよびAT / GC比(Dolezel et al、2012)を用いて目的の染色体を分離することができる分離プロセスである。フローソーティングは、ゴリラY(Tomaszkiewiczら、2016)およびブタの性染色体(Skinnerら、2016)の配列決定およびアセンブリの前に適用されている。

  フローソーティングはY染色体のデータ量を増加させるが、試料中に非Y配列が有意な量、依然として残る。これは、大きな染色体由来の破片、またはYと同様のサイズおよびGC含量を有する染色体(例えば、ヒトおよび大型類人猿の染色体21および22)由来であり得る。結果として、シーケンスデータは、フローソーティングの効率によって決まるレベルでのYリードのエンリッチメントを伴うゲノムワイドシーケンスからなる。この非特異性は、アセンブリの前または後にデータから除去しなければならない。これは、オスのシーケンスで上記のディープNGSシーケンスしたように、既知のメスアセンブリに対してリードまたはコンティグをアライメントさせることによって達成できる。しかしながら、すでに述べたように、そのような戦略は女性のリファレンスを必要とし、X染色体および常染色体と高い相同性を有する領域を保持することにおいて非効率的である。

 Y特異的なリードを分離する別の戦略は当初、Tomaszkiewiczらのグループによって提案された (2016)。これは、Yから始まるリードを識別するためk-mer(長さkの部分文字列)カバレッジを使用する戦略である。基本的な原則は、豊富なデータセットの場合、 Y染色体は、ゲノムの他の場所からの非反復性のk-mersよりも高い存在量で生じる。ユーザが選択した閾値を上回る存在量を有するK-merは、Y-mer(Y特異的k-mer)として選択される。続いて、ユーザ定義の閾値(すなわち、Y-mer一致閾値)を上回る多数のY-mersを含むリードは、Yからのものとして識別され、アセンブリを実行するために使用される。この戦略は、フローソートされたデータからゴリラY染色体アセンブルするために使用されたが、時間がかかり、潜在的に偏った推測 - チェック手法を使用して、豊富性とY-merの一致閾値の両方を手動で見つける必要があった。

本論文では、Y染色体特異的リードを単離するための改良された自動化ツールを紹介する。 Tomaszkiewicz et al(2016)で提案された元の戦略に対する主な改善は、存在量閾値の選択とY-mer一致閾値の計算を自動化するアルゴリズムである。これらは、RecoverYのフィルタリングステップの前に実行される。これらの閾値の精度と、ヒトYからシミュレートされたデータを使用したkの選択の効果を評価する。最後に、シミュレートされたデータおよび実際のフローソートされたゴリラYデータに対してRecoverYを実行する。得られたアセンブリは、リードを事前にフィルタリングするか、コンティグをポストフィルタリングするために、メスにアライメントを使用する代替アプローチで得られたアセンブリより優れていることを示す。これらの改良は、Tomaszkiewicz et al(2016)で提案されたパラメータを使用するよりも、ゴリラYでより長い連続したアセンブリにつながる。

 

インストール

mac os 10.13 python2.7.14環境でテストした。

依存

  • numpy
  • biopython
  • matplotlib
  • seaborn

いずれもpipで導入できる。

本体 Github

https://github.com/makovalab-psu/RecoverY

git clone https://github.com/makovalab-psu/RecoverY
cd RecoverY/

> python recoverY.py -h

$ python recoverY.py -h

usage: recoverY.py [-h] [--read_length READ_LENGTH] [--kmer_size KMER_SIZE]

                   [--Ymer_match_threshold YMER_MATCH_THRESHOLD]

                   [--abundance_threshold ABUNDANCE_THRESHOLD]

                   [--threads THREADS] [--plots]

 

RecoverY selects Y-specific reads from an enriched data set.

 

optional arguments:

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

  --read_length READ_LENGTH

                        Set read length (defaults to 150)

  --kmer_size KMER_SIZE

                        Set kmer size (defaults to 25)

  --Ymer_match_threshold YMER_MATCH_THRESHOLD

                        Set Y-mer match threshold (default is calculated by

                        formula : 0.4(l-k+1-2kl/100))

  --abundance_threshold ABUNDANCE_THRESHOLD

                        Set abundance threshold for classifying kmers as

                        Y-mers (default is calculated by using Trusted k-mers)

  --threads THREADS     Set number of threads for RecoverY (defaults to 2)

  --plots               Use this to generate k-mer abundance plots if you have

                        matplotlib & seaborn (defaults to False)

 

 

ラン

付属のテストデータで動作確認する。

cd data/
tar xf kmers_from_reads.tar.xz
cd ../
python recoverY.py # "data"ディレクトリのr1.fastqとr2.fastqを自動認識してランする設計になっている

元データは以下のようになっている。r1とr2がエンリッチされてシーケンスされたシーケンスリードになる。

  • r1.fastq : Enriched raw reads (first in pair)
  • r2.fastq : Enriched raw reads (second in pair)
  • kmers_from_reads : kmer counts for r1.fastq
  • trusted_kmers : kmers occuring in one copy on the Y chromosome. 

Y由来と考えられるリードはoutput/に出力され、op_r1.fastqとop_r2.fastqができる。

 

 

 自分のデータを使ってランするには、dataディレクトリにr1.fastqとr2.fastq、trusted_kmersを配置する。またはRecoverY.pyのコードをその都度修正する。ランにはリードのk-mer配列リストも必要になる。raw fastqのk-merをDSKで分析して、k-merリストを作成しておく。dependencyにDSKのバイナリ(linuxで動作)が入っているので、DSKを動かす。k-merサイズ25なら

cd dependency
./run_dsk_Linux.sh input.fastq 25

 

 kmersファイルが準備できたら、dataディレクトリに配置し、recoverY.pyを実行する。

python recoverY.py --read_length 250 --kmer_size 31 --Ymer_match_threshold 50 --threads 12
  • --kmer_size     kmer size used for classifying reads. This must be the same as DSK's k-mer-size. We recommend a value between 25 and 31 for Illumina 150x150 bp reads. This value is used for calculating Ymer_match_threshold below  (default: 25).
  • --read_length     The length of longest un-trimmed read input to RecoverY. This value is used for calculating Ymer_match_threshold below (default:150) 
  • --Ymer_match_threshold     the number of k-mers a read must match to the Ymer table in order to be classified as coming from the Y. The default value is calculated by the formula : 0.4 * (read_len - kmer_size + 1 - (2*kmer_size*read_len/100)). User may change this, but we recommend a value between 20 and 50 for Illumina 150x150 bp reads.
  • --abundance_threshold    the abundance threshold beyond which all k-mers are assumed to be Ymers. Specify this only if user does not have access to a trusted k-mers file and would like to manually set an abundance threshold.
  • --threads     The number of threads that RecoverY should use (default: 2).

 

引用

RecoverY: k-mer-based read classification for Y-chromosome-specific sequencing and assembly.

Rangavittal S, Harris RS, Cechova M, Tomaszkiewicz M, Chikhi R, Makova KD, Medvedev P

Bioinformatics. 2018 Apr 1;34(7):1125-1131.