IGVのtips 分割表示やGC変化の表示

banファイルを読み込んでリード情報を可視化できるフリーソフト。多様なオプションが用意されており、使いこなすにはかなり勉強が必要。   例を挙げて説明する。 

1、Split view 表示

下のようにウィンドウを分けて表示できる。本来はmate-pairの相方リード表示用だが、複数領域を同時に比較することにも利用できる。   やり方はMiなどのテキストエディタにあらかじめタブで区切った二つのポジションを書いておき、それをウィンドウ内にペーストするだけでOK。

例えばchrの471900-474900と519,700-522,701を同時にみたいなら、間をタブで区切り以下のように記載する。

chr:471900-474900 chr:519,700-522,701  

二つの領域が表示される。以下は４つ表示した例。

split viewを解除するには、上のウィンドウの空間でダブルクリックする。

２、変異のコール結果のVCFファイルをIGVに取り込ませる

 gatkの検出結果を読み込ませてみる。例えば

gatk -T UnifiedGenotyper -R $f -I R1R2_sorted.bam -o UnifiedGenotyper.vcf

出力されるvcfファイルを読み込ませる。

indel検出ソフトの出力はVCF形式に対応しているものが多い。変異をコールして気になる部位があればIGVで見てみればいい。主観なので危険性は伴うが、偽陽性かどうか推測できる場合も多い。

３、ゲノムのGC含量をIGVで表示。

例えば100bpのウィンドウサイズで描画したいとするなら、以下のように行う（1-bpのサイズは意味がないので、ある幅のウィンドウサイズで計算する）

解析に必要なもの

・samtools 

・bed tools

・gawk　（macOSXで計算する時） 

macOSXへのgawkの導入は

brew install gawk

（homebrewがない人は前もって導入しておく。ターミナルでbrewと打ってコマンドが見つからなければ入ってない）

gawkと同様にsamtoolsとbedtoolsの導入は

brew install samtools

brew install bedtools

ゲノム(test.faとする）のGC%を100bpのウィンドウサイズで計算してIGVで可視化するには、以下のようなフローで行なう。

１、fastaのindex作成

samtools test.fa

２、カラム１、２を抽出

cut -f 1,2 test.fasta.fai > test.sizes

３、GC計算幅を指定するファイルを作成

bedtools makewindows -g test.sizes -w 100 > test_100bps.bed

（w=100ならこんなのができる）

chr 0 100

chr 100 200

chr 200 300

chr 300 400

chr 400 500

….

４、ATGC分布を計算。

bedtools nuc -fi test.fasta -bed test_100bps.bed > test_nuc.txt

こんなのができる

head -5 test_100bps.bed

#1_usercol 2_usercol 3_usercol 4_pct_at 5_pct_gc 6_num_A 7_num_C 8_num_G 9_num_T 10_num_N 11_num_oth 12_seq_len

chr 0 100 0.390000 0.610000 20 25 36 19 0 0 100

chr 100 200 0.560000 0.440000 20 18 26 36 0 0 100

chr 200 300 0.580000 0.420000 31 15 27 27 0 0 100

chr 300 400 0.500000 0.500000 23 23 27 27 0 0 100

５、必要なカラムを編集しながら抽出。gawkを使っている。

gawk -v w=100 'BEGIN{FS=" "; OFS=" "}

{if (FNR>1) {print $1,$2,$3,"GCpc_"w"bps",$5}}' test_nuc.txt > test_100bp.igv

成功するとこんなファイルができる。

head -4 test_100bp.igv

chr     0       100     GCpc_100bps     0.610000

chr     100     200     GCpc_100bps     0.440000

chr     200     300     GCpc_100bps     0.420000

chr     300     400     GCpc_100bps     0.500000

６、test_100bp.igvをIGVに読み込ませる。

おまけ；
IGVtoolsでバイナリー化しておけば軽くなる。

IGVtools -> Run igvtools

Command -> toTDF を選択

input file Browse -> .igvファイルを選択。

Runをクリック

test_100bp.igvができる。Load from genomesから読み込ませる。igvtoolsは、この例のようにIGVのGUIから呼び出して計算させることができる。

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

参考ページ

http://wiki.bits.vib.be/index.php/Create_a_GC_content_track

InDel_Hunterのマッピングソフト検討

ARTで250bpでカバレッジ100のシングルfastqを生成。マッピングソフトによるカバレッジの差を調べる。

まずはfastqのジェネレート。

art_illumina -ss MSv3 -sam -i input.fasta -p -l 250 -f 100 -s 10 -o single-read

マッピングソフトデフォルト条件での平均カバレッジ

エラー許容マッピング条件での平均カバレッジ

bwa alnは感度の面から使えないように見える。ただし、リアルデータだと少し話が変わってくることに気づいた。

次はARTのリファレンスにしたゲノムを実際にMiseqで読んだ時のデータ。表にはゲノム全体の平均カバレッジを記載した。

(カバレッジ ＝ アップされたリード数 x 平均リード長 / ゲノムサイズ)

以前読んだシーケンスデータも３つ調べてみた。

 

 bwa alnのマッピング率がARTのデータで極端に落ち込んだ原因は、ARTの fastqが5'側にlow qualityの部位をつくりやすいため。これは MiSeqのデータで5'側のクオリティが落ち込むことを反映しているためらしい。

 ARTの人口データでも、bra alnのシード領域のミスマッチを2まで許容させるとマッピング率は7.1から63.5まで増加した。

変異部位

一番落ち込んでいるのはbwa alnの標準マッピング条件（中央が変異）。エラー許容条件ではindelでもカバレッジが落ち込まない。

最後にmismatch permitting valuの大きさの分析、数値が大きければミスマッチがより大量に蓄積していることを意味する

打ったコマンドは以下の通り。

---------------------------------------------------------------------------------------------------------------

デフォルト条件

１、bwa aln

bwa index -p $f -a is $f

bwa aln -t 20 $f $p > R1.sai

bwa samse -r "@RG\tID:HSS\tSM:WT\tPL:illumina" $f R1.sai $p > R1.sam

2、bwa mem

bwa index -p $f -a is $f

bwa mem -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" -t 22 $f $p $q > bwa_mem_R1.sam

３、smalt

smalt index -k 14 -s 8 k14s8 $f

smalt map -n 22 -o R1.sam k14s8 $p

４、bowtie2

bowtie2-build $f R1_index

bowtie2 -x R1_index -U $p -S R1.sam -k 100 -p 20

---------------------------------------------------------------------------------------------------------------

エラー許容条件

１、bwa aln

bwa index -p $f -a is $f

bwa aln -l 25 -k 0 -n 225 -i 225 -t 20 $f $p > R1.sai

bwa samse -r "@RG\tID:HSS\tSM:WT\tPL:illumina" $f R1.sai $p > R1.sam

２、bowtie2

bowtie2-build $f R1_index

bowtie2 -x R1_index -U $p -L 25 --gbar 225 --mp 2 --end-to-en --np 0 -N 0 -S R1.sam -k 100 -p 24

--------------------------------------------------------------------------------

ゲノム比較ビューア Artemis comparison tool (ACT)

2019 2/13 condaインストール追記

2020 2/25 コメント追加、3/9 インストール方法変更、5/1 使用例追記

2021 1/8 インストール方法変更(blastを追加)、5/23 インストール手順の誤字修正

2023/10/24 biopythonのインストール方法変更

　Artemis comparison tool (ACT)は２つ以上のゲノムを比較して、塩基配列同一性の高い領域を描画するツール。解析には、比較する生物種ごとのfastaファイルと、遺伝子のアノテーションのgenbankフォーマットなどのファイルが必要（遺伝子は任意）。ACTは予めblast検索を行っておき、それからjavaのプログラムを立ち上げる。そのため、ACTを起動する前に、blastサーバーにアクセスして総当たりblastするか、ローカルblastコマンドで総当たりのblast検索を行う必要がある（具体的にはローカルブラストを -outfmt 6をつけて実行し、タブ形式の結果ファイルを出力する）。

 下記に、blastサーバーを利用したやり方とローカルblastを行う方法をまとめる。 

Document

インストール

#bioconda
conda install artemis -c bioconda




#2020 3/9 condaでpython2の環境を作ってbwast.pyを動かせるようにしておくと楽にゲノム比較できます。
mamba create -n artemis -y python=2.7
conda activate artemis
mamba install -c bioconda blast -y
mamba install -c bioconda artemis -y
mamba install -c anaconda biopython -y

#bwastの導入
git clone https://github.com/bawee/bwast.git

#genbank
python bwast/bwast.py genome1.gbk genome2.gbk genome3.gbk -a
#fasta
python bwast/bwast.py genome1.fa genome2.fa genome3.fa -a

 手順 

１-B、blastサーバーを使ったblast解析を行う場合（bwastを使うなら不要）

ACTの解析専用に、クラウドでブラスト解析を行ってくれるサーバーが下記URLである。

http://www.hpa-bioinfotools.org.uk/pise/double_actv2.html

top画面

2020 2/25追記

既にサーバーは停止しているようです。ローカルマシンでblastn検索を行ってください。

２つある大きなウィンドウに比較するゲノムのfasta形式のファイルを２つ入れる。テストなら、OR select buttonを押して生物種を選ぶ。下の方のblastnかblastpを選び、メールアドレスを記入してRun_genome_blastボタンを押せば解析が始まる。blastは精度より速度重視の手法のため、プラスミドのような小さいゲノムだと瞬時に終わり下の画面になる。

ブラスト終了後の画面

上のgenome_blast.result を右クリックして "リンク先のファイルをダウンロード"でダウンロードする。中身はこのようなファイル。

user$ head -4 genome_blast.result.txt 

39750 99.00 3214 24043 M0_2_50000_blastdb_WT_pilon.fasta.00000001.out 3214 24042 1 dna:chromosome chromosome:ASM83294v1:1:1:4109373:1 REF

20400 99.00 39426 50000 M0_2_50000_blastdb_WT_pilon.fasta.00000001.out 39408 49980 1 dna:chromosome chromosome:ASM83294v1:1:1:4109373:1 REF

1219 97.00 31977 32642 M0_2_50000_blastdb_WT_pilon.fasta.00000001.out 31976 32641 1 dna:chromosome chromosome:ASM83294v1:1:1:4109373:1 REF

1154 100.00 36841 37422 M0_2_50000_blastdb_WT_pilon.fasta.00000001.out 36840 37421 1 dna:chromosome chromosome:ASM83294v1:1:1:4109373:1 REF

１-B、local blastを行う場合（blastn）

makeblastdb -in reference.fa -out blastDB -dbtype nucl #データベース作成
blastn -db blastDB -query input.fasta -out blast_result.txt -outfmt 6

 比較元をデータベースにし、その比較対象をクエリにしてblast解析を行う。この出力をACTのcomparison fileに指定する。

次にACTのダウンロードを行う。 

２、ACTのインストール と起動

http://www.sanger.ac.uk/science/tools/artemis-comparison-tool-actにアクセスする。

上の画面中央付近の"Mac OS"をクリックしてダウンロードする。自己解凍されてできた４つのアプリをApplications/にコピーする。これでインストールは完了。

ただし最新のmac OSだと起動できない。起動できない人は、上のリンクから、java web start版をダウンロードする。 ダウンロードしたアイコンを叩くと起動するが、mac OS sierraではApp store以外から落としたファイルを起動できない。 怒られたら、リンゴマークから環境設定を起動し、セキュリティとプライバシの項目のこのまま開くをクリック、を選択する。   

パスを入れて起動する。それでも怒られる場合、jabaのセキュリティに例外URLを登録する必要がある。リンゴマーク -> 環境設定 -> javaを選択。

java-> セキュリティタブをクリック。サイトリストの編集をクリック

http://www.genedb.org を追加して保存。もう一度ダウンロードしたファイルを叩くと、また文句を言われるが起動できるようになる。

３、gbkファイルの読み込み

起動したら、File -> Open ...を選択。

下のウィンドウが出てくるので、ブラスト解析に使ったデータのgenbankファイル２つを選ぶ。comparison fileにはブラスト解析でダウンロードしたテキストファイルを選択。

ゲノムが似たもの同士を比較すると、初期条件ではこのように蜘蛛の巣を張ったような感じになる。

閾値を上げて完全マッチだけ取り出す。画面で右クリックしてSet Score Cutoffを選択

Minimum cutoffゲージを少し上げる。

すると下のように変化する。上のウィンドウを閉じると元に戻る。 

比較には、ゲノムの中央で大きなインバーション（逆位）が起きている２ゲノムを使っている。上の図を見るとそれが一目瞭然である。

次のデータは、ほぼ同じゲノム構造を持っているが、ざっくり1kbpに１つくらい塩基置換やindelが起きている生物種同士の比較（E.coliのK12とO157のような同じ種間の生物）。

端からみていくと、例えば下の領域には片方の生物しかない遺伝子が見つかる。

 次は生物としてはあまり近縁でないが、特殊な代謝系遺伝子セットをどちらの生物も持っている例。同調して機能する特定の遺伝子セット自体はまとまって存在しているが、下の図を見るとその部分だけ線が描かれ相同性が高くなっていることが分かる。

 もっと例を見たけば、Slide shareのリンクの108-109スライド目あたりが参考になる。 http://www.slideshare.net/leightonp/plant-pathogen-genome-data-my-life-in-sequences

実はACTはサポートするラッパーツールと組み合わせることで劇的に使いやすくなる。下のエントリーでそのことを紹介しています。ラッパーツールを入れて、blastからACTの比較までをほぼ自動化するやり方を紹介してますので、興味があればご覧ください。

2020 5/1追記

artemisはリピート配列がゲノム上にどのくらい広がっているか視覚化するにも役立つ。その場合は自身のゲノムと比較する。上のラッパーツールを使うなら

conda activate artemis
python bwast/bwast.py genome.fasta genome.fasta -a

 起動後、セルフマッチを非表示にする。

repetitive region由来と思われるマッチが見えるようになった。

備考: http://www.webact.org/WebACT/generateというACTの解析を自動化するサービスもあるが、こちらはうまくワークしないので検証できなかった。開設から時間が経っているのでもしかするとサーバーがゾンビ化してるのかもしれない。

参考 ---------------------------------------------------------------------------------------------------------------------------

１、Artemis-ACT instructions for NOVA cours　　

http://www.pseudomonas-syringae.org/Artemis-ACT-NOVA.html  

２、ACT User manual ftp://ftp.sanger.ac.uk/pub/project/pathogens/workshops/BioLinux_Artemis_Data/Artemis-BioLinux.pdf    

引用

ACT: the Artemis Comparison Tool.
Carver TJ, Rutherford KM, Berriman M, Rajandream MA, Barrell BG, Parkhill J

Bioinformatics. 2005 Aug 15;21(16):3422-3. Epub 2005 Jun 23.

追記

関連

Artemis使用法

2020 3/1 リンク追加

Artemisはサンガー研で開発されたゲノムデータをビジュアルで見れるソフト。ゲノムのデータ（genebakファイルなど）を読み込ませると、遺伝子単位まで拡大して確認できる。ダウンロードはここから http://www.sanger.ac.uk/science/tools/artemis

Macなら、写真のMac OSをクリック。自動的にダウンロードが開始される。ダウンロードした.dmgファイルを解凍するとArtemisが表示される。Applicationにコピーしてインストール完了。 起動したらFile -> open... を選択し、

genebakフォーマットのファイル(.gbk)を選択する。

読み込むと以下のようにorfを可視化して表示できる。

拡大縮小方法はやや特殊。

上の右端に小さな灰色ゲージが見えるが、これをマウスで上下に引っ張れば拡大縮小される。ゲノム全体を表示することもできる。

 まとめ直しました。

真核生物のRNA-seqデータ解析

今回はhumanのデータを使う。重いのでchromosome19に限定して解析する。

------------準備------------

1、データ

fastaファイル。 FTP Downloadからhumanのfastaファイルをダウンロードする。

Homo_sapiens.GRCh38.dna_rm.chromosome.19.fa

rRNAのfastaファイル。ここではgitに上がっているのをダウンロード。

IlluminaPE/humRibosomal.fa at master · Magdoll/IlluminaPE · GitHub

配列をコピペしてrRNA.faとして保存。

gtfファイル。 FTP Downloadからhumanのgtfファイルをダウンロードする。

Homo_sapiens.GRCh38.87.gtf

fastqファイル 。http://trace.ddbj.nig.ac.jp/DRASearch/experiment?acc=DRX015048に行って、右端のfastqを選択。ウィンドウを開くか聞いてくるのでゲストを選択して開く。

DRR016695_1.fastq.bz2とDRR016695_2.fastq.bz2を解析フォルダにコピペしてダウンロード（bz2圧縮で2GBx2）。

２、解析プログラム

homebrewで実行環境を整える。homebrewはこちら

brew install tophat

brew install bowtie2

brew install samtools

brew install cufflinks

------------解析------------

はじめにrRNAにmapされたリードを除く。それからgtfに記載されている既存exsonにマッピングして、fpkmを算出する。

１、rRNA.faのindex作成 (bowit2を使う）

bowtie2-buid humRibosomal.fa humRibosomal

２、chr19.faのindex作成 (bowit2を使う）

bowtie2-buid Homo_sapiens.GRCh38.dna_rm.chromosome.19.fa chr19

３、rRNAへのマッピング解析 (tophatを使う)

bowtie2 -p 8 --un DRR016695_1_without_rRNA.fastq -x humRibosomal -U DRR016695_1.fastq -S DRR016695_1_rRNA_mapped.sam

-p　スレッド数

--un <path>  write unpaired reads that didn't align to <path>

-x　rRNAインデックス作成時の名前

-U　インプットに使うunpairのfastq

-S　マッピング結果出力名

もう片方のfastqも同様に処理する。

bowtie2 -p 8 --un DRR016695_2_without_rRNA.fastq -x humRibosomal -U DRR016695_2.fastq -S DRR016695_2_rRNA_mapped.sam

rRNAにマッピングされなかったリードとして出力されたDRR016695_1_without_rRNA.fastqとDRR016695_2_without_rRNA.fastqを以降の解析に使う。

４、chr19へのマッピング解析

tophat -G Homo_sapiens.GRCh38.87.gtf -p 8 -o human_output --no-novel-juncs chr19 DRR016695_1_without_rRNA.fastq DRR016695_2_without_rRNA.fastq

-p　たぶんスレッド数。増やしすぎるとエラーになる。CPUが多くても10以内に止める。

-o　出力フォルダ名

-G　エキソンの識別に使うgtfファイル。上でダウンロードしたファイルを使う。

-no-novel-juncs　アノテーションにない新規のジャンクションを破棄

chr19　上のbowtie2のステップで作成したインデックス名。

-Gの代わりに-gを使うと、既存のエキソン情報を使いながら、新規スプライシング産物も探索するようになる。 

解析が終わると指定した出力フォルダの中にaccepted.bamができる。これを次の定量に使う。

5、fpkmの算出

cufflinks -G Homo_sapiens.GRCh38.87.gtf -N --compatible-hits-norm -M rRNA.gtf -o output_fpkm human_output/accepted_hits.bam

-o　出力フォルダ名

-G　エキソンの識別に使うgtfファイル。

-N --compatible-hits-norm　既存のアイソフォームにマッピングされるものだけでFPKMを計算

-M rRNA.gtfにマッピングされるリードを除く。

解析が終わると

 genes.fpkm_tracking

 isoforms.fpkm_tracking

transcripts.gtf

が指定フォルダに出力される。gene.trackは遺伝子単位の発現量、isoform.trackは愛想フォーム単位の発現量のファイル。中身を見ると

>head -3 pkm human/ genes.fpkm_tracking

となっている。

tracking id　遺伝子ID。

class code　アイソフォームのタイプ。基本的に-

nearest ref id　近傍の遺伝子ID。

gene id　遺伝子ID（１列目と同じ）。

gene short name　遺伝子名。

tss id　転写開始点ID。

locus　ポジション。

length　遺伝子の全長。

coverage　カバレッジ。

FPKM　サンプルFPKM。

FPKM conf lo　下限95%区間のサンプルFPKM。

 FPKM conf hi　上限95%区間のサンプルFPKM。

 FPKM status　サンプルのFPKMステータス。OKかLOW、Hi、FAIL。

 複数サンプルがある場合は

 q0_FPK（サンプル１）　q1_FPKM（サンプル２）などの名前で区別される。

一般的な 全体のフローは以下で説明しています。

featurecountsとedgeRを使うワークフロー。

ggplot2によるグラフ作成

Contig_Hunter_version_0.1.plで解析すると、-R_ggplot2というフォルダの中に４つのファイルが出力される。それをRに読み込ませる。*1

coverage.txt GC.txt

length.txt name.txt

Rのワーキングディレクトリに４つのファイルをコピーし、Rのターミナル環境で以下のコマンドを実行（ggplot2とscalesライブラリがインストールされてない場合は、予めインストールしておくこと）。

library(ggplot2)

library(scales)

( x <- read.table("coverage.txt", header=T) ) #カバレッジ

( y <- read.table("length.txt", header=T) ) #x軸も同様に入力

( z <- read.table("GC.txt", header=T) ) #3番目の要素

df <- data.frame(x=x, y=y,z)

g <- ggplot(df,aes(x = x,y = y,z))

g <- g + xlab("coverage")    # x 軸ラベル

g <- g + ylab("length")    # y 軸ラベル

#グラフ6 目盛りを統一

xbreaks <- c(10,100,500)

ybreaks <- seq(100,10000,1000000)

#最後に以下の１行をペースト

g + geom_point(aes(colour=z,size = GC),alpha = 0.5) + scale_size_area(name = "GC", max_size = 5) + scale_x_log10(labels=trans_format("log10",math_format(10^.x))) + scale_y_log10(labels=trans_format("log10",math_format(10^.x))) + scale_colour_gradientn(colours=rainbow(3)) + theme(axis.title.x = element_text(size=20)) + theme(axis.title.y = element_text(size=20)) + theme(axis.text.x = element_text(size=15)) + theme(axis.text.y = element_text(size=15)) + coord_cartesian(xlim=c(0.1,5000),ylim=c(100,2000000))

データに応じたグラフが得られるはず。このデータでは以下のようになった。最後に保存する。

ggsave(file="graph2.png")#セーブ

繰り返し配列やコンタミを含むより複雑なゲノムだと

これはSRAに登録されている海洋性のバクテリアのデータ。純化できていないとこのようなグラフになる。

最後はSRAに登録されているmetagenomeデータ。fastq １つで15GBある！。samを作ると100GBを超えてしまった。メタゲノムおそるべし。

様々な生物種のDNAが混じっている。ゲノムを決めるにはプロットをまとめてbinningしていく必要がある。

メモ

シミュレーションデータ

Escherichia_coli_O157-H7  short_read x100

Escherichia_coli_O157-H7  short_read x100 + pacbio x 30

Escherichia_coli_O157-H7  short_read x100 + pacbio x 100

*1:ここに脚注を書きます