複数のトランスポゾン検出ツールをまとめてインストールして、ランするスクリプト

Githubで公開されているmcclintockは複数のトランスポゾン検出ツールをまとめて走らせることができるツールである。以下の６つのツールを走らせてくれる。

コードはGithubで公開されている。

GitHub - bergmanlab/mcclintock: Meta-pipeline to identify transposable element insertions using next generation sequencing data

６つはいずれもpaired-end情報を用いてトランスポゾン配列（TE）を予測するツールである。ngs_te_mapperは、挿入されるTEの末端の複製配列部分にアライメントされるリード情報からTEを検出する。RelocaTEはTE配列を指定して検出するタイプのツールで、挿入がヘテロかホモかの判定を行うことができる。TEMPはpaired-endリードがTE配列とゲノムにヘテロにマップされることや、paired-endリードのインサートサイズなどからTE挿入部位を検出する（URL）。RetroSeqは、 discorandant mate pairsからTEの挿入位置を予測し、ユーザー定義の情報に従ってTEにアライメントすることでTEを予測するらしい。PoPoolationTEは、TE挿入部位のリードのアライメント状況から向きまで考慮しTEを予測する（論文　Fig. 1）。TE-locateは、paired-endリード情報とTEのデータベースからTEの位置と種類を予測する。いずれも原理は似たようなものであり、TE挿入部分でpaired-endのアライメントが変化することを検出の拠り所としている。

これら６つをまとめて動かすのが本ツール"mcclintock"である。名前はトウモロコシの動く遺伝子Ds/Acを発見したノーベル賞受賞者のBarbara McClintockにそのまま由来していると思われる。公開されているGihubのページにはBarbara McClintockがmacbookを使った写真が公開されている。ユーモラアスな合成写真で中身が頭に入ってこないが、使えれば便利そうである。

依存するもの

FastQC (Command line installation v0.11.2)
RepeatMasker (v.4.0.2) (Necessary if no GFF is supplied)
R (v.3.0.2)
BWA (v.0.7.4-r385)
Perl (v.5.14.2)
BioPerl (v.1.006901)
SAMtools (v.0.1.19-44428cd)
Blat (v.35x1)
Bowtie (v.1.0.0)

まず６つのツールを自動でインストールしてもらう。gitでダウンロードしたフォルダを解凍し、その中で以下のコマンドを打つ。

sh install.sh

６つのツールがインストールされ、依存するものが入ってるかどうかチェックされる。

sed: 1: "RelocaTE/scripts/run_co ...": invalid command code R

sed: 1: "RelocaTE/scripts/sam2ba ...": invalid command code R

sed: 1: "RelocaTE/scripts/sam2fq.pl": invalid command code R

sed: 1: "RelocaTE/scripts/splitS ...": invalid command code R

Testing dependencies...

R... FOUND

RepeatMasker... NOT FOUND

bedtools... FOUND

samtools... FOUND

bcftools... FOUND

bwa... FOUND

CAUTION McClintock requires version 0.7.4-r385 to work correctly with all methods

You currently have Version: 0.7.15-r1140 installed.

exonerate... NOT FOUND

bowtie... FOUND

blat... NOT FOUND

twoBitToFa... NOT FOUND

java... FOUND

perl... FOUND

BioPerl... NOT FOUND

Testing optional dependencies...

fastqc... FOUND

RepeatMasker以外は導入されているようだ。ただしbwaが古いと言われている。

brew upgrade bwaしたが、brewで提供されているbwaの最新版は0.7.15のようだ。手動でインストールする必要があるが、他のツールとの整合性も考えないといけない。今回は気にせずランする。 RepeatMaskerはbrewで素早くインストールした。

テストラン

cd test/
sh runtest.sh

テスト用のイーストのシーケンスデータがダウンロードされ、テストが始まる。

一度readがないというエラーを起こしたが、テストで作られるsacCer2/SRR800842/reads/にwgetされたシーケンスデータのコピーに失敗しているだけだった。自動解凍されたSRR800842_1.fastqとSRR800842_2.fastqをsacCer2/SRR800842/reads/にコピーすると正常に動作した。

ランには少なくとも数時間かかる。終わると、results/originalmethodresults/に各ツールの出力が保存される。さらに、重複を消してフィルタリングもかけられる。最終的にはnonredundantと名前が付く以下の６つのファイルが出力される。ngs_te_mapperだけはフィルター機能がないのでそのままの出力となる。

SRR800842_popoolationte_nonredundant.bed
SRR800842_ngs_te_mapper_nonredundant.bed
SRR800842_relocate_nonredundant.bed
SRR800842_temp_nonredundant.bed
SRR800842_retroseq_nonredundant.bed
SRR800842_telocate_nonredundant.bed

中身を確認する。popoolationte_nonredundant.bed、ngs_te_mapper_nonredundant.bed、retroseq_nonredundant.bedはコール部位がゼロになったので割愛する。

SRR800842_relocate_nonredundant.bed

user$ head /Users/user/local/mcclintock-master/sacCer2/SRR800842/results/SRR800842_relocate_nonredundant.bed

track name="SRR800842_RelocaTE" description="SRR800842_RelocaTE"

chrI 189419 189754 TY2_reference_SRR800842_relocate_sr_1 0 +

chrII 9092 9424 TY2_reference_SRR800842_relocate_sr_2 0 +

chrII 197716 198057 TY3_reference_SRR800842_relocate_sr_3 0 +

chrII 266177 266256 TY1_reference_SRR800842_relocate_sr_4 0 -

chrII 643482 643853 TY4_reference_SRR800842_relocate_sr_5 0 -

chrIII 1178 4322 TY5_reference_SRR800842_relocate_sr_6 0 +

chrIII 124131 124463 TY1_reference_SRR800842_relocate_sr_7 0 -

chrIII 151518 151849 TY1_reference_SRR800842_relocate_sr_8 0 +

chrIII 168386 169064 TY3_reference_SRR800842_relocate_sr_9 0 +

SRR800842_temp_nonredundant.bed

user$ head SRR800842_temp_nonredundant.bed

track name="SRR800842_TEMP" description="SRR800842_TEMP"

chrI 22231 22552 TY1_reference_SRR800842_temp_nonab_1 0 +

chrI 138753 138990 TY1_reference_SRR800842_temp_nonab_2 0 -

chrI 160105 160237 TY1_reference_SRR800842_temp_nonab_3 0 -

chrI 160238 166163 TY1_reference_SRR800842_temp_nonab_4 0 -

chrI 182613 182953 TY3_reference_SRR800842_temp_nonab_5 0 +

chrI 183135 183468 TY1_reference_SRR800842_temp_nonab_6 0 +

chrI 189419 189754 TY2_reference_SRR800842_temp_nonab_7 0 +

chrI 209459 209769 TY1_reference_SRR800842_temp_nonab_8 0 -

chrII 8847 9518 TY1_reference_SRR800842_temp_nonab_9 0 +

telocate_nonredundant.bed

user$ head SRR800842_telocate_nonredundant.bed

track name="SRR800842_TE-locate" description="SRR800842_TE-locate"

2micron 714 715 TY1_non-reference_SRR800842_telocate_rp_1 0 -

2micron 4989 4990 TY1_non-reference_SRR800842_telocate_rp_2 0 +

chrI 3164 3165 TY1_non-reference_SRR800842_telocate_rp_3 0 .

chrI 33938 33939 TY1_non-reference_SRR800842_telocate_rp_4 0 .

chrI 48287 48288 TY4_non-reference_SRR800842_telocate_rp_5 0 -

chrI 55021 55022 TY1_non-reference_SRR800842_telocate_rp_6 0 .

chrI 63043 63044 TY1_non-reference_SRR800842_telocate_rp_7 0 .

chrI 92846 92847 TY1_non-reference_SRR800842_telocate_rp_8 0 .

chrI 100187 100188 TY1_non-reference_SRR800842_telocate_rp_9 0 .

ランするにはmcclintock.shを使う。例えば以下の様なコマンドをうつ。

sh mcclintock.sh -m "RelocaTE TEMP ngs_te_mapper" -r reference.fasta -c te_consensus.fasta -g te_locations.gff -t te_families.tsv -1 sample_1.fastq -2 sample_2.fastq -p 2 -i -b

-m　A string containing the list of software you want the pipeline to use for analysis e.g. "-m relocate TEMP ngs_te_mapper" will launch only those three methods [Optional: default is to run all methods]

-r　A reference genome sequence in fasta format.

-c　The consensus sequences of the TEs for the species in fasta format. [Required]

-g　The locations of known TEs in the reference genome in GFF 3 format. This must include a unique ID attribute for every entry. [Optional]

-t　A tab delimited file with one entry per ID in the GFF file and two columns: the first containing the ID and the second containing the TE family it belongs to. The family should correspond to the names of the sequences in the consensus fasta file. [Optional - required if GFF (option -g) is supplied].

-1　Forward reads. This should be named ending _1.fastq. [Required]

-2　Reverse reads. This should be named ending _2.fastq. [Optional]

-p　The number of processors to use for parallel stages of the pipeline [default = 1]

-i　If this option is specified then all sample specific intermediate files will be removed, leaving only the overall results. The default is to leave sample specific intermediate files (may require large amounts of disk space)

-b　Retain the sorted and indexed BAM file of the paired end data aligned to the reference genome.