macでインフォマティクス

macでインフォマティクス

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

bedtools

BEDファイルのオーバーラップ領域を抽出したり、マージしたりできるツール。BED以外にGFF、VCFも扱うことができる。

 

オーバーラップする領域の抽出。

aのfeatureのうちbの遺伝子とオーバーラップする領域だけをコールする。

bedtools intersect -a reads.bed -b genes.bed -sorted > output.bed
  • -wa Write the original entry in A for each overlap.
  • -wo Write the original A and B entries plus the number of base pairs of overlap between the two features.
  • -c For each entry in A, report the number of overlaps with B.
  • -f Minimum overlap required as a fraction of A. - Default is 1E-9 (i.e., 1bp). - FLOAT (e.g. 0.50)

-waをつけるとbとオーバーラップするaの領域がそのまま出力され、-waをつけないとaとbでオーバーラップする領域が出力される(参考)。-woをつけると右端カラムにオーバーラップする領域のサイズ付きで出力される。-cをつけると右端カラムにオーバーラップする領域数付きで出力される。デフォルトでは1bp以上のオーバーラップを重複とみなすが、-fで指定することでその閾値をコントロールできる。厳しくするなら、例えば-f 0.5にすれば50%以上オーバーラップを重複に設定できる。

 

(複数ファイル同時)オーバーラップする領域の抽出。

aのfeatureのうち3つのファイルのオーバーラップする領域をコールする。

bedtools intersect -a reads.bed -b genes1.bed genes2.bed genes3.bed -sorted > output.bed

-wa -wb をつけると重複したfeature数がわかる。-namesをつけると、どのデータ由来かわかる(-names gene1 gene2 gene3 )。

 

オーバーラップしない領域の抽出。

aのfeatureのうちbの遺伝子領域にない領域だけをコールする。

bedtools intersect -a reads.bed -b genes.bed -v -sorted > output.bed

 

Aとオーバーラップし、Bとオーバーラップしない領域の抽出。

aのfeatureのうちLINEとオーバーラップしており、かつSINEとはオーバーラップしていないfeatureを抽出。

bedtools intersect -a genes.bed -b LINES.bed | \ 
bedtools intersect -a stdin -b SINEs.bed -v > output.bed

パイプと標準入力(stdin)を使い2回連続比較している。順次実行してもよい。

 

aのそれぞれの遺伝子について、最も近いALUをコール。

bedtools closest -a genes.bed -b ALUs.bed > output.bed

 

オーバーラップしているfeatureを統合して返す。

bedtools merge -i reads.bed > output.bed
  • -d Maximum distance between features allowed for features to be merged. - Def. 0. That is, overlapping & book-ended features are merged.

-dを1以上に指定することで、オーバーラップしてないが近縁にあるfeatureを統合できる(-d 1000なら1000bp以内のfeatureを統合できる)。 -c 4 -o count,collapseをつければ、どのfeatureがどのfeatureと何回オーバーラップしたのかも出力される。

 

マージコマンドは、chr順、スタートポジション順にソートされていないとエラーを起こす。unixのsortコマンドでソートしてから使う(sort -k1,1 -k2,2n リンク)。

 

 

 

 

 

 

 

 

IGVが重い時のtips

 

IGVにbamファイルを読み込むと重くて困ることがある。表示を軽くするためのtipsを書いておく。

 

bamを読み込んでリードが見えるまで拡大した状態。

f:id:kazumaxneo:20170722115958j:plain

 

view -> Preference

f:id:kazumaxneo:20170722120151j:plain

 

赤枠部分を以下のように修正。

f:id:kazumaxneo:20170722122106j:plain

 

リードがダウンサンプリングされ、動作が軽くなった。

f:id:kazumaxneo:20170722122822j:plain

 

まずvisiblity range threshold (kb) だが、この値はリードが表示されるようになるサイズを規定している。値を小さくすると、メモリーにキャッシュされるリード量が減り動作が高速化する (10kbなら10kb以下の縮尺にならないとリードを読み込まない)。またリードの表示数を調整するため、Downsample readsの項目を"max read count => 1"、"per window size (bases) => 100にして、100bpに1リードの割合までデータ減らして表示している。

ただしリードを減らしすぎるとレアな変異などのイベントが見えなくなってしまう。減らす場合、10リードくらいにとどめておく。

max read count を10まで増やした。

f:id:kazumaxneo:20170722123546j:plain

このようになった。

f:id:kazumaxneo:20170722123643j:plain

このくらいの数なら、複数のbamファイルを読み込んでも動作すると思われる。

 

そのほか下のオプション設定を追加しても良いかもしれない。

f:id:kazumaxneo:20170722131439j:plain

Filter duplicate readsにチェックをつけると、PCR duplicateと思われるリードが表示されなくなる。Filter secondary alignmentsにチェックをつけると、複数のアライメント部位を持つリードについては、アライメントスコアがベスト以外のアライメントが表示されなくなる。

show soft-clipped basesとshow center lineなどのチェックを外す。indel周辺などミスアライメントが補正できなかった部位はミスアライメントが多発している。そういった部位のリードはsoft-clipしてマッピングされている。ふつうミスアライメントが多発した部位の情報は必要ないのでチェックを外す。center lineも必要なければチェックを外す。

 

 

 

 

 

 

 

 

 

 

 

 

 

bamからbigWigとWiggle Formatに変換するツール

bamからwiggleファイルに変換してカバレッジのtrackをviewerのトラックに取り込みことができる。ただしそれにはsamtoolsのpileupを使いbamからwiggleファイルを作る必要があり、作り方がやや面倒だった。現在では、ありがたいことにコマンド一発でwiggleファイルを作るツールが公開されている。このツールを使ってbigwig/wiggleに変換し、IGVにカバレッジtrackを追加する流れを書いておく。

 

 

インストー

wiggleとbigwigに変換するツール

https://github.com/MikeAxtell/bam2wig

バイナリなのでpathを通して実行権を与えるだけでランできる。ただしbigwigに変換するには、下のリンクからwigtobigwigツールをダウンロードしてパスを通しておく必要がある。wigtobigwigがないと、wiggleファイルだけ作られる。

 

wigtobigwigのダウンロード - unix

Index of /admin/exe/linux.x86_64

 

wigtobigwigのダウンロード - mac OSX

 

 

ラン

bam2wig input.bam 

使いそうなオプション

-r  : Limit analysis to specified read group

-l : Minimum size of read to report. Default: no minimum.

-L : Maximum size of read to report. Default: no maximum.

-c  : Limit analysis to the indicated locus.

-m : Count each as (1/(n mapped reads)) * 1E6. In other words, normalize to reads per million.

ランが終わると、bamのあるディレクトリにbamファイル名のディレクトリが追加される。その中にbigwigとwig形式のファイルが作られる。

 

IGVで表示する。まずbamを取り込む。

igv -g input.fasta input.bam

f:id:kazumaxneo:20170720151839j:plain

IGVが起動してbamが表示された。さらにFIle -> Load From filesからbigwigファイルを追加する。

 

 

カバレッジが深い部位があると変化が見えにくい。その場合、パネルの左側で右クリックし、scaleをautoに変更する。

f:id:kazumaxneo:20170720151202j:plain

 

ついでにtrack のheightも変更する。ここでは250にした。

f:id:kazumaxneo:20170720151301j:plain

カバレッジ変化が見やすくなった。

f:id:kazumaxneo:20170720151927j:plain

 

 

bigwigはbainaryだが、wiggleファイルはテキストファイルである。下はwiggle変換ツールで作ったwiggleファイルの中身を見たものである。

variableStep chrom=AP017303.1 span=90

390 2

variableStep chrom=AP017303.1 span=2

522 2

variableStep chrom=AP017303.1 span=87

524 4

variableStep chrom=AP017303.1 span=1

611 5

variableStep chrom=AP017303.1 span=2

612 3

variableStep chrom=AP017303.1 span=4

614 1

variableStep chrom=AP017303.1 span=83

618 3

variableStep chrom=AP017303.1 span=7

IGVでこの部位のリードを見てみる。

f:id:kazumaxneo:20170720144929j:plain矢印部分がwiuggleファイルで表記されている部位である。

 

左側を拡大する。

f:id:kazumaxneo:20170720145200j:plain

wiglgleでは

variableStep chrom=AP017303.1 span=90

390 2

となっている。つまり390からspanの90ntだけカバレッジ1の領域が続くことを表している。

  

wiggle formatの詳細はUCSCのリンクを確認してください。

Genome Browser Wiggle Track Format

 

 

引用

https://github.com/MikeAxtell/bam2wig

 

 

 

 

eukaryotesのアノテーションツール; Augustus

 

Augustusはblast2goでも使われているeukaryotesのアノテーションツール。既存の他の手法と比較しても精度が高い手法と述べられている(検証リンク)。高速なwebサーバー版と、RNA-seqのbamファイルを指定してexon-intron情報を与え、予測精度を上げるlocal版がある。このbamによる支援機能はblast2goにも取り込まれている。

 local版はbinaryのみ提供されている。ただしmac OSXでは動かない。unixマシンやunixサーバーを用意する必要がある。

 

 

webサーバー版

Augustus: job submission

 

local版(breakerとaugustus両方必要)

http://bioinf.uni-greifswald.de/augustus/binaries/

 

ポスター

http://bioinf.uni-greifswald.de/bioinf/publications/pag2015.pdf

 

トレーニング済みゲノム一覧 (README.TXT)

http://bioinf.uni-greifswald.de/augustus/binaries/README.TXT

 近縁なゲノムであれば使えるらしい。例えば哺乳類であればヒトのトレーニングデータを使えると書かれている。

 

bamはtophatなどで作る。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

引用

 

アセンブル結果をCore gene setの検出数で評価する BUSCO

 

ゲノムのアセンブルやde novo transcriptomeの評価手法の1つに、Core gene setがアセンブルされた配列の中にあるかどうか調べるやり方がある。そのようなツールとしてCEGMAがよく知られており、version2.5までアップーデートされていた。が、2015年にはその後継のBUSCOが発表され、CEGMAはサポートが中止されている。BUSCOはCEGMA3.0とも言うべきツールで、CEGMAより随分高速化されているらしい。開発者のブログにその経緯が詳しく書かれている。リンクを貼っておきます。

Goodbye CEGMA, hello BUSCO! — ACGT

 

BUSCOは、90%以上の種でsingle copy orthologousな遺伝子をcore geneとして定義し、それをデータベースにしてクエリー配列を検索する。90パーセント以上と定義しているので、core geneヒット数が90パーセントを大きく割るようならおそらくアセンブルが不完全と判断できる。

 

マニュアルは公式HPからダウンロードできる。

BUSCO

 

インストー

brew install busco

*CEGMAもbrewでインストール可能。

 

ラン

解析には、Core geneのデータベースを用意する必要がある。データセットは以下のリンクからダウンロードする。バクテリアも用意されている。

 

シロイヌナズナのデータなら陸上植物のデータベースを使う。wgetでダウンロードする。

wget http://busco.ezlab.org/v2/datasets/embryophyta_odb9.tar.gz

ダウンロードが終わったら解凍する。

 

Genome assembly

BUSCO -g scaffolds.fasta -o OUTPUT -l embryophyta_odb9 -m genome -c 8
  • -g Input file in fasta format. Can be a genome, proteome or transcriptome.
  • -o output
  • -l lineage,lineage Which BUSCO lineage to be used.(ここではembryophyta_odb9を指定)
  • -m mode which module to run the analysis to run, valid modes are 'all'(genome assembly), 'OGS' (gene set / proteome) and 'trans' (transcriptome). (Defaults all)
  • -c Number of threads/cores to use.

 

Transcriptome

BUSCO -g transcripts.fasta -o OUTPUT -l embryophyta_odb9 -m trans -c 8

アラビのRNAをde novoでアセンブルした配列(fasta)をクエリーとしている。

  • -l ここではembryophyta_odb9を指定

アラビのtranscriptomeデータをde novo assemblyして作った配列(54000 transcrpts)を解析したところ、10分くらいかかった。

 

full_table_OUTPUT~が詳細なデータとなる。またshort_summary_OUTPUT_PLANTに結果がまとめられている。summaryには以下のような情報が載っている。

#Summarized BUSCO benchmarking for file: transcripts.fasta 

#BUSCO was run in mode: trans

Summarized benchmarks in BUSCO notation

 : C:0%[D:0%],F:0%,M:0%,n:1440

 1357 Complete BUSCOs

 896 Complete and single-copy BUSCOs

 461 Complete and duplicated BUSCOs

 17 Fragmented BUSCOs

 66 Missing BUSCOs

 1440 Total BUSCO groups searched

 

core geneのカタログは合計1440遺伝子あるが、そのうち66が見つからないと出ている。missing_buscos_list~には見つからなかった66のIDが載っている。orthoDBを検索して、どんな遺伝子がsingle copy orthologousとして見つからなかったか確認する。

Ortholog Search | cegg.unige.ch Computational Evolutionary Genomics Group

 

 今回の検証ではシロイヌナズナのリファレンスcDNAを使っているので、アセンブルの漏れは無いはずである。見つからなかった66遺伝子は2コピー以上存在するか、ゲノムに存在しないかのどちらかと考えられる。

 

引用

BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs.
Felipe A. Simão, Robert M. Waterhouse, Panagiotis Ioannidis, Evgenia V. Kriventseva, and Evgeny M. Zdobnov
Bioinformatics, published online June 9, 2015 | Abstract | Full Text PDF | doi: 10.1093/bioinformatics/btv351
Supplementary Online Materials: SOM

シロイヌナズナのRNA seq解析

 

勉強会用資料。

 

 

論文

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3566969/

 

シーケンスデータ

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP010481

 

シロイヌナズナfastaファイルとgtfファイル

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

 

 

1、index

bowtie2-build -f A_thaliana_TAIR10.fa chr

 

2、mapping

tophat2 -p 4 -I 50000 -i 20 -o tophat_WT1-1 -G Arabidopsis_thaliana.TAIR10.36.gtf chr WT_15flower_receptacles1-1/input.fastq
  • -p  num-threads  
  • -I  --max-intron-length (default: 500000)
  • -i  --min-intron-length  (default: 50)
  • -G --GTF
  • -o --output-dir (default: ./tophat_out)

 

3、rename

mv tophat_WT1-1/accepted_hits.bam tophat_WT1-1/WT1-1.bam

 

4、index

samtools index tophat_WT1-1/WT1-1.bam

step1-4を全てのサンプルに対して実行(またはスクリプトで処理)。

 

5、IGVでアライメントを確認

igv -g A_thaliana_TAIR10.fa 
Arabidopsis_thaliana.TAIR10.36.gtf,WT1-1/WT1-1.bam,WT1-2/WT1-2.bam,WT1-3/WT1-3.bam,MT1-1/MT1-1.bam,MT1-2/MT1-2.bam,MT1-3/MT1-3.bam

 

6、featurecountsでカウント。

featureCounts -T 8 -t exon -g gene_id -a Arabidopsis_thaliana.TAIR10.36.gtf 
-o counts.txt
WT1-1/WT1-1.bam WT1-2/WT1-2.bam WT1-3/WT1-3.bam
WT2-1/WT2-1.bam WT2-2/WT2-2.bam WT2-3/WT2-3.bam
WT3-1/WT3-1.bam WT3-2/WT3-2.bam WT3-3/WT3-3.bam
MT1-1/MT1-1.bam MT1-2/MT1-2.bam MT1-3/MT1-3.bam
MT2-1/MT2-1.bam MT2-2/MT2-2.bam MT2-3/MT2-3.bam
MT3-1/MT3-1.bam MT3-2/MT3-2.bam MT3-3/MT3-3.bam

 

7、データの分析(R)

library("ggplot2") 
library("reshape2")
 
rawCount = read.table("featureCounts_counts_trimmed.txt", header = TRUE, row.names = 1)
head(rawCount) #必ず毎回確認
artificialCount = log2(rawCount + 1) 
head(artificialCount) #必ず毎回確認
ggplot(artificialCount, aes(x = WT_15flower_receptacles1.1.accepted_hits.bam)) +
ylab(expression(log[2](count + 1))) +
geom_histogram(colour = "white", fill = "#525252", binwidth = 0.4)

f:id:kazumaxneo:20170715153703j:plain

 #サンプル1のraw count dataのヒストグラム

 

df <- melt(logdata) #data.frame形式に変換。 
head(df) #必ず毎回確認
g <- ggplot(df, aes (x = variable, y = value )) + geom_boxplot()
plot(g)

f:id:kazumaxneo:20170715154018j:plain

18サンプルのraw count dataのbox plot。  縦軸はlog2。

 

8、edgeRで正規化、DEG検出(R) 

次にedgeRでTMM正規化してbox plotを書く。

library("edgeR") 
count <- as.matrix(rawCount)
dim(count)
group <- factor(c("WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT","MT", "MT", "MT", "MT", "MT", "MT", "MT", "MT", "MT"))
d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)  #全遺伝子commonの分散を計算
d <- estimateCommonDisp(d) #moderated tagwise dispersionの計算)
d <- estimateTagwiseDisp(d) #exact test
result <- exactTest(d) 
topTags(result)
table <- as.data.frame(topTags(result, n = nrow(count))) #toptagsで全行指定してデータフレームに変換
write.table(table, file = "result.txt", col.names = T, row.names = T, sep = " ")
is.DEG <- as.logical(table$FDR < 0.05) #logical クラス(true、false)に判定付きで変換。
DEG.names <- rownames(table)[is.DEG] #trueだけDEG.namesに収納
plotSmear(result, de.tags = DEG.names, cex=0.3)

 FDR < 0.05

f:id:kazumaxneo:20170715154708j:plain

 

FDR < 0.01

f:id:kazumaxneo:20170715154833j:plain

 正規化係数とライブラリサイズの確認

正規化係数

> d$samples$norm.factors

 [1] 0.9888439 0.9993429 0.9877061 0.9953806 1.0011457 0.9931222 0.9891827

 [8] 0.9974448 0.9891228 0.9994781 1.0106604 1.0006687 1.0045275 1.0151970

[15] 0.9991256 1.0073002 1.0186367 1.0037970

ライブラリサイズ

> d$samples$lib.size

 [1]  8621707  6139876 10119340 10551274  7571024 12623646  8430112  5948333

 [9]  9897606  8022571  5732382  9460971  6524094  4626588  7563974  9157905

[17]  6514360 10694322

 TMM正規化してbox plotを描画。

f:id:kazumaxneo:20170715182417j:plain

TMM補正するとboxの高さがほぼ揃った。 縦軸はlog2。

 

参考

FPKM補正

f:id:kazumaxneo:20170716102942j:plain

 

f:id:kazumaxneo:20170716103527j:plain

 左からraw count、FPKM、TMM。

 

9、クラスタリング(R) 

クラスタリング

rawCount = read.table("featureCounts_counts_trimmed.txt", header = TRUE, row.names = 1)
count <- as.matrix(rawCount)
rho <- cor(count, method = "spearman")
d <- dist(1 - rho)
h <- hclust(d, method = "ward.D2")
plot(h)

f:id:kazumaxneo:20170715185844j:plain

 

WTとMTで仕分けられた。また容器ごとに仕分けられている。

 

10、ヒートマップ(R) 

FDR < 0.01のtrranscriptsを使いheat-mapを描く。

rawCount = read.table("result_FDR0.01.txt", header = TRUE, row.names = 1)
heatmap(as.matrix(data), margin=c(6,4), main="Heat Map 1 (Raw Data)")

f:id:kazumaxneo:20170715192734j:plain

 

 

data.z <- genescale(rawCount, axis=1, method="Z")
heatmap(as.matrix(data.z), margin=c(18,4), main="Heat Map 2")

 

f:id:kazumaxneo:20170715193312j:plain

 

heatmap.2(as.matrix(data.z), col=greenred(75), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=1, margin=c(6,8), main="Heat Map 3", Colv=NA)

f:id:kazumaxneo:20170715193602j:plain

 

 11、 主成分分析(R)

count <- as.matrix(count) 
count.log <- log10(count + 1)
count.log.var <- apply(count.log, 1, var)
count.log.pca <- count.log[order(count.log.var, decreasing = TRUE), ][1:200, ]
res <- prcomp(count.log.pca)
col <- c("red", "red", "red", "red", "red", "red", "red", "red", "red", "blue", "blue", "blue", "blue", "blue", "blue", "blue", "blue", "blue")
time <- c("1", "1", "1", "2", "2", "2", "3", "3", "3", "1", "1", "1", "2", "2", "2", "3", "3", "3")
pc1 <- res$rotation[, 1]
pc2 <- res$rotation[, 2]
pc3 <- res$rotation[, 3]
plot(pc1, pc2, type = "n", asp = TRUE)
text(pc1, pc2, labels = time, col = col)

f:id:kazumaxneo:20170715220318j:plain

青がWT、赤がMT。

plot(pc2, pc3, type = "n", asp = TRUE) 
text(pc2, pc3, labels = time, col = col)

f:id:kazumaxneo:20170715220344j:plain

 

12、DEGの保存。(R)

 log Fold changeが0以上のtranscripts(fold change >1:誘導)と、 log Fold changeが0以下のtranscripts(fold change <1: 抑制)を取り出す。edgeRのところで出力したresult.txtを使う。

input = read.table("result.txt", header = TRUE, row.names = 1)
head(input) #確認
significant <- input[(input$FDR <= 0.01),]
head(significant) #確認
write.table(induction, file = "significant.txt", col.names = T, row.names = T, sep = "¥t") #FDR<=0.01を保存しておく

#誘導遺伝子を抽出
induction <- significant[(significant$logFC >= 0),] #誘導
head(induction) #確認
write.table(induction, file = "induction.txt", col.names = T, row.names = T, sep = "¥t") #保存

抑制遺伝子を抽出
repression <- significant[(significant$logFC < 0),] #抑制
head(repression) #確認
write.table(repression, file = "repression.txt", col.names = T, row.names = T, sep = "¥t") #保存

 

 13、GO解析1(R)

library(org.At.tair.db) #リンク
library(goProfiles) #リンク

x <- read.csv("induction.txt", header=TRUE)#上の方でlogFC>=0条件で抽出したデータを読み込む。
head(x, n=3)
# name logFC logCPM PValue FDR
# 1 AT5G03820 3.811659 6.257744 0.00e+00 0.00e+00
# 2 AT2G02010 3.590547 5.904101 0.00e+00 0.00e+00
# 3 AT1G71380 2.365841 5.759736 0.00e+00 0.00e+00


induction <- basicProfile(genelist=x$name, onto='ANY', level=2, orgPackage="org.At.tair.db") #x$nameが遺伝子リスト

 basicProfileリンク)。

  • genelist List of genes on which the Profile has to be based
  • onto Ontology on which the profile has to be built( ANY, ’MF’,’BP’,’CC’).
  • Level Level of the ontology at which the profile has to be built
  • orgPackage Name of a Bioconductor's organism annotations package ('org.Xx-eg-db').

 

plotProfiles(induction) #plotProfilesリンク

 出力が以下になる。

f:id:kazumaxneo:20170718174229j:plain

ANYだと一部のカテゴリが削られている。 ’MF’,’BP’,’CC’に分けた方が良い。下はANYと同じパラメータで ’MF’,’BP’,’CC’を指定してグラフを書いている。

f:id:kazumaxneo:20170718175228j:plain

 

 抑制

y <- read.csv("repression.txt", header=TRUE)#上の方でlogFC<0条件で抽出したデータを読み込む。
head(y, n=3)
# name logFC logCPM PValue FDR
# 1 AT1G75945 -7.543735 2.836073 1.53e-272 8.35e-269
# 2 AT3G23460 -1.906906 4.256721 4.08e-145 8.37e-142
# 3 AT4G05715 -5.862613 2.371176 2.38e-97 2.37e-94

repression <- basicProfile(genelist=y$name, onto='MF', level=2, orgPackage="org.At.tair.db")
plotProfiles(repression)

f:id:kazumaxneo:20170718174618j:plain

 誘導、抑制を同時に書く。

result <- mergeProfilesLists(induction, repression, profNames = c("induction", "repression"))
plotProfiles(result, legendText = c("induction", "repression"), percentage = TRUE)

f:id:kazumaxneo:20170717124336j:plain

保存  

write.table(result, file = "GO_profile.txt")

テキストファイルが保存される。

MF.Description "MF.GOID" "MF.induction" "MF.repression"
13 "antioxidant activity" "GO:0016209" 33 12
9 "binding" "GO:0005488" 674 492
4 "catalytic activity" "GO:0003824" 852 609
1 "electron carrier activity" "GO:0009055" 13 4
14 "metallochaperone activity" "GO:0016530" 0 1
20 "molecular function regulator" "GO:0098772" 18 18
19 "molecular transducer activity" "GO:0060089" 7 18
3 "nucleic acid binding transcription factor activity" "GO:0001071" 62 96
6 "nutrient reservoir activity" "GO:0045735" 4 5
10 "protein tag" "GO:0031386" 1 0
5 "signal transducer activity" "GO:0004871" 14 25
7 "structural molecule activity" "GO:0005198" 238 10
2 "transcription factor activity, protein binding" "GO:0000988" 2 0
8 "transporter activity" "GO:0005215" 115 171

数が多い少ないでは意味が薄い。割合を出し検定を行い、エンリッチメントされた可能性がある系を絞り込む必要がある。

 

pathway

#先ほど読み込んだx(誘導された遺伝子を示したdata.frame)を使う。
library(pathview) #(pathviewインストール情報リンク

path.id <- "00010" #解糖系
path <- pathview(gene.data = as.matrix(x$name), pathway.id = path.id, species = "ath", gene.idtype = "KEGG")
  • gene.data either vector (single sample) or a matrix-like data (multiple sample).
  • pathway.id character vector, the KEGG pathway ID(s), usually 5 digit, may also include the 3 letter KEGG species code. 
  • species character, either the kegg code, scientific name or the common name of the tar- get species.
  • gene.idtype character, ID type used for the gene.data, case insensitive. Default gene.idtype="entrez", i.e. Entrez Gene, which are the primary KEGG gene ID for many common model organisms. For other species, gene.idtype should be set to "KEGG" as KEGG use other types of gene IDs.

 pathwayのpngがワーキングディレクトリに保存される。

f:id:kazumaxneo:20170718123936j:plain

解糖系のpathway。左はアラビでアノテーションされている遺伝子(緑)の図。右が誘導されていた遺伝子(赤)の図。

 

最初は全てのpathwayを見た方が良いかもしれない。しかし全部調べるのは面倒臭い。それならpathwayのIDを一括取得して、ループで回した方が便利なので、pathwayのIDを取得して全pathwayを描画する流れを書いておく。

 

1、Rを中断(crt + z)、または新しいターミナルのタブを開く。

2、まずkeggに登録された生物名をこのリンクから確認する。Arabidopsis属だと2つ登録されている。

T00041 ath Arabidopsis thaliana (thale cress) Eukaryotes;Plants;Eudicots;Mustard family 
T01578 aly Arabidopsis lyrata (lyrate rockcress) Eukaryotes;Plants;Eudicots;Mustard family

2、A.thaliana とA.lyrataが見つかるが、ここではA.thalianaを使う。左端にathと略号が書かれていることに注目。これをURL(http://rest.kegg.jp/list/pathway/)につけて、http://rest.kegg.jp/list/pathway/ath と検索する。

3、A.thalianaのpathway一覧が表示されている。

f:id:kazumaxneo:20170718144601j:plain

#参考  http://rest.kegg.jp/list/athだと遺伝子リスト。

 

4、このページ全体をコピペするか、wgetを叩きダウンロードする。wgetなら下のように打つ。

wget http://rest.kegg.jp/list/pathway/ath -O - |sed -e 's/path://g' > pathway_list

sedで先頭のpath:ath00010のpath:を消しながらダウンロードしている。

5、A.thalianaのkegg pathway IDのリストが入手できた。

user$ head -3  pathway_list

ath00010 Glycolysis / Gluconeogenesis - Arabidopsis thaliana (thale cress)

ath00020 Citrate cycle (TCA cycle) - Arabidopsis thaliana (thale cress)

ath00030 Pentose phosphate pathway - Arabidopsis thaliana (thale cress)

 

 

入手したIDリストを使い、全pathwayを描く。

cut -f 1 pathway_list > pathway_ID

Rに戻る(一時停止させているならfg %job番号で再開)

IDlist <- read.table("pathway_list2", header=FALSE) #pathway IDを読み込み

path <- pathview(gene.data = as.matrix(x$name), pathway.id = unlist(IDlist), species = "ath", gene.idtype = "KEGG")

IDがデータフレームのままだと認識してくれないので、unlist(IDlist)

でベクトルに変換している。使った全てのpathwayのグラフが描画される。アラビだと10分くらいかかる。

f:id:kazumaxneo:20170718144141j:plain

 全部は載せられないが、このように一括してpathwayを取得してマップを描画できる。もちろん、エンリッチされたpathwayがGO解析で分かっていれば、それだけ書いてもいい。

 

誘導、抑制をpathwayで表示できた方が便利そうなので、少し改良する。先ほど出力したFDR0.01以下のリストに遺伝子名とfold_changeが記載されている。これをRに読み込む。

significant <- read.table("significant.txt", header=TRUE)
head(significant, n=3) #確認
name logFC logCPM PValue FDR
1 AT5G03820 3.811659 6.257744 0.000000e+00 0.000000e+00
2 AT2G02010 3.590547 5.904101 0.000000e+00 0.000000e+00
3 AT1G71380 2.365841 5.759736 0.000000e+00 0.000000e+00

genelist <- as.vector(significant$logFC)
names(genelist) <- as.vector(significant$name)
path2 <- pathview(gene.data = genelist, pathway.id = unlist(IDlist), species = "ath", gene.idtype = "KEGG")

誘導、抑制が色でわかるようになった。

f:id:kazumaxneo:20170718153016j:plain

f:id:kazumaxneo:20170718152951j:plain

全部では無いが、その他の系。

 

 

library("clusterProfiler") #リンク
library(org.At.tair.db)
ids <- bitr(x$name, fromType="TAIR", toType="ENTREZID", OrgDb="org.At.tair.db") #ID変換にはclusterProfilerのbitrを使用(リンク
head(ids, n=3)
TAIR   ENTREZID
1 AT5G03820 831713
2 AT2G02010 814732
3 AT1G71380 843479

ggo <- groupGO(gene = ids$ENTREZID, OrgDb = org.At.tair.db, ont = "CC", level = 3, readable = TRUE)
barplot(ggo, drop=TRUE, showCategory=20, font=12)

f:id:kazumaxneo:20170722234219j:plain

 

all = read.table("result.txt", header = TRUE) #全遺伝子(バックグラウンドとして使う)
head(all, n=3)
name logFC logCPM PValue FDR
1 AT5G03820 -3.967569 6.623068 0 0
2 AT2G02010 -3.770837 6.244671 0 0
3 AT4G23690 -2.719579 6.855122 0 0

allgene <- bitr(all$name, fromType="TAIR", toType="ENTREZID", OrgDb="org.At.tair.db") #比較できるようにIDをTAIRからENTREZIDに変換。
ego <- enrichGO(gene = ids$ENTREZID, universe = allgene$ENTREZID, OrgDb = org.At.tair.db, ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE)
barplot(ego, showCategory=20, font=12)

f:id:kazumaxneo:20170723002756j:plain

dot plotで表現することも可能。

dotplot(ego,showCategory=20, font=12)

f:id:kazumaxneo:20170723002934j:plain

 

plotGOgraph(ego)

f:id:kazumaxneo:20170723004242j:plain

 

enrichMap(ego, n=10) 

f:id:kazumaxneo:20170723005837j:plain

 

 

 

 

 

 

Bioconductor のpathway解析ツール(117種類)

Bioconductor - BiocViews

 

 Bioconductor のenrichment解析ツール(82種類)

Bioconductor - BiocViews

 

引用 

TMM正規化

http://www.nathalievilla.org/doc/pdf/tutorial-rnaseq.pdf

 

edgeR

二群間比較(edgeR) | RNA-Seq による比較トランスクリプトーム解析

 

クラスタリング

階層クラスタリング | スピアマンの順位相関係数を利用したクラスタリング

 

成分解析

主成分分析 | RNA-Seq リードデータをクラスタリング

 

GO profile

goProfiles | 遺伝子セットの GO term の出現頻度を調べ棒グラフで視覚化

org.At.tair.db | シロイヌナズナのアノテーション Bioconductor

 

pathview

pathview | KEGG パスウェイの視覚化パッケージ

 

RNA seqのリードカウント featureCounts

RNA reqのリードカウントツールを紹介する。

featureCounts

 

ダウンロード

sourceforgeリンク

https://sourceforge.net/projects/subread/files/subread-1.5.2/

 

インストー

ソースオードをダウンロードして解凍し、/srcに移動。macでは以下のようにしてビルドする。

cd subread-1.4.6-source/src
make -f Makefile.MacOS

linux公式マニュアル参照。 

 

マニュアル

featureCounts - quick guide

WEHI Bioinformatics - featureCounts

  

ラン

inputはsam/bamファイルと、リファレンスのカウントするfeatureの場所を記したgtfファイルか、よりシンプルなSAFというフォーマットである。タブ区切りの5フィールドの形式で、公式サイトに例がある。

GeneID	Chr	Start	End	Strand
497097	chr1	3204563	3207049	-
497097	chr1	3411783	3411982	-
497097	chr1	3660633	3661579	-

 

シングルエンドのデータをカウント。

featureCounts -T 8 -t exon -g gene_id -a annotation.gtf -o counts.txt <input.bam>
  • -T Number of the threads. 1 by default.
  • -t Specify the feature type. Only rows which have the matched matched feature type in the provided GTF annotation file will be included for read counting. `exon' by default.
  • -g  Specify the attribute type used to group features (eg. exons) into meta-features (eg. genes), when GTF annotation is provided. `gene_id' by default.
  • -a Give the name of the annotation file. The program assumes that the provided annotation file is in GTF format. Use -F option to specify other annotation formats.
  • -o Give the name of the output file.

6サンプルカウントした結果のキャプチャが以下となる。

f:id:kazumaxneo:20170710110348j:plain

先頭部分数十行を載せている。出力はタブ区切りのテキストファイルで、1行目はコメント行(#)である。2行目以下は以下のような並びになっている。

‘Geneid’, ‘Chr’, ‘Start’, ‘End’, ‘Strand’,‘Length’,‘sample1_count’,‘sample2_count’ ...

2列目に1;1;1;1;1;1とあるのは、exsonが6あることを意味している。lengthは1つの遺伝子のオーバーラップをのぞいたexonの合計サイズ (bp) である。右端にあるのがカウントの列である。サンプル数分だけできる。

 

 

strand specificにカウント。

featureCounts -s 1 -T 8 -t exon -g gene_id -a annotation.gtf -o counts.txt <input.bam>
  • -s  Indicate if strand-specific read counting should be performed. It has three possible values: 0 (unstranded), 1 (stranded) and 2 (reversely stranded). 0 by default.

"-s 1"はfeatrureと同じ向きにアライメントされたリードだけカウント。"-s 2"はfeatureの反対向きにアライメントされたリードだけカウント。

 

ペアリードのフラグメントをカウント。

featureCounts -p -T 8 -t exon -g gene_id -a annotation.gtf -o counts.txt <input_PE.bam>
  • -p If specified, fragments (or templates) will be counted instead of reads. This option is only applicable for paired-end reads. The two reads from the same fragment must be adjacent to each other in the provided SAM/BAM file.

 

複数ファイルを同時にカウントして出力する。

featureCounts -T 8 -t exon -g gene_id -a annotation.gtf -o counts.txt library1.bam library2.bam library3.bam

bam/samをスペース区切りで記載する。

 

paired-endの両方がマッピングされたものだけカウント。

featureCounts -p -B -t exon -g gene_id -a annotation.gtf -o counts.txt <input_PE.bam>
  • -B If specified, only fragments that have both ends successfully aligned will be considered for summarization. This option is only applicable for paired-end reads.

 

紹介したオプション以外に、インサートサイズに下限、上限を設けるオプション(-P-d, -Dや、異なるクロモソームにアライメントされたpaired-endを排除するオプション(-C)、マッピングクオリティの閾値-Q)、複数箇所にマッピングされたリードのカウント(-M)、など多様なオプションがあります。詳細は公式quickマニュアルを確認してください。

 

 

 

 

引用

featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.

Liao Y, Smyth GK and Shi W (2014).  

Bioinformatics, 30(7):923-30.