library(Seurat)
## Registered S3 method overwritten by 'spatstat.geom':
## method from
## print.boxx cli
## Attaching SeuratObject
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(presto)
## Loading required package: Rcpp
## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:purrr':
##
## transpose
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.6.2
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(circlize)
## ========================================
## circlize version 0.4.14
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================
pbmc <- readRDS("/mnt/nectar_volume/home/eraz0001/KELLY 2020/E11.5/pbmc_E11.5_final.rds")
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
## Calculating cluster C1
## Calculating cluster C2
## Calculating cluster C3
## Calculating cluster C4
## Calculating cluster C5
## Calculating cluster C6
## Calculating cluster C7
## Calculating cluster C8
## Calculating cluster C9
## Calculating cluster C10
top200 = head(pbmc.markers$gene, 200)
top200
## [1] "Itm2a" "Capn6" "Cdkn1c" "Gpc3"
## [5] "Mfap4" "Sfrp2" "Runx1t1" "Fibin"
## [9] "Col3a1" "Tbx18" "Shox2" "Sfrp1"
## [13] "Col2a1" "Meg3" "6330403K07Rik" "Col9a2"
## [17] "Cpe" "Epha7" "Vcan" "Col9a3"
## [21] "Mfap2" "Tox" "Lect1" "Lgals1"
## [25] "Eid1" "Eif3e" "Matn4" "Osr2"
## [29] "Ldhb" "Cxcl12" "Ftl1" "Tbx15"
## [33] "Rcn3" "Tgfb2" "Mest" "Crabp2"
## [37] "Igf2" "Lhfp" "Wwp2" "Mab21l2"
## [41] "Tsc22d1" "Pkdcc" "Pitx1" "Col1a2"
## [45] "Gas2" "Dact1" "Gas1" "Robo2"
## [49] "Hes1" "Peg3" "Maged1" "Nnat"
## [53] "Mdk" "Hoxc10" "Pbx1" "Leprel2"
## [57] "Sox4" "Zfhx4" "Hmga2" "Grb10"
## [61] "Pebp1" "Scarf2" "Tcf7l2" "Cd63"
## [65] "Ttc3" "Aldh2" "Samd14" "Mxd4"
## [69] "Eif3f" "Sepw1" "Asb4" "Nsg1"
## [73] "Gstm7" "Gpc6" "Six1" "Meis1"
## [77] "Sox9" "Foxp2" "Cthrc1" "H19"
## [81] "Meis2" "Lrig3" "Rarg" "Fstl1"
## [85] "Eif4a2" "Nrk" "Enpp2" "Col27a1"
## [89] "Aes" "Arf4" "Creb5" "Sox5"
## [93] "2810055G20Rik" "Rhoc" "Snhg18" "Ptov1"
## [97] "Gsn" "Ndn" "Mcm6" "Wbp5"
## [101] "Sh3bgrl" "Serpinf1" "Fgfr1" "Mex3b"
## [105] "Col11a1" "Tspan3" "Peg10" "Tshz2"
## [109] "Pfn2" "Igsf3" "Malat1" "Ak1"
## [113] "Marcks" "Pmp22" "Igfbp5" "Rcn1"
## [117] "Pgrmc1" "Atxn7l3b" "Map1lc3a" "Arhgef25"
## [121] "Col25a1" "Papss1" "Gabarap" "Cyr61"
## [125] "Mageh1" "Id2" "Acot1" "Nedd4"
## [129] "Fzd2" "Dap" "Dtx3" "Sept6"
## [133] "Chd3" "Ugdh" "Scx" "Cd24a"
## [137] "Amotl2" "Mpped2" "Tgfbi" "Nktr"
## [141] "Mtch1" "Sox11" "Kdelr1" "Rian"
## [145] "Col9a1" "Twist2" "Trps1" "Efemp2"
## [149] "Jam3" "Ccnd2" "Jun" "Tmem263"
## [153] "Mmp11" "Bnip3l" "Clip3" "Dnm1"
## [157] "Mbnl2" "Ccdc80" "Ypel3" "Ghr"
## [161] "Podxl2" "Fgfrl1" "Il11ra1" "Hoxd9"
## [165] "Tead2" "2700089E24Rik" "Lix1" "Aplp1"
## [169] "Eno3" "Uchl1" "Frmd4a" "Celf2"
## [173] "Mcm2" "Smim14" "Ssr3" "Fbn2"
## [177] "Lpar1" "Mecom" "Grina" "Hoxa10"
## [181] "Rab34" "Gas5" "Tnfaip6" "Evl"
## [185] "Synpo" "Tcp11l2" "Tmem167" "Cdh2"
## [189] "Nek6" "Dab2" "Col1a1" "Xrcc5"
## [193] "Rcn2" "Frzb" "Jund" "Cdo1"
## [197] "Ccnl2" "Ogfr" "Sirt2" "Phactr2"
p <- DotPlot(object = pbmc, features = top200)
p
df<- p$data
exp_mat<-df %>%
dplyr::select(-pct.exp, -avg.exp) %>%
pivot_wider(names_from = id, values_from = avg.exp.scaled) %>%
as.data.frame()
row.names(exp_mat) <- exp_mat$features.plot
exp_mat <- exp_mat[,-1] %>% as.matrix()
head(exp_mat)
## C1 C2 C3 C4 C5 C6
## Itm2a 2.474792 -0.5915715 -0.5846374 0.9833787 -0.43956429 -0.4914915
## Capn6 1.800184 -0.3088892 -0.9512784 0.9399522 0.10904952 -0.9154756
## Cdkn1c 1.773028 -0.4505345 -0.7644985 0.7334957 0.27150138 -0.8714400
## Gpc3 1.611815 -0.3257414 -1.3003161 0.7587943 0.13689252 -0.2907668
## Mfap4 1.271478 -0.1255730 -1.0279182 0.6771833 0.21448980 -0.8258094
## Sfrp2 2.196049 -0.6240812 -0.7757256 0.9296428 -0.09148524 -0.7488678
## C7 C8 C9 C10
## Itm2a -0.23589316 0.08376446 -0.6171997 -0.5815780
## Capn6 0.09414872 1.11512517 -0.9236566 -0.9591598
## Cdkn1c 0.23728233 1.18376686 -0.9905409 -1.1220600
## Gpc3 0.05733375 1.33695143 -0.7049859 -1.2799764
## Mfap4 0.12875346 1.73507778 -1.0033164 -1.0443653
## Sfrp2 -0.17108688 0.77935219 -0.7335596 -0.7602380
percent_mat<-df %>%
dplyr::select(-avg.exp, -avg.exp.scaled) %>%
pivot_wider(names_from = id, values_from = pct.exp) %>%
as.data.frame()
row.names(percent_mat) <- percent_mat$features.plot
percent_mat <- percent_mat[,-1] %>% as.matrix()
head(percent_mat)
## C1 C2 C3 C4 C5 C6 C7
## Itm2a 91.58654 48.00000 49.529781 92.62821 71.31783 65.82915 57.78894
## Capn6 96.15385 83.46667 7.210031 94.87179 86.82171 15.57789 80.40201
## Cdkn1c 100.00000 98.40000 94.043887 99.35897 99.61240 89.44724 98.49246
## Gpc3 99.03846 97.86667 15.987461 99.67949 95.73643 93.46734 93.46734
## Mfap4 97.83654 92.00000 9.090909 98.39744 91.47287 43.21608 82.41206
## Sfrp2 81.97115 52.00000 5.642633 79.80769 71.31783 16.08040 64.82412
## C8 C9 C10
## Itm2a 72.17391 30.392157 47.368421
## Capn6 91.30435 8.823529 10.526316
## Cdkn1c 100.00000 85.294118 52.631579
## Gpc3 100.00000 91.176471 28.947368
## Mfap4 98.26087 14.705882 7.894737
## Sfrp2 80.00000 11.764706 13.157895
range(percent_mat)
## [1] 0 100
dim(exp_mat)
## [1] 200 10
dim(percent_mat)
## [1] 200 10
library(viridis)
## Loading required package: viridisLite
library(Polychrome)
## get an idea of the ranges of the matrix
quantile(exp_mat, c(0.1, 0.5, 0.9, 0.99))
## 10% 50% 90% 99%
## -1.0954770 -0.1407776 1.3706449 2.1905934
## any value that is greater than 2 will be mapped to yellow
col_fun = circlize::colorRamp2(c(-1, 0, 2), viridis(20)[c(1,10, 20)])
cell_fun = function(j, i, x, y, w, h, fill){
grid.rect(x = x, y = y, width = w, height = h,
gp = gpar(col = NA, fill = NA))
grid.circle(x=x,y=y,r= percent_mat[i, j]/100 * min(unit.c(w, h)),
gp = gpar(fill = col_fun(exp_mat[i, j]), col = NA))}
Heatmap(exp_mat,
heatmap_legend_param=list(title="expression"),
column_title = "clustered dotplot",
col=col_fun,
rect_gp = gpar(type = "none"),
cell_fun = cell_fun,
row_names_gp = gpar(fontsize = 5),
row_km = 4,
border = "black")
colnames(exp_mat)
## [1] "C1" "C2" "C3" "C4" "C5" "C6" "C7" "C8" "C9" "C10"
library(RColorBrewer)
cluster <- c("C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8", "C9", "C10")
column_ha<- HeatmapAnnotation(
cluster = cluster,
col = list(cluster = setNames(brewer.pal(10, "Paired"), unique(cluster))
),
na_col = "grey"
)
layer_fun = function(j, i, x, y, w, h, fill){
grid.rect(x = x, y = y, width = w, height = h,
gp = gpar(col = NA, fill = NA))
grid.circle(x=x,y=y,r= pindex(percent_mat, i, j)/100 * unit(2, "mm"),
gp = gpar(fill = col_fun(pindex(exp_mat, i, j)), col = NA))}
lgd_list = list(
Legend( labels = c(0,0.25,0.5,0.75,1), title = "pt",
graphics = list(
function(x, y, w, h) grid.circle(x = x, y = y, r = 0 * unit(2, "mm"),
gp = gpar(fill = "black")),
function(x, y, w, h) grid.circle(x = x, y = y, r = 0.25 * unit(2, "mm"),
gp = gpar(fill = "black")),
function(x, y, w, h) grid.circle(x = x, y = y, r = 0.5 * unit(2, "mm"),
gp = gpar(fill = "black")),
function(x, y, w, h) grid.circle(x = x, y = y, r = 0.75 * unit(2, "mm"),
gp = gpar(fill = "black")),
function(x, y, w, h) grid.circle(x = x, y = y, r = 1 * unit(2, "mm"),
gp = gpar(fill = "black")))
))
set.seed(123)
hp<- Heatmap(exp_mat,
heatmap_legend_param=list(title="expression"),
column_title = "clustered dotplot",
col=col_fun,
rect_gp = gpar(type = "none"),
layer_fun = layer_fun,
row_names_gp = gpar(fontsize = 5),
row_km = 4,
border = "black",
top_annotation = column_ha)
draw( hp, annotation_legend_list = lgd_list)
