Ma et al., 2021を参考に試行錯誤しながらRNA-seqの解析を行ったのでその方法を簡単に説明していきます！独学でなんとか終えた解析なので、もしここ変じゃない？？という箇所を見つけた方がいらっしゃいましたらぜひ改善しながら次の世代へと繋いでいただけると幸いです。
※この工程はノートパソコンだと耐えきれないかもしれません。また、これはRではなくコマンド（ターミナル）に打って操作するものです。
$ java -jar /Users/tomokahiguchi/Downloads/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 4 -trimlog log.txt SRR6050414_1.fastq SRR6050414_2.fastq SRR6050414_1_clean.fastq unpaired_SRR6050414_1.fastq SRR6050414_2_clean.fastq unpaired_SRR6050414_2.fastq ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:20 TRAILING:20 MINLEN:36-jar /Users/tomokahiguchi/Downloads/Trimmomatic-0.39/trimmomatic-0.39.jar
Trimmomatic-0.39.jarのパス。要はTrimmomaticをダウンロードし保存してあるフォルダのなかの実行ファイルであるtrimmomatic-0.39.jarがパソコン内のどこにあるのか示しています。もしパスを通してあれば実行ファイル名を記載するだけでよいです。Trimmomaticのバージョンが更新されている場合があるので、常に新しいのが出ていないか確認した方がいいと思います。
-trimlog log.txt
ログファイルの出力先名。カレントディレクトリ以外の場所に出力させたければファイル名のみならずフォルダ名も明記すること。
SRR6050414_1.fastq SRR6050414_2.fastq
インプットファイル名。ペアエンドなら同じサンプルに対して2つあるはず。
SRR6050414_1_clean.fastq unpaired_SRR6050414_1.fastq SRR6050414_2_clean.fastq unpaired_SRR6050414_2.fastq
ペアエンドリードの両方が残っているリード（paired output）とどちらか一方しか残らなかったリード（unpaired output）の出力先名。今後必要なのはpaired outputなのでunpairedのほうはすぐ消去してしまってよいです。
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10
TruSeq3-PE.faは除去対象となるアダプター配列をFASTA形式で記載してあるファイル名。主なものはusadellab/Trimmomaticにあります。私の場合は解析していただいた会社から教えてもらったアダプター配列をもとに独自のfastaファイルを作成しました。
LEADING:20 TRAILING:20 MINLEN:36
これはオプションなので参考に挙げたurlを見て各人の好きなようにどうぞ。
※もしファイルが沢山あっていちいち回すのがめんどくさい！という私みたいな人は「bashシェルスクリプト」を調べてみてください。以下のような感じでbashを用いてfor文で回したら早いです。
SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)
for seqlib in ${SEQLIBS[@]}; do
java -jar trimmomatic-0.38.jar PE -phred33 -threads 4 \
${seqlib}_1.fastq.gz ${seqlib}_2.fastq.gz \
${seqlib}_1.clean.fastq.gz ${seqlib}_1.unpaired.fastq.gz \
${seqlib}_2.clean.fastq.gz ${seqlib}_2.unpaired.fastq.gz \
ILLUMINACLIP:adapters.fa:2:30:10 LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:60
done※この工程はノートパソコンだと耐えきれないかもしれません。また、これはRではなくコマンド（ターミナル）に打って操作するものです。
まず、レファレンスにしたいトランスクリプトーム配列をダウンロードします。たいていはEnsembl Plantsからダウンロードするのかなと思いますが、N. benthamianaの場合はそこには載っていないのSolgenomicsのFTP SiteからcDNA配列をダウンロードしました。これを元にインデックスを作成します。
インデックスができたら疑似的マッピングと定量を行います。Trimmomaticで得たpaired outputをインデックスに貼り付けていくイメージです。
#インデックス作成
$ salmon index -t Niben101_annotation.transcripts.fasta -i Niben101_index
#疑似的マッピングと定量
$ salmon quant -i Niben101_index -l A -1 SRR6050414_1_clean.fastq -2 SRR6050414_2_clean.fastq -o salmon_quant_result --gcBias–gcBias
このフラッグは場合によって付けなくてもいいみたいですが、常に付けても悪いことはないと公式サイトに書いてあったので付けています。（参考：Salmon 1.10.0 documentation）
※ここからRでの作業になります。Ma et al., 2021をベースにしているのでそちらも参考に。
以降使う関数の定義とライブラリの読み込みを行います。
checkZeros <- function(v, threshold) {
res <- ifelse(sum(v == 0) > threshold, FALSE, TRUE)
return(res)
}
check <- function(v, threshold) {
require('magrittr')
#split内の数値はサンプル数によって異なるので適宜変更してください。私の場合はn=4の3処理区があったのでこのように指定しています。1,1,1,1,2,2,2,2,3,3,3,3と同じ処理区のサンプルたちには同じ数値がふられ、1,2,3それぞれに対してcheckZeros関数を適用しているイメージです。
res <- v %>%
split(c(rep(1, 4), rep(2, 4), rep(3, 4))) %>%
sapply(checkZeros, threshold) %>%
all
return(res)
}
library('DESeq2')
library('tidyverse')
library("tximport")
library("dplyr")
library('ggplot2')
library('limma')
library('pvclust')
library('readr')
library('RColorBrewer')
library('gridExtra')
library('ComplexHeatmap')tximportにSalmonの定量結果を読み込ませます。
samples <- read.table("Treatment_groups.txt", header = TRUE)
samples$RNA_extraction <- factor(samples$RNA_extraction) #have to change the class of this column from int to factor. See https://support.bioconductor.org/p/78516/
salmon.files <- file.path(list.files(".", pattern="salmon"), "quant.sf")
names(salmon.files) <- samples$sample
tx2gene <- read.table('Niben101_tx2gene.txt', header = TRUE)
txi <- tximport(salmon.files, type="salmon", tx2gene=tx2gene)RLE正規化を行い、DEGsを検出します。これについては色々な人が分かりやすく記述してくれていると思うので、都度検索してみてください。 もしサンプル間に、RNA抽出した日の違いなどのバッチエフェクトが存在するならばdesign式に組込むようにしてください。 バッチエフェクトがあるかどうか分からないという場合はsvaというライブラリーで自動検出する方法もあります。
まずは正規化します。
#designにRNA_extractionを入れておくことでバッチエフェクトとして認識してもらう
ddsTxi <- DESeqDataSetFromTximport(txi,
colData = samples,
design = ~ RNA_extraction + treatment)
nrow(ddsTxi) #56883
ddsTxi$treatment <- relevel(ddsTxi$treatment, ref = "GUS")
ddsTxi %<>%
estimateSizeFactors %>%
counts(normalized = TRUE) %>%
apply(1, check, 1) %>%
ddsTxi[., ]
ddsTxi <- DESeq(ddsTxi)
# estimating size factors
# using 'avgTxLength' from assays(dds), correcting for library size
# estimating dispersions
# gene-wise dispersion estimates
# mean-dispersion relationship
# final dispersion estimates
# fitting model and testing
rld <- rlog(ddsTxi)
ntd <- normTransform(ddsTxi)
cond <- list(c('WT', 'GUS'),
c('WT', 'GAA'),
c('GAA', 'GUS'))
resRaw <- lapply(cond,
function(x) {
ddsTxi %>%
results(contrast = c('treatment', x), alpha=0.05) %T>%
summary %>%
as_tibble %>%
dplyr::select(pvalue, padj, log2FoldChange, lfcSE) %>%
rename_all(.funs = list(~paste0(paste(x, collapse = '_vs_'), '_', .)))
}) %>%
bind_cols
# out of 37021 with nonzero total read count
# adjusted p-value < 0.05
# LFC > 0 (up) : 5472, 15%
# LFC < 0 (down) : 5209, 14%
# outliers [1] : 9, 0.024%
# low counts [2] : 1436, 3.9%
# (mean count < 5)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results
#
#
# out of 37021 with nonzero total read count
# adjusted p-value < 0.05
# LFC > 0 (up) : 2240, 6.1%
# LFC < 0 (down) : 1733, 4.7%
# outliers [1] : 9, 0.024%
# low counts [2] : 2871, 7.8%
# (mean count < 7)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results
#
#
# out of 37021 with nonzero total read count
# adjusted p-value < 0.05
# LFC > 0 (up) : 2562, 6.9%
# LFC < 0 (down) : 2195, 5.9%
# outliers [1] : 9, 0.024%
# low counts [2] : 2154, 5.8%
# (mean count < 6)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results
res <- cbind.data.frame(as.matrix(mcols(ddsTxi)[, 1:10]), counts(ddsTxi, normalize = TRUE), stringsAsFactors = FALSE) %>%
rownames_to_column(., var = 'ID') %>%
as_tibble %>%
bind_cols(resRaw) %>%
dplyr::select(ID, E1.LDB9130 : GAA_vs_GUS_lfcSE)
write_csv(res, 'HEL_WT_vs_GAA_Niben101.csv')
#こんな感じで時々結果を出力しておくとのちのち楽になることがあります。何度も繰り返して試行錯誤するのでいちいち頭からRを動かすのは非常に骨が折れます…。limmaを用いてバッチエフェクト除去を行います。
dat <- rld %>%
assay %>%
{.[rowMeans(.) > 1, ]}
group <- samples$treatment
design <- model.matrix(~ group)
rldData <- dat %>%
removeBatchEffect(batch = samples$RNA_extraction,
design = design)
write_csv(as.data.frame(dat), 'HEL_WT_vs_GAA_wo_lowcount_Niben101.csv')
save(ddsTxi, rldData, file = 'ddsTxi_HEL_vs_GAA_Niben101.RData')最後に発現変動遺伝子(DEGs)のみ取り出し、heatsig変数に入れます。
scaleC <- rldData %>%
t %>%
scale %>%
t %>%
as.data.frame %>%
rownames_to_column('ID') %>%
as_tibble %>%
rename_at(-1, .funs = list(~paste0('Scale_', .)))
rawC <- rldData %>%
as.data.frame %>%
rownames_to_column('ID') %>%
as_tibble %>%
rename_at(-1, .funs = list(~paste0('Raw_', .)))
degresC <- res %>%
select(ID, WT_vs_GUS_pvalue : GAA_vs_GUS_lfcSE)
heatPlot <- rawC %>%
inner_join(scaleC) %>%
inner_join(degresC)
fcsig <- heatPlot %>%
dplyr::select(ends_with('FoldChange')) %>%
transmute_all(list(~ case_when(. > log2(1) ~ 1,
. < -log2(1) ~ -1,
TRUE ~ 0)))
padjsig <- heatPlot %>%
dplyr::select(ends_with('padj')) %>%
`<`(0.05) %>%
as_tibble %>%
transmute_all(list(~ if_else(is.na(.), FALSE, .)))
heatsig <- (padjsig * fcsig) %>%
as_tibble %>%
abs %>%
rowSums %>%
{. >= 1} %>%
which %>%
dplyr::slice(res, .)
heatsig %>%
mutate_at(c('ID'), .funs = list(~if_else(is.na(.), '', .))) %>%
write_csv('Niben101_sig.csv')主成分分析（PCA）では実験計画がうまくいっているか確認できます。処理区ごとにクラスターを作ってくれれば、それはその処理の影響が一番大きく表れているということになります。
cols <- colData(rld)[, 1] %>% factor(levels = c('#a6cee3', '#1f78b4', '#e31a1c'))
pca <- prcomp(t(rldData))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
percentVar <- round(100 * percentVar)
pca1 <- pca$x[,1]
pca2 <- pca$x[,2]
pcaData <- data.frame(PC1 = pca1, PC2 = pca2, Group = colData(rld)[, 1], ID = rownames(colData(rld))) %>%
mutate(Treatment = rep(c('GUS', 'HEL-GAA', 'HEL-WT'), c(4, 4, 4)) %>% factor)
ggplot(pcaData, aes(x = PC1, y = PC2, colour = Treatment)) +
geom_point(size = 4) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
scale_colour_manual(values = c('#000000', '#4DBBD5B2', '#E64B35B2')) +
stat_ellipse(aes(x = PC1, y = PC2, group = Treatment), type = 't', level = 0.95) + #95%信頼楕円
coord_fixed(1) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 12, face = 'bold'),
legend.text.align = 0,
axis.text = element_text(size = 13),
axis.title = element_text(size = 14),
legend.text=element_text(size= 13),
legend.title = element_text(size = 14))
ggsave('PCA_wo_outliers_Niben101.pdf', width = 5, height = 5)
ggsave('PCA_wo_outliers_Niben101.jpg', width = 5, height = 5)各処理と相関を持って発現する遺伝子群を見つけます。方法は参考に挙げたものとほぼ同じです。
library(genefilter)
library(WGCNA)
library(tidyr)
library(igraph)
library(limmaDE2)
library(DescTools)
library(goseq)
datTraits <- read.csv("WGCNA_datTrait.csv", header = TRUE, row.names=1)
vsd <- varianceStabilizingTransformation(ddsTxi)
wpn_vsd <- getVarianceStabilizedData(ddsTxi)
dim(wpn_vsd) #[1] 37021 12
rv_wpn <- rowVars(wpn_vsd)
summary(rv_wpn)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.00183 0.04950 0.08282 0.16054 0.15446 6.82883
q75_wpn <- quantile( rowVars(wpn_vsd), .75) # <= original
q95_wpn <- quantile( rowVars(wpn_vsd), .95) # <= changed to 95 quantile to reduce dataset
expr_normalized <- wpn_vsd[ rv_wpn > q95_wpn, ]
allowWGCNAThreads()
input_mat = t(expr_normalized)
powers = c(c(1:10), seq(from = 12, to = 50, by = 2))
sft = pickSoftThreshold(
input_mat, # <= Input data
#blockSize = 30,
powerVector = powers,
networkType = "signed",
verbose = 5
)
# pickSoftThreshold: will use block size 1851.
# pickSoftThreshold: calculating connectivity for given powers...
# ..working on genes 1 through 1851 of 1851
# Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
# 1 1 0.628000 3.48000 0.5560 1070.0 1130.0 1290
# 2 2 0.312000 0.86700 0.3280 750.0 804.0 1030
# 3 3 0.169000 0.43600 0.1230 580.0 603.0 883
# 4 4 0.098900 0.26000 -0.0761 475.0 474.0 789
# 5 5 0.038600 0.13500 -0.2080 404.0 384.0 722
# 6 6 0.000114 -0.00625 -0.2620 352.0 321.0 672
# 7 7 0.041000 -0.10900 -0.2070 312.0 274.0 631
# 8 8 0.143000 -0.19400 -0.0697 279.0 240.0 597
# 9 9 0.260000 -0.27200 0.0822 253.0 211.0 567
# 10 10 0.409000 -0.34800 0.2750 231.0 188.0 541
# 11 12 0.591000 -0.46000 0.5090 195.0 152.0 497
# 12 14 0.731000 -0.53800 0.6720 168.0 124.0 459
# 13 16 0.793000 -0.60900 0.7450 147.0 105.0 427
# 14 18 0.818000 -0.66100 0.7730 129.0 88.7 399
# 15 20 0.836000 -0.70200 0.7940 115.0 76.1 374
# 16 22 0.893000 -0.74500 0.8640 103.0 64.9 351
# 17 24 0.915000 -0.77300 0.8910 92.7 56.1 331
# 18 26 0.931000 -0.80500 0.9120 83.8 48.7 313
# 19 28 0.932000 -0.83700 0.9130 76.2 41.9 297
# 20 30 0.938000 -0.86600 0.9220 69.5 37.1 282
# 21 32 0.943000 -0.89600 0.9310 63.6 32.5 268
# 22 34 0.933000 -0.92200 0.9220 58.4 28.9 256
# 23 36 0.924000 -0.94700 0.9140 53.8 26.0 244
# 24 38 0.918000 -0.97000 0.9100 49.7 23.1 233
# 25 40 0.913000 -0.99700 0.9070 46.0 20.6 223
# 26 42 0.909000 -1.02000 0.9060 42.7 18.4 213
# 27 44 0.916000 -1.04000 0.9180 39.7 16.6 205
# 28 46 0.920000 -1.05000 0.9230 37.0 15.0 196
# 29 48 0.922000 -1.06000 0.9290 34.5 13.4 188
# 30 50 0.926000 -1.07000 0.9340 32.3 12.4 181
png("SoftThreshold_powers.png", width = 200, height = 100)
par(mfrow = c(1,2));
cex1 = 0.9;
#ここからは一応確認するためのグラフ作図
dev.new()
plot(sft$fitIndices[, 1],
-sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
xlab = "Soft Threshold (power)",
ylab = "Scale Free Topology Model Fit, signed R^2",
type="n",
main = paste("Scale independence")
)
text(sft$fitIndices[, 1],
-sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
labels = powers, cex = cex1, col = "red"
)
abline(h = 0.90, col = "red")
plot(sft$fitIndices[, 1],
sft$fitIndices[, 5],
xlab = "Soft Threshold (power)",
ylab = "Mean Connectivity",
type = "n",
main = paste("Mean connectivity")
)
text(sft$fitIndices[, 1],
sft$fitIndices[, 5],
labels = powers,
cex = cex1, col = "red")
dev.off()
#scale free性の確認
Adj <- adjacency(input_mat, power = sft$powerEstimate)
row.names(Adj) = colnames(input_mat)
colnames(Adj) = colnames(input_mat)
thre <- 0.75
subAdj <- (Adj>thre)*Adj
k <- Adj %>%
apply(2, sum, na.rm=T) %>%
as.vector()
pdf("CheckScaleFree.pdf")
par(mfrow=c(1,2))
hist(k)
scaleFreePlot(k, main="Check scale free topology\n")
dev.off()上記で得られたsoft-thresholding powerを用いて実際にネットワークを構築していきます。
picked_power = sft$powerEstimate
temp_cor <- cor
cor <- WGCNA::cor # Force it to use WGCNA cor function (fix a namespace conflict issue)
set.seed(123)
netwk <- blockwiseModules(input_mat, # <= input here
# == Adjacency Function ==
power = picked_power, # <= power here
networkType = "signed",
# == Tree and Block Options ==
deepSplit = 2,
pamRespectsDendro = F,
# detectCutHeight = 0.75,
minModuleSize = 30,
maxBlockSize = 4000,
# == Module Adjustments ==
reassignThreshold = 0,
mergeCutHeight = 0.25,
# == TOM == Archive the run results in TOM file (saves time)
saveTOMs = T,
saveTOMFileBase = "ER",
# == Output Options
numericLabels = T,
verbose = 3)
cor <- temp_cor
pdf("WGCNA_Cluster_Dendrogram.pdf")
# Convert labels to colors for plotting
mergedColors = labels2colors(netwk$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(netwk$dendrograms[[1]], mergedColors[netwk$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()
module_df <- data.frame(
gene_id = names(netwk$colors),
colors = labels2colors(netwk$colors)
)
write.table(module_df, file = "WGCNA_gene2module.csv", row.names = F, sep = ",", quote = FALSE)
#Display the correlation values within a heatmap plot~~~~~~~~~~~~~~~~~~~~~~~~~~~
nGenes = ncol(input_mat)
nSamples = nrow(input_mat)
MEs0 <- moduleEigengenes(input_mat, module_df$colors)$eigengenes
MEs0 <- orderMEs(MEs0)
module_order = names(MEs0) %>% gsub("ME","", .)
moduleTraitCor = cor(MEs0, datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
pdf("WGCNA_Module-trat_Relationships.pdf")
par(mar = c(6, 8.5, 3, 3))
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = c("GUS", "HEL-GAA", "HEL-WT"),
yLabels = names(MEs0),
ySymbols = names(MEs0),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.8,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
dev.off()構築したネットワークを実際に作図していきます。 まずはlimmaDE2のwgcna2igraph関数を用いて、WGCNAで構築したネットワークをigraphオブジェクトに変換します。(参考:Demonstration of limmaDE2)
graph<-wgcna2igraph(net = netwk, datExpr = input_mat,
modules2plot = c(1), #turquoise
colors2plot = c("gray"),
kME.threshold = 0.9, adjacency.threshold = 0.9,
adj.power = sft$powerEstimate, verbose = T,
node.size = 1, frame.color = "black", node.color = "turquoise",
edge.alpha = .5, edge.width =1)
# using KME values to find genes with high module membership
# subsetting to modules: 1
# culling edges by adjacency
# removing unconnected nodes
# coverting to igraph format
# returning a network with 234 nodes and 1428 edgesigraphを用いて作図します。
library(tibble)
hs <- hub_score(graph)$vector
fine = 500
palette = colorRampPalette(c('magenta','green'))
graphCol = palette(fine)[as.numeric(cut(hs, breaks = fine))]
V(graph)$color <- graphCol
set.seed(123) #seedを設定しないと再現性がなくなります
#遺伝子名なしバージョン
pdf("WGCNAnet_turquoise.pdf")
plot(graph, vertex.size=hs*10, main="Module turquoise")
ColorLegend(x = "topright", width=0.4, height=0.05, col=palette(fine), horiz=TRUE, labels = c("Min", "Max"), title = "Hub Score")
dev.off()
#遺伝子名ありバージョン
pdf("WGCNAnet_turquoise_w_name.pdf")
plot(graph, vertex.size=hs*10, vertex.label = V(graph)$name, vertex.label.cex=0.05, main="Module turquoise")
ColorLegend(x = "topright", width=0.4, height=0.05, col=palette(fine), horiz=TRUE, labels = c("Min", "Max"), title = "Hub Score")
dev.off()
#一応変数を保存しておきます
MEs = netwk$MEs
geneTree = netwk$dendrograms[[1]]
save(netwk, graph, MEs, geneTree,
file = "WGCNA_NetworkConstruction.RData")
HubScore_df <- inner_join(rownames_to_column(data.frame(hs)), module_df, by = c("rowname" = "gene_id")) %>%
column_to_rownames("rowname")
write.csv(HubScore_df, "turquoise_HubScore.csv")ここからは着目したいモジュールのGOおよびKEGGパスウェイ解析を行います。 まずは、padj < 0.05 & |log2FC| > log2(1)の遺伝子について、発現上昇遺伝子を1、減少遺伝子を-1、変動していない遺伝子を0と置き換えます。(sig変数に代入)
fcsig <- heatsig %>%
select(ends_with('FoldChange')) %>%
transmute_all(list(~ case_when(. > 1 ~ 1,
. < -1 ~ -1,
TRUE ~ 0)))
padjsig <- heatsig %>%
select(ends_with('padj')) %>%
abs %>%
`<`(0.05) %>%
as_tibble %>%
transmute_all(list(~ if_else(is.na(.), FALSE, .)))
sig <- (padjsig * fcsig) %>%
as_tibble %>%
mutate(ID = heatsig$ID) %>%
select(ID, everything()) %>%
dplyr::rename(WT_vs_GUS = WT_vs_GUS_padj,
WT_vs_GAA = WT_vs_GAA_padj,
GAA_vs_GUS = GAA_vs_GUS_padj) %>%
tibble::column_to_rownames('ID')
Module_sum <- group_by(module_df, colors) %>%
summarise(., uu = n_distinct(gene_id))
sig <- inner_join(rownames_to_column(sig), module_df, by = c("rowname" = "gene_id")) まずはGO解析についてです。これは遺伝子長も加味した解析を行うことのできるgoseqを用いて行います。各遺伝子の長さの情報が必要です。
A. thalianaなどの主要な植物についてはデフォルトでgoseqに遺伝子長が登録されているようですが、N. benthamianaについては登録がないので自分で与えてやる必要があります。
そのため、GenomicFeaturesを用いて遺伝子の長さをSolgenomicsに公開されているgffから取り出してNiben101_genelength.csvに追加しました。(参考:Non-Overlapping Exon Length の計算)
txdb <- makeTxDbFromGFF('Niben101_annotation.gene_models.gff', format = 'gff')
exons_list_per_gene <- exonsBy(txdb, by = 'gene')
exonic_gene_sizes <- lapply(exons_list_per_gene, function(x){sum(width(reduce(x)))})
df_exonic_gene_sizes <- as.data.frame(exonic_gene_sizes)
write.csv(as.data.frame(t(df_exonic_gene_sizes)), "Niben101_genelength.csv")また、遺伝子とGOタームの組み合わせについてもN. benthamianaについては登録がないので、自分で与えてやる必要があります。(参考:Retrieving GO mappings for GOseq analysis with non-native organism)
ここではSolgenomicsからダウンロードしたNiben101_annotation.functional.txtというファイルを元に1つのGOタームごとに対応する遺伝子を列挙する形でgene2cat変数に入れました。
annotations <- read.delim("Niben101_annotation.functional.txt", skip=2)
annotations %>%
select(Protein.Accession, Gene.Ontology.ID..Name.) -> annotations
annotations[annotations$Gene.Ontology.ID..Name. == "", ]$Gene.Ontology.ID..Name. <- NA
annotations <- na.omit(annotations)
annotations %>%
separate_rows(Gene.Ontology.ID..Name., sep="\\), ") %>%
mutate(Gene.Ontology.ID..Name. = substring(Gene.Ontology.ID..Name., 1, 10)) %>%
mutate(Protein.Accession = gsub(Protein.Accession, pattern = "\\.1$", replacement = "")) -> annotations
gene2cat <- split(annotations$Protein.Accession, annotations$Gene.Ontology.ID..Name.)goseqによりGO解析を行います。
vsGroup <- c("WT_vs_GUS", "WT_vs_GAA", "GAA_vs_GUS")
gene_length <- read.csv('Niben101_genelength.csv')
colnames(gene_length) <- c('ID', 'length')
for(j in vsGroup){
for (i in c(4, 7)){
degVec <- sig %>%
transmute((!!as.name(j)) != 0 &
colors == Module_sum$colors[[i]]) %>%
unlist %>%
as.integer
names(degVec) <- sig$rowname
bias.data <- inner_join(sig, gene_length, by = c("rowname" = "ID"))
bias.data <- bias.data$length
pwf <- nullp(degVec, bias.data = bias.data)
GOTestWithCat <- goseq(pwf, gene2cat = gene2cat, use_genes_without_cat = TRUE) %>%
as_tibble %>%
filter(!is.na(ontology))
write.csv(GOTestWithCat, paste0('Mo_', Module_sum$colors[[i]],'_', j, '_GO_DEGs.csv'))
}
}GO解析結果の作図とKEGGパスウェイ解析を行います。 パスウェイ解析はclusterProfilerを用いますが、ここに登録されている生物以外はKOを用いて解析する必要があります。 KEGGMatには、KofamKOALAにアミノ酸配列を投げた結果を載せました。
library(foreach)
KEGGMat <- read.csv("Niben101_gene2KO.csv", header = T)
keggsig <- dplyr::inner_join(sig, KEGGMat, by = c("rowname" = "ID"))
for (i in c(4, 7)){
vsGroup <- c("WT_vs_GUS", "WT_vs_GAA", "GAA_vs_GUS")
pathPlot <- foreach(k = seq_along(vsGroup), .combine = bind_rows) %do% {
vsGroup[k] %>%
{paste0('Mo_', Module_sum$colors[[i]], '_', ., '_GO_DEGs.csv')} %>%
read.csv %>%
select(term, over_represented_pvalue, numDEInCat, numInCat, ontology) %>%
rename(pvalue = over_represented_pvalue) %>%
mutate(group = vsGroup[k], GeneRatio = numDEInCat / numInCat) %>%
filter(pvalue < 0.05 &
numDEInCat >= 1)
}
pathPlot %>%
mutate(group =ordered(group, levels = unique(group))) %>%
ggplot(aes(x = group, y = term)) +
geom_point(aes(size = GeneRatio, colour = pvalue)) +
scale_color_gradient(low = "red", high = "blue") +
#scale_colour_gradientn(name = '-log10(P-value)', limits=c(0, max(-log10(pathPlot$pvalue))), colours = colorPal) +
scale_x_discrete(labels = c('HEL-WT vs GUS', 'HEL-WT vs HEL-GAA', 'HEL-GAA vs GUS')) +
ylab('') +
xlab('') +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_grid(ontology ~ ., scales = "free")
ggsave(paste0('Mo_', Module_sum$colors[[i]], '_GO_dotplot.pdf'), width = 6, height = 5)
ggsave(paste0('Mo_', Module_sum$colors[[i]], '_GO_dotplot.jpg'), width = 6, height = 5)
keggID <- list(
keggsig[keggsig$WT_vs_GUS != 0 & keggsig$colors== Module_sum$colors[[i]],]$KO,
keggsig[keggsig$WT_vs_GAA != 0 & keggsig$colors== Module_sum$colors[[i]],]$KO,
keggsig[keggsig$GAA_vs_GUS != 0 & keggsig$colors== Module_sum$colors[[i]],]$KO
)
names(keggID) <- c("HEL-WT vs GUS", "HEL-WT vs HEL-GAA", "HEL-GAA vs GUS")
#clusterProfiler
x <- compareCluster(keggID, fun = "enrichKEGG", organism='ko', keyType='kegg', pvalueCutoff = 0.05, minGSSize = 10, maxGSSize = 500, pAdjustMethod = "BH", qvalueCutoff = 0.05)
dotplot(x) + theme(axis.text.x=element_text(angle=45, hjust=1)) #https://github.com/YuLab-SMU/clusterProfiler/issues/80
ggsave(paste0("Mo_", Module_sum$colors[[i]],"KEGG_dotplot.pdf"), width = 6, height = 8)
}