macでインフォマティクス

macでインフォマティクス

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

IGVをコマンドラインから起動する

2019/12/28 インストール追記

 

IGVはいくつかのバージョンが提供されいる。デスクトップにショートカットを作ってjavaのwebスタート版を使っている人もいるかもしれないが、コマンドラインから叩くやり方も知っておくと便利である。

 

まずigvを導入。igvtoolsも入れておくとよい。

#igv
conda install -c bioconda -y igv
conda install -c bioconda -y igvtools

#homebrew
brew install igv
brew install igvtools

 

ファンレスファイルを指定してigvを立ち上げる。

fastaを指定。

igv -g input.fa

-g リファレンス配列の指定

f:id:kazumaxneo:20170705133032j:plain

IGVのウィンドウが出現する。

 

genbankファイルを指定してigvを立ち上げる。

igv -g input.gbk

f:id:kazumaxneo:20170705133144j:plain

終了はIGVのウィンドウを閉じるだけでよい。

 

すでにゲノムが登録してある時。

igv -g <name>

すでにIGVに登録済みのゲノムであれば、登録名でも配列の呼び出しはできる。

例えば、下の様に登録済みの7002という配列を呼び出すとする。

f:id:kazumaxneo:20170705134248j:plain

以下の様に打つ。

igv -g 7002

7002の配列を呼び込んだ状態でigvが起動する。

 

ファンレスとそこから作ったbamファイルを指定してIGVを起動する。

igv -g input.fa sample1.bam

 

複数のbamファイルを指定してIGVを起動する。

igv -g input.fa sample1.bam,sample2.bam,sample3.bam

 

複数のbamファイルとgtfファイルを指定して起動する。

igv -g input.fa sample1.bam,sample2.bam,sample3.bam,input.gtf

 

 複数のbamファイルと変異コール結果のVCFファイルを指定して起動する。

igv -g input.fa gatk.vcf,sample1.bam,sample2.bam,sample3.bam,input.gtf

 

さらにgtfとbedファイルも指定して起動する。

igv -g input.fa sample1.bam,sample2.bam,input.gtf,input.bed

f:id:kazumaxneo:20170705163110j:plain

 gtfとbedの両方の表示は、例えばRNAの解析でfusion transcriptのアノテーションファイルを作り、オリジナルのアノテーションと比較したいときなど便利である。

 

注意点

ファイル間にズレがあったり、ヘッダ名が違うとエラーが起こります。こういった時は、genbankfastaなど起点となるファイルを決め、そこから全てのファイルを作り出すとズレが起こりにくくなります。ゲノム解析が進んでいる種なら、もう一度全てのデータを同じデータベースからダウンロードしてみましょう。例えばEnsembl plantだと、各種フォーマットを選んでダウンロードできるようになっています。Ensemblのリンクを貼っておきます。

=> Ensembl ftp top

=> Ensembl Plant

=> Ensembl bacteria

plantはGTFやGFF3もダウンロードできるようです。

 

 

IGVが重い場合は以下を参考にしてください。

 

wiggle形式のcoverage trackを追加する流れは以下にまとめています。


 

RNA seqのワークフロー タキシードプロトコル

 

2012年にnatureでtophatとcufflinksを使ったRNA-seq解析の手法が発表されている。

http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html

tophat、bowtie、cufflinksなどを使ったいわゆるタキシードプロトコルRNA seqの1つのワークフローとして知られており、検索すると多くのプロトコルが見つかる(検索例)。

その1つが以下になる。

上記の手法にのっとり、RNA seq解析を行ってみよう。使えそうなRのプログラムも試してみることにする。

生き物はなんでもよいが、今回は高精度なgtfファイルが手に入るシロイヌナズナRNA seqデータを使うことにする。SRAからアラビドプシスRNA-seqデータを探すと、以下の2群比較のデータが見つかった。これを使ってみる。 

Transcriptome analysis of arabidopsis kamchatica under normal and stress conditions

アラビドプシス属ではあるが、シロイヌナズナではない。Arabidopsis kamchatica(画像)を実験材料に使っている。このシーケンスデータをナズナゲノムに当てて定量していく。

 

シロイヌナズナfastaファイルとgtfファイルはEnsembl genomeの

http://plants.ensembl.org/info/website/ftp/index.html

からダウンロード。

 

nature protocolに載ったワークフローは以下になる。

f:id:kazumaxneo:20170703213533j:plain

 

Prestep、Quality check

 まずデータのクオリティを確認しておく。末端でクオリティが極端に悪くなる部位があればこの段階でクオリティトリムする。

ペアリードのR1を調べる。

fastqc -t 12 -o fastx_test_R1/ -f fastq DRR013381_1.fastq 
  • -t Specifies the number of files which can be processed simultaneously.
  • -o Create all output files in the specified output directory.
  • -f Bypasses the normal sequence file format detection and forces the program to use the specified format. Valid formats are bam,sam,bam_mapped,sam_mapped and fastq
  • -k Specifies the length of Kmer to look for in the Kmer content module. Specified Kmer length must be between 2 and 10. Default length is 7 if not specified.

出力されたhtmlを開く。 

f:id:kazumaxneo:20170630180215j:plain

トリミングは必要なさそうである。マッピングに移ろう。

 

step1、Mapping

cufflinksまでは、biological replicatesがあっても独立したサンプルとして扱っていく。今回はテストなので、chr1だけに限定してマッピングを行う(本当は良くない。アライメントの仕方次第でアーティファクトが出る可能性もゼロではない)。

bowtie2のidnex作成

bowtie2-build -f Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa chr1

最後のchr1はユーザーが決めた名前である。この名前を使いtophatを走らせる。

 

tophat2でpaired-endリードをmapping

tophat2 -p 4 -I 50000 -i 20 --library-type fr-firststrand \
-o tophat_control1 -G Arabidopsis_thaliana.TAIR10.36.gtf chr1 \
DRR013381_1.fastq DRR013381_2.fastq
  • -p  num-threads  
  • -I  --max-intron-length (default: 500000)
  • -i  --min-intron-length  (default: 50)
  • -G --GTF
  • -o --output-dir (default: ./tophat_out)

(gtfをtophatの時点で指定してmappingする方法もある。その場合-Gをつける)。

アライメントファイルが出力される。今回はリファレンスとデータで種が違うので、アライメントは非常に困難になる。エラーは覚悟の上、--read-mismatches 4 --read-gap-length 4 --read-edit-dist 4も上のtohat2につけてランした(How to control the alignment of reads in terms of number of mismatches, gap length etc. ? を参考。edit distはwikiに説明がある)。ナズナはエコタイプの違いでもたくさん塩基置換が起きている。このパラメータでマッピング感度は上がるが、ランダムマッチも増える。どれくらい有効 or 無謀な方法なのかは分からない。または、bowtie2のvery-sensitiveオプションの方がいいかもしれない。tophatはbowtie2のオプションを --b2-のprefixをつけると認識できるようになっている。bowtie2には以下のような一括設定のフラグがある。

--very-fast Same as: -D 5 -R 1 -N 0 -L 22 -i S,0,2.50

--fast Same as: -D 10 -R 2 -N 0 -L 22 -i S,0,2.50

--sensitive Same as: -D 15 -R 2 -L 22 -i S,1,1.15 (default in --end-to-end mode)

--very-sensitive Same as: -D 20 -R 3 -N 0 -L 20 -i S,1,0.50

なのでtophatで--very-sensitive をつけるなら--b2-very-sensitive というフラグをつければいい(リンク)。 

 

 

追記;マッピング率は40%ほどだった(default parameterで20%ほど)。

 

 tophatの出力はaccepted_hits.bamという定型名のため、サンプルの区別がつかなくなる恐れがある。出力されるbamをサンプル名が分かる名前にリネームする。

mv tophat_output/accepted_hits.bam tophat_output/control1.bam

 

のちにIGVで可視化するために、bamに indexをつけておく。

samtools index tophat_output/control1.bam

 gtfファイルを指定して解析した場合、次の2-4は飛ばしてstep5に進む。gtfを指定していない場合、cufflinksにアライメント状況からfeatureファイル(gtf)を作ってもらう必要がある。step2-4を実行する。(詳細については、topにリンクを貼ったnature protocolの論文box1. Cを参照)。

 今回はcuffmergeしてアライメント領域を1から再定義した方が賢明と思われるが、時間があまりない。でナズナのgtfをそのまま使う。

 

2、Cufflinks

アライメント部位の情報に基づき、orfを予測してgtfファイルを出力する。

A、 mRNAを再構成して、gftファイルを出力する。

cufflinks --no-update-check --overlap-radius 1 --library-type fr-firststrand -o cufflinks control1.bam
  • –overlap-radius <int> Transfrags that are separated by less than this distance get merged together, and the gap is filled. Default: 50 (in bp).
  • –no-update-check Turns off the automatic routine that contacts the Cufflinks server to check for a more recent version.

B、出力されるgtfファイルをサンプル名がわかる名前にリネームする。

mv cufflinks/transcripts.gtf cufflinks/control1.gtf 

 これをサンプル全てに対して行う。サンプルが多いならシェルスクリプトを書いて自動処理を考える。

 

3-4、Cuffmerge

このstepもgtfファイルをアライメント時に指定している場合は必要ない。cufflinksに読み込んで全サンプルのgtfを統合し、1つのgtfファイルを出力する。最初に読み込むファイルのリストを作る(サンプルが少ないならコマンド実行時に明示してもよい)。

echo cufflinks/control1.gtf > assemblies.txt
echo cufflinks/control2.gtf >> assemblies.txt
echo cufflinks/control3.gtf >> assemblies.txt
echo cufflinks/stress1.gtf >> assemblies.txt
echo cufflinks/stress2.gtf >> assemblies.txt
echo cufflinks/stress3.gtf >> assemblies.txt

 catで中身を確認しておく。

 

gtfのマージを実行する。

cuffmerge -s Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa  assemblies.txt

merged_asm/merged.gtfができる。このmeggeしたgtfファイルを使って以下の計算を進める。

   

5、Cuffdiff

cuffdiffコマンドで サンプル間の発現量差を計算する。

使い方は

cuffdiff [options] <transcripts.gtf> <sample1_hits.sam> <sample2_hits.sam> [... sampleN_hits.sam]

実際のランは以下のように行った。

cuffdiff -o cuffdiff_output -p 12 --library-type fr-firststrand \
-u Arabidopsis_thaliana.TAIR10.36.gtf -b \
Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa \
-N --compatible-hits-norm -L control,stress \
tophat_control1/control1.bam,tophat_control2/control2.bam,tophat_control3/control3.bam \
tophat_stress1/stress1.bam,tophat_stress2/stress2.bam,tophat_stress3/stress3.bam
  • -o --output-dir
  • -p --num-threads
  • -b use bias correction - reference fasta required
  • -u use 'rescue method' for multi-reads
  • --library-type Library prep used for input reads 
  • -N --compatible-hits-norm count hits compatible with reference RNAs only 
  • -L comma-separated list of condition labels

-Lでつけた名前を今後使うことになる。区別がつくシンプルな名前が良いn > 2の実験でreplicatesがある場合、replicatesは","でつなぐ。異なるサンプルはスペースでつなぐmanual上のデータはcontrolとstressのtriplicateのデータなので、上のboxの様に記述している。そのほか、末端にある""はコマンドを複数行にわたって記述する時のマークである。

 

出力のcuffdiff_outputにはたくさんのファイルができる。

f:id:kazumaxneo:20170704170316j:plain

発現量の差を調べたファイルの中を見てみる。

gene_exp.diff

f:id:kazumaxneo:20170704150718j:plain

最初の10行だけ表示。 value(FPKMではない)とp_valueが表示されている。これだけでも十分解析できそうである。

FPKMはgenes.fplm.trackingの中にある。

 

 

おまけ

igvでbamを開いて、データを確認してみる。igvがない人は

brew install igv

で入る。

fastaとgtf、bamファイルをコマンドから指定して開く。

igv -g Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa 
Arabidopsis_thaliana.TAIR10.36.gtf,tophat_control1/control1.bam,tophat_control2/control2.bam,tophat_control3/control3.bam,tophat_stress1/stress1.bam,tophat_stress2/stress2.bam,tophat_stress3/stress3.bam

gtfファイル以降は ","でつないで記載する(gtf,bam1,bam2,bam3...)。

gtfがソートされてないと怒られるが、そのまま開く。(gtfのソートについてはバイオインフォ道場さんの記事が参考になります http://bioinfo-dojo.net/2017/04/01/igv_gtf_bed/)。

f:id:kazumaxneo:20170704190723j:plain

bedも選択した場合、gtfの下に追加表示される。

 

pop upが邪魔なら上のボタンからnever showをチェック。

f:id:kazumaxneo:20170704190809j:plain

 

 

6-18、CummeRbund

Rにデータを読み込む。ターミナルでRとタイプ

uesaka$ R

 

R version 3.3.0 (2016-05-03) -- "Supposedly Educational"

Copyright (C) 2016 The R Foundation for Statistical Computing

Platform: x86_64-apple-darwin13.4.0 (64-bit)

 

R は、自由なソフトウェアであり、「完全に無保証」です。 

一定の条件に従えば、自由にこれを再配布することができます。 

配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 

 

R は多くの貢献者による共同プロジェクトです。 

詳しくは 'contributors()' と入力してください。 

また、R や R のパッケージを出版物で引用する際の形式については 

'citation()' と入力してください。 

 

'demo()' と入力すればデモをみることができます。 

'help()' とすればオンラインヘルプが出ます。 

'help.start()' で HTML ブラウザによるヘルプがみられます。 

'q()' と入力すれば R を終了します。 

 

ターミナル内でRを起動する。3.3が起動している。

 

 

library("cummeRbund")を打ってcummeRbundを読み込む。cummeRbundを入れてない人はBiocunducterのパッケージをインストールする。

source("http://bioconductor.org/biocLite.R") #時間がかかります 
biocLite("cummeRbund")
biocLite("org.At.tair.db") #GO解析用
library("cummeRbund")
library("org.At.tair.db")

 cuffdiff出力ディレクトリを読み込む。

cuff = readCufflinks('diff_output')

 

cummeRbundのグラフコマンドを一通り実行してみる。

csDensity(genes(cuff))

 

csScatter(genes(cuff), 'stress1', 'control1') 

 

csScatterMatrix(genes(cuff))

 

csVolcanoMatrix(genes(cuff))

 

 

 

 cufflinks2.2.1にバグがあるのか、cummeRbundで正しく読み込めない。というかsqliteのデータのファイルサイズが妙に小さいので、正しくデータベースが構築できていない気がする。以前cuffldiffで解析したデータは正常にcummeRbundに読み込めたため、cummeRbundは問題ない可能性が高い。正しくcummeRbundに読み込め次第追記します。

 

cufflinksを使わず、featurecountsでリードをカウントし、edgeRで正規化、検定してDEGを出す流れは以下の2記事を参考にしてください。

 

リードのカウント

edgeRで検定


全体の流れを再現するならこちらをどうぞ。2017年夏に少人数で勉強会を行った時の資料です。コピペで再現できるはずです。


 

 

引用

 RNA seq workflow

https://en.m.wikibooks.org/wiki/Next_Generation_Sequencing_(NGS)/RNA

 

 

 

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

 

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

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

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

 

6つはいずれも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のアライメントが変化することを検出の拠り所としている。

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

 

依存するもの

 

 

 

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

sh install.sh

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

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と名前が付く以下の6つのファイルが出力される。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.

 

 

トランスポゾン検出ツール5 RelocaTEとRelocaTE2

 

RelocaTE

RelocaTEはゲノム中のトランスポゾンを検出する手法。トランスポゾンの配列を入力してランする。 検出するトランスポゾンの配列、ターゲット配列、などがわかっていないと正しく機能しない。

 

依存するもの

  • Blat
  • Bowtie 1
  • BioPerl
  • SAMtools
  • BWA Recommeded for the creation of the BAM file needed by CharacTErizer
  • Blast (Legacy) formatdb and fastacmd are used for indexed sequence retrieval in an additional companion tool, ConstrucTEr, more info coming soon.

 

Github

GitHub - srobb1/RelocaTE: Find the locations of TEs using the TSD in unassembled short reads by comparing to a closely related reference genome assembly

 

script/relocaTE.pl がTE検出ツールの本体である。

perl relocaTE.pl 

-t  TE FASTA File (Required).

検出したいトランスポゾン配列をfastaで定義する。マルチファスタで複数記入する事もできる。TE配列は末端の terminal inverted repeats (TIRs) [or LTR]も含めて入力してやる必要がある。また、TEのターゲットサイト (TSD) の配列についても塩基で書く必要がある。sampleデータでは">mping TSD=TTA" と3塩基書かれている。GithubのHPには、以下の例が書かれている。

  • TSD=(A|T)GCC => A or Tの後にGC
  • TSD=CGA.A(CT|G) => CGAの後にいずれか1塩基、次がAで、最後がCT or G。

-d Directory of fq files  (Required).

fastqのディレクトリを指定。pairでもpairでなくても入力可能だが、paired-endならばpaired _p1.fq & _p2.fqのような名前にする。拡張子はfq、fastqに対応。

-g Reference genome fasta file (Optional).

任意だが、入力されるとTE挿入位置をコールしてくれる。入力がなければ、TEにアライメントされたリードや、TE末端でトリムされたリードの配列、トリム位置などの情報が出力される。

-e Sample identifier (Optional). 

出力のID名などに使われる。

-o Output directory name (Optional).

デフォルトはoutdir_teSearch。

  -1 Unique mate/pair 1 string (Optional, Recommended). mate-pairのファイル。

  -2 Unique mate/pair 2 string (Optional, Recommended). mate-pairのファイル。

 

他にもいくつかオプションがある。例えば動作速度を上げるためにparallelで並列化するオプションなどは使えそうである。

 

テストランを実行する。

cd sample_relocaTE_run/
sh run_relocaTE.sh #

 

RelocaTE2

RelocaTE2はRelocaTEのバージョンアップ版である。RelocaTEがparallelで並列実行して他のに対し、RelocaTE2はゲノム全体のTEを1サイクルで検索する。これによって、ラージゲノムに数千飛んでいるようなTEのcopy number variationなども現実的な時間で分析できるようになったとされる。

 

Github

GitHub - JinfengChen/RelocaTE2: RelocaTE2

 

 

 

 

 

 

 

 

 

 

バクテリアのIS検出ツール IS_mapper

2019 2/19 インストールの流れを修正 

2021 8/11 condaインストール追記, help更新

 

見つけたいIS配列や抗生物質耐性カセット配列をあらかじめ入力することで、ペアエンド情報を使いISの位置を検出してくれるツール。バクテリア用に設計されており、macbook airなどのlaptopでも高速に動作する。トランスポゾンやマーカー遺伝子でタギングした系の位置の検出にも力を発揮する。論文では予測結果をPCRで確認しており、false callが非常に少ないことも主張されている。

 

原理

BWA-MEMを用いて、ペアエンドのイルミナリードをISクエリ配列にマッピングする 。得られたアラインメントファイル(SAM形式)から、SAMtools view(v0.1.19以降)を用いて、ISクエリ配列の末端にペアがマッピングされていないリード(すなわち、ISの直下に位置する配列を表すリード)を抽出し、SAMフラグに基づいてリードを検索する。具体的には、フラグ「-f 36」を用いて左フランキングリード(入力配列を左から右とする)を、フラグ「-F 40 -f 4」を用いて右フランキングリードを抽出し、別々のBAMファイルに格納した後、BedToolsを用いてFASTQ形式に変換する。Samblaster (v0.1.21) [30] を用いて、SAMファイルから、ISの末端にマッピングされ、隣接する配列にまで及んでいるリード(すなわち「ソフトクリップ」リード、図1b)を抽出する。結果として得られたFASTQファイルを、BioPythonを用いてフィルターをかけ、指定したサイズ範囲(デフォルトでは5~30bp)に収まるリードのうち、ソフトクリップされた部分を抽出する。得られた配列を左右のフランキング配列に分け、それぞれをBWA-MEMを用いて参照ゲノムまたはアセンブリマッピングすることで、解析対象のゲノムにおけるクエリISの位置を特定する。

 

インストール

依存

  • Python v2.7.5
  • BioPython v1.63
  • BWA v0.7.5a
  • Samtools v0.1.19
  • Bedtools v2.20.1 - 
  • BLAST+ v2.2.28 

依存ツールはバージョンアップしているが、最新バージョンでも問題なく動作するようである。全てpipとbrewでインストール可(例: brew install samtools & pip install pysam)。

本体 Github

git clone https://github.com/jhawkey/IS_mapper/
pip install IS_mapper/

#conda
mamba install -c bioconda ismapper -y

ismap --version

$ ismap --version

ismap 2.0

ismap -h

$ ismap -h

usage: ismap [--version] --reads READS [READS ...] --queries QUERIES

             [QUERIES ...] --reference REFERENCE [REFERENCE ...]

             [--output_dir OUTPUT_DIR] [--log LOG] [--help_all HELP_ALL]

 

Basic ISMapper options:

  --version             show program's version number and exit

  --reads READS [READS ...]

                        Paired end reads for analysing (can be gzipped)

  --queries QUERIES [QUERIES ...]

                        Multifasta file for query gene(s) (eg: insertion

                        sequence) that will be mapped to.

  --reference REFERENCE [REFERENCE ...]

                        Reference genome for typing against in genbank format

  --output_dir OUTPUT_DIR

                        Location for all output files (default is current

                        directory).

  --log LOG             Prefix for log file. If not supplied, prefix will be

                        current date and time.

  --help_all HELP_ALL   Display extended help

compiled_table.py -h

$ compiled_table.py -h

usage: compiled_table.py [-h] --tables TABLES [TABLES ...] --reference

                         REFERENCE --query QUERY [--gap GAP] [--cds CDS]

                         [--trna TRNA] [--rrna RRNA] [--imprecise IMPRECISE]

                         [--unconfident UNCONFIDENT] --out_prefix OUT_PREFIX

 

Create a table of IS hits in all isolates for ISMapper

 

optional arguments:

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

  --tables TABLES [TABLES ...]

                        tables to compile

  --reference REFERENCE

                        gbk file of reference

  --query QUERY         fasta file for insertion sequence query for

                        compilation

  --gap GAP             distance between regions to call overlapping, default

                        is 0

  --cds CDS             qualifier containing gene information (default

                        product). Also note that all CDS features MUST have a

                        locus_tag

  --trna TRNA           qualifier containing gene information (default

                        product). Also note that all tRNA features MUST have a

                        locus_tag

  --rrna RRNA           qualifier containing gene information (default

                        product). Also note that all rRNA features MUST have a

                        locus_tag

  --imprecise IMPRECISE

                        Binary value for imprecise (*) hit (can be 1, 0 or

                        0.5), default is 1

  --unconfident UNCONFIDENT

                        Binary value for questionable (?) hit (can be 1, 0 or

                        0.5), default is 0

  --out_prefix OUT_PREFIX

                        Prefix for output file

 

 

テストラン 

オーサーらが準備したテストデータをダウンロードする。

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR225/ERR225612/ERR225612_*.fastq.gz

Githubにテストデータのリファレンス (S_suis_P17.gbk) とIS配列 (ISSsu3.fasta) もアップされている。それもダウンロードする。

 

中身を確認する。

$ cat ISSsu3.fasta

>ISSsu3

TAGTTAAATGAAACAAAAACAGTACATTTATGATATAATGTATTTATGGCATATTCATTA

GATTTTCGTAAAAAAGTTCTCGCATACTGTGAGAAAACCGGCCGTATTACTGAAGCATCA

GCTATTTTCCAAGTTTCACGTAACACTATCTATCAATGGCTAAAATTAAAAGAGAAAACC

GGCGAGCTTCATCACCAAGTTAAAGGAACCAAGCCAAGAAAAGTGGATAGAGATAAATTA

AAGAATTATCTTGAAACTCATCCAGATGCTTATTTGACTGAAATAGCTTCTGAATTTGAC

TGTCATCCAACAGCTATTCATTACGCCCTCAAAGCTATGGGATATACTCGAAAAAAAAGA

GCTGTACCTACTATGAACAAGACCCTGAAAAAGTAGAACTGTTCCTTAAAGAATTGAATA

ACTTAAGCCACTTGACTCCTGTTTATATTGACGAGACAGGGTTTGAGACATGTTTTCATC

GAGAATATGGTCGCTCTTTGAAAGGTCAGTTGATAAAAGGTAAGGTCTCTGGAAGAAGAT

ACCAGCGGATATCTTTAGTGGCAGGTCTCATAAATGGTGCGCTTATAGCCCCGATGACAT

ACAAAGATACTATGACGAGTGGCTTTTTCGAAGCTTGGTTCAAAATATTCTTACTACCCA

CTTTAGGTAAACCATCTGTTATCATCATGGACAATGCAAAGTTTCATAGGATGAGTAAGC

TAAAAGATTTATGCGAGGAGCAGGGACATAGACTTTTACCACTTCCTCCTTACTCACCGG

AATATAATCCCATTGAGAAAATATGGGCTCACATCAAAAAACACCTCAGAAGAGTATTGC

CAAATTGCGATACTTTTCTTGAGGCACTTTCGTCCTGCTCTTGTTTCAGTTGACTA

 

$ head -40 S_suis_P17.gbk

LOCUS       AM946016             2007491 bp    DNA     circular BCT 14-JUL-2009

DEFINITION  Streptococcus suis P1/7 complete genome.

ACCESSION   AM946016

VERSION     AM946016.1  GI:251819067

DBLINK      BioProject: PRJNA352

KEYWORDS    complete genome.

SOURCE      Streptococcus suis P1/7

  ORGANISM  Streptococcus suis P1/7

            Bacteria; Firmicutes; Bacilli; Lactobacillales; Streptococcaceae;

            Streptococcus.

REFERENCE   1

  AUTHORS   Holden,M.T., Hauser,H., Sanders,M., Ngo,T.H., Cherevach,I.,

            Cronin,A., Goodhead,I., Mungall,K., Quail,M.A., Price,C.,

            Rabbinowitsch,E., Sharp,S., Croucher,N.J., Chieu,T.B., Mai,N.T.,

            Diep,T.S., Chinh,N.T., Kehoe,M., Leigh,J.A., Ward,P.N.,

            Dowson,C.G., Whatmore,A.M., Chanter,N., Iversen,P., Gottschalk,M.,

            Slater,J.D., Smith,H.E., Spratt,B.G., Xu,J., Ye,C., Bentley,S.,

            Barrell,B.G., Schultsz,C., Maskell,D.J. and Parkhill,J.

  TITLE     Rapid evolution of virulence and drug resistance in the emerging

            zoonotic pathogen Streptococcus suis

  JOURNAL   PLoS ONE 4 (7), E6072 (2009)

   PUBMED   19603075

  REMARK    Publication Status: Online-Only

REFERENCE   2  (bases 1 to 2007491)

  AUTHORS   Holden,M.T.G.

  TITLE     Direct Submission

  JOURNAL   Submitted (10-MAR-2008) Holden M.T.G., Pathogen Genomics, Sanger

            Institute Wellcome Trust, Wellcome Trust Genome Campus, Hinxton,

            Cambridge, CB10 1SA, UNITED KINGDOM

FEATURES             Location/Qualifiers

     source          1..2007491

                     /organism="Streptococcus suis P1/7"

                     /mol_type="genomic DNA"

                     /strain="P1/7"

                     /db_xref="taxon:218494"

     gene            1..1374

                     /gene="dnaA"

                     /locus_tag="SSU0001"

                     /gene_synonym="dnaH"

     CDS             1..1374

 

準備ができたらランする。gzip圧縮にも対応している。".fastq.gz”になっている必要がある。

ismap --reads ERR225612_*.fastq.gz --queries ISSsu3.fasta --reference S_suis_P17.gbk --output_dir outdir
  • --reads    Paired end reads for analysing (can be gzipped)
  • --queries   Multifasta file for query gene(s) (eg: insertion sequence) that will be mapped to.
  •   --reference    Reference genome for typing against in genbank format
  •   --output_dir   Location for all output files (default is current directory).
  •   --log    Prefix for log file. If not supplied, prefix will be current date and time.

 

 

テストデータの出力

f:id:kazumaxneo:20170703113539j:plain

gbkファイルを入力すると、このようにIS挿入位置の遺伝子をコールしてくれる。

 

  • ISのクエリ配列は--readsを何度も書くことで複数入力できる。
  • リファレンスはfastaかマルチfasta、またはgenbank形式に対応している。テストデータではgenbank形式のリファレンスを用意している。
  • paired-endシーケンスデータはfastqの非圧縮、またはgzip圧縮に対応している。ただし名前はテストのようにペアとわかる構造の名前でなくてはならない。
  • アセンブルしたcontigや、それにアノテーションをかけgbkにしたファイルを使ってもよい。

 

本手法はゲノムにない新規の挿入でも、配列がわかれば挿入位置を拾ってくれる。バクテリアのゲノムだとランニングタイムも数分で快適に使える。トランスタギングしたような株からIS挿入位置を探したいならまずこのツールを試し、それからペアードエンドデータから手動で探す方法と比較してみるとよいと思われる。

 

* 他の手法と同様に、日本語のフォルダパスの下の方にあると動作しないので気をつける。一度エラーになりました。samtoolsは度々仕様が変わりますが、1.4.1でも動作します。

 

引用

ISMapper: identifying transposase insertion sites in bacterial genomes from short read sequence data

Elizabeth Hénaff, Luís Zapata, Josep M. CasacubertaEmail author and Stephan Ossowski

BMC Genomics. 2015 Sep 3;16:667. 

トランスポゾン検出ツール3 Jitterbug

 

ショートリードのアライメントデータから、トランスポゾン挿入位置を検出するツール。入力はリファレンスにアライメントしたbamファイルで、トランスポゾン配列を準備してアライメントする必要はない。配列の位置がgff3で入力されていればよい。その代わりに、Jitterbugはbamファイルのsoft-clipされたリードの情報を使いトランスポゾンを予測する。人とナズナでテストされ、long readで結果は検証されている。

高い感受性を持つとされ、論文中では特定のアレルの挿入がホモかヘテロかも判断できると主張されている。

 

インストール

依存するpytonモジュール

  • python2.7以上
  • pysam (0.7.5 or 0.8.1) 
  • pybedtools
  • psutil
  • matplotlib
  • matplotlib-venn
  • memory_profiler

pipが入っていれば、上記は全てpip install <module名>で導入できる。

本体のダウンロードはGithubから行う。

GIthubリンク

 

pysamの最新版だとcsamtoolsがないというエラーが出る。

Google Code Archive - Long-term storage for Google Code Project Hosting.から、pysam0.7.5をダウンロードし、python2 setup.py installでpysam0.7.5を手動インストールして対応した。

 

実行方法

python jitterbug.py -n 8 -b 50000000 -o /path/to/my/dir/prefix sample.bam te_annot.gff3

 でランする。オプションの後に、bamファイルと、トランスポゾン配列をgff3形式にしたものを入力する。python3以降とpython2.7をいれている人はpython2と明示してランする。

-n Number of CPUs to use

-b If parallelized, size of bins to use, in bp. 

-o Prefix of output files.

-hでヘルプを表示する。

 

data/には2014年にplant journalで発表されたシロイヌナズナゲノムのTAIR10のトランスポゾン遺伝子のgff3ファイルが保存されている。

f:id:kazumaxneo:20170703015328j:plain

27506行ある。

f:id:kazumaxneo:20170703015805j:plain

拡大表示。

ナズナならばこれを使いすぐランできると思われる。

 

 

引用

Jitterbug: somatic and germline transposon insertion detection at single-nucleotide resolution.

Elizabeth Hénaff, Luís Zapata, Josep M. CasacubertaEmail author and Stephan OssowskiEmail author

BMC Genomics 2015 16:768 DOI: 10.1186/s12864-015-1975-5© Hénaff et al. 2015

 

トランスポゾン検出ツール2 ngs_te_mapper

 

 ショートリードをリファレンスゲノムにアライメントし、de novoでトランスポゾン挿入部位を検出する。論文ではBLATをアライメントに使っていたが、gitでダウンロードできる現バージョンはbwaでアライメントを行うようになっている。トランスポゾン挿入時にトランスポゾン末端に複製されてできるtarget site dupplicationを指標にしてトランスポゾン挿入部位を探すらしい。

 

 依存するもの

  • BWA
  • R

本体はGithubからダウンロードできる。

リンク

 

ダウンロードしたディレクトリにテスト用のシェルスクリプトが含まれている。これをランしてみる。

git clone https://github.com/bergmanlab/ngs_te_mapper.git 
cd ngs_te_mapper
bash sourceCode/run_ngs_te_mapper.sh

 上のシェルスクリプトを実行すると、まずwgetでリファンレス配列とシーケンスデータがダウンロードされ、整形後、ngs_te_mapper.Rが実行される。終わると続いてngs_te_logo.Rが実行される。テストだがかなり重たいデータなので、実行時は環境に注意する。

 

実行方法

テスト中

 

 

 

 

引用

Transposon Insertion Finder (TIF): a novel program for detection of de novo transpositions of transposable elements.

Mariko Nakagome, Elena Solovieva, Akira Takahashi, Hiroshi Yasue, Hirohiko Hirochika, and Akio Miyao

BMC Bioinformatics. 2014; 15: 71.