macでインフォマティクス

macでインフォマティクス

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

シロイヌナズナの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 パスウェイの視覚化パッケージ