【R】遺伝子発現データのPCAとかヒートマップとかしたい

・お題：遺伝子発現データセットを入手した。遺伝子発現データのパターンを俯瞰するために、PCAで次元削減してサンプル間の関係を可視化したり、ヒートマップを作成して遺伝子発現のパターンを可視化したい。

・今回題材にしたのは、GSE269178という遺伝子発現データセット。

www.ncbi.nlm.nih.gov

・遺伝子組換えマウス/野生型マウス × 試験物質/媒体、の4条件で、各群n = 3でRNA-Seqしたものらしい。

・前準備として、先のリンクのカウントデータ（GSE269178_BMDM.gene_counts.tsv.gz）をダウンロードして、カレントディレクトリに保存しておく。

# ライブラリ読み込み

library(tidyverse) # いろいろ便利になる
library(magrittr) # パイプなど
library(GEOquery) # GEOからデータ釣ってくる
library(edgeR) # 遺伝子発現データのオブジェクトを作ったり
library(RColorBrewer) # きれいな色付けに使用
library(matrixStats) # row varianceなど

library(pheatmap) # ヒートマップ

library(psych) # ヒストグラム

# カウントデータ読み込み

df_counts <-
read.delim("GSE269178_BMDM.gene_counts.tsv.gz")

df_counts %>% head()

# メタデータ読み込み

gsem <- getGEO(GEO = "GSE269178", GSEMatrix = TRUE, AnnotGPL = FALSE)
phenoData <- pData(gsem[ [ 1 ] ] ) # [ ]にスペースは不要だが、ブログの仕様で表示がおかしくなるのでスペースをいれた

df_group <-
phenoData %>%
dplyr::select(characteristics_ch1.1, characteristics_ch1.2) %>%
set_colnames(c("genotype", "treatment")) %>%
mutate(genotype = str_replace_all(genotype, c("genotype: " = "")) %>% make.names()) %>%
mutate(treatment = str_replace_all(treatment, c("treatment: " = ""," 6 hours" = ""))) %>%
mutate(group = paste0(genotype, "_", treatment))
df_group

# 遺伝子発現データセットのオブジェクト作成

dge <-
DGEList(
counts =
dplyr::select(df_counts, str_detect(colnames(df_counts), pattern = "Sample") %>% which()) %>%
set_colnames(rownames(df_group)),
genes =
dplyr::pull(df_counts, ensembl_gene_id)
)

dge$samples$group <- as.factor(df_group$group)
dge$samples$genotype <- as.factor(df_group$genotype)
dge$samples$group <- as.factor(df_group$treatment)

# 遺伝子発現量の分布見てみる

p <-
cpm(dge, log = TRUE) %>%
as.data.frame() %>%
pivot_longer(everything(), names_to = "sample_id", values_to = "log_cpm") %>%
ggplot(mapping = aes(x = log_cpm, color = sample_id))

p +
geom_density() +
theme_bw()

# 低発現量の遺伝子を削ってNormalize。

dge <- dge[filterByExpr(dge), ]

dge <- calcNormFactors(dge)

# 対数化したカウントをlcpmとしてとってきて、PCAする

lcpm <- cpm(dge, log = TRUE)

pca <- prcomp(t(lcpm), scale = FALSE)
summary(pca)

# importanceをまとめる

df_importance <-
summary(pca)$importance %>% # $importanceでとってこれる
t() %>%
as.data.frame() %>%
set_colnames(colnames(.) %>% make.names()) %>%
mutate(PC = rownames(.) %>% as.factor())

df_importance

# 図示する。

p <-
df_importance %>%
ggplot(mapping = aes(x = PC))

p +
geom_bar(mapping = aes(y = Proportion.of.Variance), stat="identity", fill = "pink") +
geom_line(mapping = aes(y = Cumulative.Proportion, group = 1), color = "grey50") +
geom_point(mapping = aes(y = Cumulative.Proportion), color = "grey60") +
scale_x_discrete(limit = df_importance$PC) +
theme_bw()

# サンプル間の関係を次元圧縮して図示する。

p <-
pca$x %>% # 主成分得点をとってくる
as.data.frame() %>%
ggplot(mapping = aes(x = PC1, y = PC2))

p +
geom_point(mapping = aes(color = df_group$group), size = 3) +
labs(title = "PCA",
x = paste0("PC1"," ", df_importance[1, 2] %>% round(digits = 2) * 100, "%"),
y = paste0("PC2", " ", df_importance[2, 2] %>% round(digits = 2) * 100, "%")) +
theme_bw()

# ヒートマップを作成する前に、分散トップ200の遺伝子をとってくる

# 200の根拠は、今回は私のノートPCが悲鳴をあげない程度。。

var_rank <- lcpm %>% rowVars() %>% rank()
var_rank <- length(lcpm %>% rowVars()) - var_rank # 降順にランキングする方法がわからなかったので、強引にいった

dge_topvar <- dge[var_rank < 200, ] # フィルタリング

# ヒートマップを描いてみる

heat_plot <-
pheatmap(
mat = cpm(dge_topvar, log = TRUE),
col = brewer.pal(9, 'YlOrRd'),
cluster_rows = TRUE, cluster_cols = TRUE,
clustering_method = 'ward.D',
annotation_col = dplyr::select(df_group, -group),
angle_col = 90,
show_colnames = TRUE, show_rownames = TRUE,
main = "Heatmap")