Zenn
🐕

【Plot】[ComplexHeatmap] UpSet Plotから組み合わせパターンごとに遺伝子を取り出す

2023/05/09に公開
R
bioinformatics
seurat
complexheatmap
upsetplot
tech

こちらの記事で書いたUpSet plotのデータを使用する。リンク元では細胞集合 (クラスター)に特徴的な遺伝子が他のクラスターとどれだけ重複しているか、独立しているのかをComplexHeatmapパッケージのUpSet plotで図示した。
https://zenn.dev/rchiji/articles/0cdf3fec190c31

組み合わせパターンごとの要素を取り出す

UpSet plotで可視化した後は、各パターン（1つのクラスターだけのもの、cluster2と3で共通しているもの、、など）ごとの遺伝子が確認したくなる。

1つのクラスターだけのものを取り出すのは次のように、deg_listのi番目 vs その他でsetdiffをしていけば良いが、その他パターンのものはコードで書くには煩雑である。

# 空のリストを用意
unique_deglist <- list()

for(i in seq(deg_list)){
  # deg_listのi番目とdeg_listのi番目以外をsetdiffして、結果をリストに保存
  unique_deglist[[i]] <- Reduce(f = setdiff, x = c(deg_list[i], deg_list[-i]))
}

names(unique_deglist) <- names(deg_list)


extract_comb

ComplexHeatmapではextract_comb()機能を使って特定のパターンの要素を取得することができる。その際にcomb_name=引数で1つのパターンを指定するが、comb_nameは01の2進数コードを指定する必要がある。

パターンごとの2進数コードはcomb_name機能で確認できる。

comb_name(mt)

readable=TRUEにすると2進数コードを可読性のあるものに変換してくれる。それでも見づらいが。

パターンコードを指定して要素を取り出してみる。"(!cluster_0)&cluster_3&cluster_1&(!cluster_2)&cluster_4&cluster_5&(!cluster_6)"の箇所である。つまりcluster3, 1, 4, 5に共通している遺伝子を見ている。

extract_comb(m = mt, comb_name = "0110110")


全てのパターンごとに要素を取り出す

パターンコードごとに遺伝子をまとめる例を示す。

# パターンコード一覧を取得
patterns <- comb_name(m = mt)

# 空のリストを用意。名前も付けておく
extract_genes <- vector(mode = "list", length = length(patterns))
names(extract_genes) <- patterns　# names(extract_genes) <- comb_name(m = mt, readable = T)でもいいかも

# パターンコードごとにextract_combで取り出して、リストに保存
for(i in patterns){
  extract_genes[[i]] <- extract_comb(m = mt, comb_name = i)
}


さらにリストを要素の長さでソートすると扱いやすい。これをunlist()してやればパターン順に並び替えた遺伝子ベクトルが作成できる。

# sapplyでリストにlength関数を繰り返し処理した結果を、orderで降順に並び替え
extract_genes <- extract_genes[order(sapply(extract_genes, length), decreasing = T)]



ちなみに次のようにしてやれば各遺伝子とパターンコード、他との重複数をデータフレームにまとめることができる。

# 遺伝子名とパターンコードについてのデータフレームを作成する機能をリストに繰り返し処理
tmp <- mapply(function(x, y){
  df <- data.frame(genes = x,  # 遺伝子ベクトルが入ったリストから、遺伝子ベクトルを取り出し
                   code = rep(y, times = length(x)) # リストの名前 (パターンコード)を遺伝子ベクトルの要素数だけ繰り返し
  )
  return(df)
  },
  x = extract_genes,
  y = names(extract_genes),
  SIMPLIFY = F # 結果をリストとして出力させる
  )

# リストの要素をrbindで縦につなげる
genes_df <- do.call(rbind, tmp)
rownames(genes_df) <- genes_df$genes

# str_counts関数でパターンコード中の1の数を調べる。
genes_df$num_overlap <- stringr::str_count(string = genes_df$code, pattern = "1")


他のクラスターとの重複度をヒートマップアノテーションに与える

各クラスターの上位DEGがどれだけクラスターに特異的なのか、UpSet plotのパターンコードの情報も使って可視化する例を示す。

# Heatmap用に発現データをスケーリングしておく
seurat_object <- ScaleData(object = seurat_object, features = rownames(seurat_object))

# unique degが多いクラスターの順番を取得
unique_deglist <- list()
for(i in seq(deg_list)){
  unique_deglist[[i]] <- Reduce(f = setdiff, x = c(deg_list[i], deg_list[-i]))
}
names(unique_deglist) <- names(deg_list)
cluster_order <- names(unique_deglist[order(sapply(unique_deglist, length), decreasing = T)])

## Annotation data
anno <- seurat_object@meta.data
anno <- anno[,"Cluster", drop = FALSE]

## color setting
library(RColorBrewer)

cluster_size <- length(unique(anno$Cluster))
cluster_color <- scales::hue_pal()(cluster_size)
names(cluster_color) <- unique(anno$Cluster)

col_fun = circlize::colorRamp2(breaks = seq(cluster_size),
                               colors = rev(brewer.pal(n = cluster_size, name = "Spectral"))
                               )

## column order (unique DEGが多いクラスター順)
anno$Cluster <- factor(anno$Cluster, levels = cluster_order)
anno <- arrange(anno, Cluster)

## Expression matrix
mat <- pbmc@assays$RNA@scale.data[rownames(genes_df), rownames(anno)]
mat <- as.matrix(mat)

## culated_gene label position
# 各クラスターDEGの平均log2FCの上位5個ずつを見てみる。
curated_genes <- 
  deg %>%
  group_by(cluster) %>%
  top_n(n = 5, wt = avg_log2FC)
curated_genes <- unique(curated_genes$gene)

lab_position <- which(rownames(mat) %in% curated_genes)
lab_final <- rownames(mat)[lab_position]

## Annotation object
rowAnno <- rowAnnotation(num_overlap = anno_simple(genes_df$num_overlap,
                                                   col = col_fun,
                                                   width = unit(2,"mm"),
                                                   ),
                         show_annotation_name = FALSE
                         )

topAnno <- HeatmapAnnotation(df = anno,
                             show_annotation_name = FALSE,
                             col = list(Cluster = cluster_color))

## Heatmap
ht <- Heatmap(matrix = mat,
              cluster_rows = F,
              cluster_columns = F,
              show_column_names = F,
              show_row_names = F,
              col = circlize::colorRamp2(breaks = c(-3,0,3), colors = c("magenta","black","yellow")),
              heatmap_legend_param = list(title = "Scaled expression", at = seq(-3, 3, 1)),
              top_annotation = topAnno,
              left_annotation = rowAnno,
              column_split = anno$seurat_clusters,
              use_raster = TRUE,
              raster_quality = 5,
              raster_by_magick = FALSE, # --> magickっていうパッケージがあると、ラスタライズされすぎる。
              )

lgd <- Legend(title = "No. overlap",
              col_fun = col_fun,
              at = seq(0, cluster_size, by = 2))

draw(ht, annotation_legend_list = list(lgd))

top DEGだからと言ってクラスター特異的とは限らないのが分かる。

Discussion