ロングリードを使ってSVを検出する Picky - macでインフォマティクス

　ゲノム構造変異の獲得（SV）は、ガンゲノムの主要な特徴であるが、ショートリードシーケンシングデータから再構成することは困難である。ここでは、カスタマイズされたパイプライン、Picky（https://github.com/TheJacksonLaboratory/Picky）を使用し、ナノポアプラットフォームのロングリードを使って乳癌モデルにおける多様なアーキテクチャのSVを明らかにした。著者らはショートリード解析に比べて優れた特異性および感度を有するSVの全スペクトルを同定し、主要な変異源としてリピートDNAを明らかにした。 SVに関連する共通の構造的特徴として、ヌクレオチド解明されたマイクロ挿入におけるゲノムワイドなブレークポイントの検討。ゲノムにわたるブレークポイント密度は、染色体間連結性の傾向と関連しており、ゲノムのプロモーターおよび転写領域が豊富であることが判明した。さらに、染色体二重交叉から段階的SVに至る相互転座の過剰発現を観察した。著者らは、Picky分析が、ロングリードデータからガンゲノムのSVを包括的に検出する有効な手段であることを示す。

Github Motivationsより

　2016年4月時点では、Nanoporeのロングリードに利用できるエンドツーエンドの構造変異検出パイプラインはなかった。著者らはbwa-memのような既存のツールを使って解析を始めたが、機能はするがいくつかの点では最適ではないことに気づいた。例えば、ショートリードの世界では、リードがtandem duplicationを完全にスパンする可能性はほとんど無視できる。従って、大部分の構造変異同定ソフトウェアはtandem duplicationを同定できない。最も重要なのは、SVが識別できるかどうかは、SVコーラーに提供されたアラインメントに依存することである（アライメントが一番大事）。ほとんどのアライナーは低いシークエンシングエラー率を仮定しており、ロングリードでは最適なアライメントを提供しない可能性がある。そこで、ロングリードアラインメントとSVコールの2つの領域を探索することにした。

　SVをコールするには適切なアライメントが重要であるため、bwa-mem、lastal、damapper、graphmap、NextGenMap-LR、blastnを使ってSV（tandem duplication, inversion, deletion, and insertion）を含むロングリードアライメント能力を探索した。

lastal、blastn、そしてdamapperは、非常に似たアライメントパフォーマンス（計算時間パフォーマンスではない）を持つことがわかった。 Bwa-memも近いが、感度の差が無視できなかった。graphmapは、「グローバル」アラインメントのフィットソートを強制的に行う傾向がある。 NextGenMap-LRは有望だが、PacBioのロングリードに焦点を当てていて、ONTのリードでは不安定だったため、使用できなかった。これは徹底的な比較ではないが、実際にアライナーを選択するには役立つ。最終的に、ゲノムとゲノムとのアライメントの文脈で使用されるlastalを、感度と計算時間両方を考慮して使うことに決定した。

　いろいろなSVコーラーを調べると、ほとんどのツールがSVの１つのサブセットに焦点を絞っていて、そこで非常にうまく機能していることが明らかである。また、ショートリードのシナリオで最適化されているため、ロングリードによるSV検出シナリオで問題になることもある。 .vcf形式のレポートは良いスタンダードだが、主にブレークポイントレベルで提示された場合、生物学者が概念レベルでSVを追跡することを難しくしている。著者はコミュニティがより多くのシナリオを想定してスタンダオードの基準を改善し続けることは疑いがない。今のところ、本ツールは都合が良いタブ区切りのテキストファイル形式でレポートし、近い将来に.vcf形式に対応する。

インストール

cent os6でテストした。

依存

LAST

LASTはbrewで導入できる（LAST紹介）。

本体　Github

GitHub - TheJacksonLaboratory/Picky: Structural Variants Pipeline for Long Reads

git clone https://github.com/TheJacksonLaboratory/Picky.git
cd Picky/src/

> perl picky.pl -h

$ perl picky.pl -h

Please specify the command.

picky.pl <command> -h

<command> [hashFq, selectRep, callSV]

hashFq : hash read uuids to friendly ids

lastParam : Last parameters for alignment

selectRep : select representative alignments for read

callSV : call structural variants

xls2vcf : convert Picky sv xls file to vcf

preparepbs: chunk last fastq file and write pbs script for cluster submission

script : write a bash-script for single fastq processing

License = JACKSON LABORATORY NON-COMMERCIAL LICENSE AGREEMENT

https://github.com/TheJacksonLaboratory/Picky/blob/master/LICENSE.md

Documentation = https://github.com/TheJacksonLaboratory/Picky/wiki

Repository = https://github.com/TheJacksonLaboratory/Picky/

ラン

faToFastqを使い、フェイクqualityのfastqを作成する。なければダウンロードする。

curl -O http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faToFastq 
chmod u+x faToFastq

ScrappieでベースコールしたONTのfastaをダウンロードし、fastqにする。

curl -O http://s3.amazonaws.com/nanopore-human-wgs/na12878.chr20ScrappieFiltered.fasta

#フェイクfastqに変換
./faToFastq -qual=H na12878.chr20ScrappieFiltered.fasta LongRead.fastq

hg19のindexを作成する。

lastdb -v -P 10 hg19.lastdb hg19.fa

#sequence dictionaryの作成
samtools dict -H hg19.fa > hg19.seq.dict

以下のコマンドでアライメントとSV検出を実行するためのシェルスクリプトが作成される。

./picky.pl script --fastq LongRead.fastq --thread 20 > LongRead.sh

実行権をつける。

chmod u+x LongRead.sh

中身を確認。

# general installation

export LASTAL=last-755/src/lastal

export PICKY=./picky.pl

# general configuration

export LASTALDB=hg19.lastdb

export LASTALDBFASTA=hg19.fa

export NTHREADS=20

# FASTQ file

export RUN=LongRead

# reads alignments

time (${LASTAL} -C2 -K2 -r1 -q3 -a2 -b1 -v -v -P${NTHREADS} -Q1 ${LASTALDB} ${RUN}.fastq 2>${RUN}.lastal.log | ${PICKY} selectRep --thread ${NTHREADS} --preload 6 1>${RUN}.align 2>${RUN}.selectRep.log)

# call SVs

time (cat ${RUN}.align | ${PICKY} callSV --oprefix ${RUN} --fastq ${RUN}.fastq --genome ${LASTALDBFASTA} --exclude=chrM --sam)

# generate VCF format

${PICKY} xls2vcf --xls ${RUN}.profile.DEL.xls --xls ${RUN}.profile.INS.xls --xls ${RUN}.profile.INDEL.xls --xls ${RUN}.profile.INV.xls --xls ${RUN}.profile.TTLC.xls --xls ${RUN}.profile.TDSR.xls --xls ${RUN}.profile.TDC.xls > ${RUN}.allsv.vcf

hg38でindexを作っているなら、hg19のとこも修正する。

LASTのパスはbrewで通していて間違っていたので、２行目の

export LASTAL=last-755/src/lastal

　↓

export LASTAL=lastal

に修正した。

スクリプトを実行。

./LongRead.sh

ランが終わると VCFが出力される（すでにVCFに対応している）。模擬fastq（合計19億bp）からhg19のアライメントには3時間ほどかかった。

$ head -n 30 LongRead.allsv.vcf |tail -n 10

##FILTER=<ID=singleton,Description="Single read">

##FILTER=<ID=CIPOS,Description="The CIPOS is greater or less than 20">

##FILTER=<ID=CIEND,Description="The CIEND is greater or less than 20">

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT LongRead

chr16 33968730 ID1 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33968756;SVTYPE=DEL;RE=3;RNAMES=nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch253_read1117_strand.fast5,DEAMERNANOPORE_20161206_FNFAB49164_MN16450_sequencing_run_MA_821_R9_4_NA12878_12_06_16_71094_ch489_read1323_strand.fast5,nanopore2_20161128_FNFAB49712_MN17633_sequencing_run_20161128_Human_Qiagen_1D_R9_4_64849_ch247_read294_strand.fast5;SVLEN=27;CIPOS=-2,434;CIEND=-3,743;NOTE=CIPOS_CIEND;ISVTYPE=INDEL(3) GT ./.

chr16 33978912 D2 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33979358;SVTYPE=DEL;RE=2;RNAMES=LomanLabz_PC_20161125_FNFAF01132_MN17024_sequencing_run_20161124_Human_Qiagen_1D_R9_4_96293_ch470_read2257_strand.fast5,PLSP61583_20161015_FNFAB42316_MN17048_sequencing_run_Hum94_ss_jt_86199_ch225_read262_strand1.fast5;SVLEN=447;CIPOS=-90,90;CIEND=-211,212;NOTE=CIPOS_CIEND;ISVTYPE=DEL(1),INDEL(1) GT ./.

chr16 33979001 D3 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33980791;SVTYPE=DEL;RE=9;RNAMES=DEAMERNANOPORE_20161206_FNFAB49164_MN16450_sequencing_run_MA_821_R9_4_NA12878_12_06_16_71094_ch74_read4369_strand.fast5,LomanLabz_PC_20161128_FNFAF01253_MN17024_sequencing_run_20161128_Human_Qiagen_1D_R9_4_50330_ch444_read8778_strand.fast5,PLSP61583_20161128_FNFAB49914_MN17048_sequencing_run_Hu_Nott_Bi_FC4_tune_45519_ch481_read148_strand.fast5,MinION2_20161115_FNFAB46664_MN20093_sequencing_run_Chip105_Genomic_R9_4_450bps_tune_16109_ch48_read301_strand.fast5,FMH_15Le080325s_20161103_FNFAB42798_MN17638_sequencing_run_161103_Human5_LSK108R9_4_13493_ch170_read2217_strand.fast5,nanopore2_20161209_FNFAF04090_MN17633_sequencing_run_20161209_Human_Rapidv2_29324_ch79_read58_strand.fast5,nanopore2_20161122_FNFAF01169_MN17633_sequencing_run_20161122_Human_Qiagen_1D_R9_4_65629_ch138_read425_strand.fast5,MinION2_20161013_FNFAB42706_MN16454_sequencing_run_Chip99_Genomic_R9_4_480bps_69747_ch509_read42_strand.fast5,nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch42_read1895_strand.fast5;SVLEN=1791;CIPOS=-5,11;CIEND=-1,5;NOTE=CIPOS_CIEND;ISVTYPE=DEL(9) GT ./.

chr16 33983845 ID4 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33985021;SVTYPE=DEL;RE=2;RNAMES=MinION2_20161013_FNFAB42706_MN16454_sequencing_run_Chip99_Genomic_R9_4_480bps_69747_ch509_read42_strand.fast5,LomanLabz_PC_20161128_FNFAF01253_MN17024_sequencing_run_20161128_Human_Qiagen_1D_R9_4_50330_ch444_read8778_strand.fast5;SVLEN=1177;CIPOS=-120,121;CIEND=-6,7;NOTE=CIPOS_CIEND;ISVTYPE=INDEL(2) GT ./.

chr16 33985139 ID5 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33985230;SVTYPE=DEL;RE=2;RNAMES=nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch121_read1284_strand.fast5,nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch42_read1895_strand.fast5;SVLEN=92;CIPOS=-21,21;CIEND=-38,38;NOTE=CIPOS_CIEND;ISVTYPE=INDEL(2) GT ./.

chr16 33985190 D6 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33986368;SVTYPE=DEL;RE=2;RNAMES=nanopore2_20161209_FNFAF04090_MN17633_sequencing_run_20161209_Human_Rapidv2_29324_ch79_read58_strand.fast5,nanopore2_20161122_FNFAF01169_MN17633_sequencing_run_20161122_Human_Qiagen_1D_R9_4_65629_ch138_read425_strand.fast5;SVLEN=1179;CIPOS=-2,3;CIEND=-10,11;NOTE=CIPOS_CIEND;ISVTYPE=DEL(2) GT ./.

(anaconda2-4.2.0) [uesaka@cyano src]$

引用

Picky comprehensively detects high-resolution structural variants in nanopore long reads.

Gong L, Wong CH, Cheng WC, Tjong H, Menghi F, Ngan CY, Liu ET, Wei CL

Nat Methods. 2018 Apr 30. doi: 10.1038/s41592-018-0002-6.