0. 數(shù)據(jù)準(zhǔn)備
- 輸入數(shù)據(jù)集的要求：已經(jīng)進(jìn)行了如下分析的Seurat對(duì)象
- 導(dǎo)入演示數(shù)據(jù)
#官方演示數(shù)據(jù)集
wget https://swaruplab.bio.uci.edu/public_data/Zhou_2020.rds
seurat_obj <- readRDS('Zhou_2020.rds')
這是一個(gè)正常的腦組織數(shù)據(jù)集，包含了使用Harmony整合的12個(gè)樣本的Seurat對(duì)象。
- 加載R包和數(shù)據(jù)集
# single-cell analysis package
library(Seurat)
# plotting and data science packages
library(tidyverse)
library(cowplot)
library(patchwork)
# co-expression network analysis packages:
library(WGCNA)
library(hdWGCNA)
# using the cowplot theme for ggplot
theme_set(theme_cowplot())
# set random seed for reproducibility
set.seed(12345)
#可視化一下看看
p <- DimPlot(seurat_obj, group.by='cell_type', label=TRUE) +
umap_theme()
p
1. Set up Seurat object for WGCNA
在進(jìn)行hdWGCNA運(yùn)算之前，我們首先需要處理Seurat object。絕大多數(shù) hdWGCNA運(yùn)算所需要的信息都儲(chǔ)存在Seurat object’s @misc slot。一個(gè)Seurat object can hold multiple hdWGCNA experiments, for example representing different cell types in the same single-cell dataset. Notably, since we consider hdWGCNA to be a downstream data analysis step, we do not support subsetting the Seurat object after SetupForWGCNA has been run.
在這一步我們使用 SetupForWGCNA 來指定 hdWGNCA experiment的名稱。這一步指定了用于WGCNA分析的基因 。關(guān)于基因的選擇，有如下三種方法：
在這個(gè)演示里，我們選擇了至少在5%的細(xì)胞中表達(dá)的基因（3857個(gè)），并將hdWGCNA experiment 命名為 “tutorial”。
seurat_obj <- SetupForWGCNA(
seurat_obj,
gene_select = "fraction", # the gene selection approach
fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included
wgcna_name = "tutorial" # the name of the hdWGCNA experiment
)
length(seurat_obj@misc$tutorial$wgcna_genes)
# [1] 12639
2. Construct metacells
準(zhǔn)備好 Seurat 對(duì)象后，首先需要從單細(xì)胞數(shù)據(jù)集中構(gòu)建metacells。 簡(jiǎn)而言之，metacells是來自相同生物樣品的一小組相似細(xì)胞的集合。 k-最近鄰 (KNN) 算法可以識(shí)別相似細(xì)胞組，然后計(jì)算這些細(xì)胞的平均或總表達(dá)量，從而產(chǎn)生metacells基因表達(dá)矩陣。 metacells表達(dá)矩陣的稀疏性與原始表達(dá)矩陣相比大大降低，因此更適用于后續(xù)分析（WGCNA 等相關(guān)網(wǎng)絡(luò)方法對(duì)數(shù)據(jù)稀疏性很敏感）。（單細(xì)胞表觀基因組方法，例如 Cicero，在構(gòu)建可訪問網(wǎng)絡(luò)之前采用了類似的元細(xì)胞聚合方法。）
hdWGCNA 使用 MetacellsByGroups函數(shù)來構(gòu)造給定單細(xì)胞數(shù)據(jù)集的metacells表達(dá)矩陣。 此函數(shù)為存儲(chǔ)在 hdWGCNA experiment中的metacells數(shù)據(jù)集構(gòu)造一個(gè)新的 Seurat 對(duì)象。 group.by 參數(shù)確定將在哪些組中構(gòu)建metacells。我們只想從來自相同生物樣本的細(xì)胞構(gòu)建metacells，因此通過 group.by 參數(shù)將該信息傳遞給 hdWGCNA 至關(guān)重要。 此外，我們通常為每種細(xì)胞類型分別構(gòu)建metacells。 因此，在這個(gè)演示中，我們按 Sample 和 cell_type 進(jìn)行分組以達(dá)到預(yù)期的結(jié)果。
k值（The number of cells to be aggregated）應(yīng)根據(jù)輸入數(shù)據(jù)集的大小進(jìn)行調(diào)整，通常較小的 k 數(shù)可用于小數(shù)據(jù)集。 我們一般使用 20 到 75 之間的 k 值。本教程使用的數(shù)據(jù)集有 40,039 個(gè)細(xì)胞，每個(gè)生物樣本中的細(xì)胞數(shù)從 890 到 8,188 不等，這里我們使用 k=25。 可以使用 max_shared 參數(shù)調(diào)整元單元之間允許的重疊量。
注意：metacells聚合方法對(duì)extremely underrepresented cell types效果不好。 例如，在這個(gè)數(shù)據(jù)集中，腦血管細(xì)胞（周細(xì)胞和內(nèi)皮細(xì)胞）是the least represented，因此將它們排除在此分析之外。 MetacellsByGroups 有一個(gè)參數(shù) min_cells 來排除小于指定單元數(shù)的組。
# construct metacells in each group
seurat_obj <- MetacellsByGroups(
seurat_obj = seurat_obj,
group.by = c("cell_type", "Sample"), # specify the columns in seurat_obj@meta.data to group by
k = 25, # nearest-neighbors parameter
max_shared = 10, # maximum number of shared cells between two metacells
ident.group = 'cell_type' # set the Idents of the metacell seurat object
)
# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)
Optional: Process the Metacell Seurat Object
由于我們將 Metacell 表達(dá)信息存儲(chǔ)為單獨(dú)的 Seurat 對(duì)象，因此我們可以對(duì) Metacell 數(shù)據(jù)運(yùn)行 Seurat 函數(shù)。 我們可以使用 GetMetacellObject 從 hdWGCNA 實(shí)驗(yàn)中獲取元單元對(duì)象。
metacell_obj <- GetMetacellObject(seurat_obj)
Additionally, we have included a few wrapper functions to apply the Seurat workflow to the metacell object within the hdWGCNA experiment. Here we apply these wrapper functions to process the metacell object and visualize the aggregated expression profiles in two dimensions with UMAP.
seurat_obj <- NormalizeMetacells(seurat_obj)
seurat_obj <- ScaleMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunPCAMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunHarmonyMetacells(seurat_obj, group.by.vars='Sample')
seurat_obj <- RunUMAPMetacells(seurat_obj, reduction='harmony', dims=1:15)
這些結(jié)果儲(chǔ)存在Seurat object@misc slot
p1 <- DimPlotMetacells(seurat_obj, group.by='cell_type') + umap_theme() + ggtitle("Cell Type")
p2 <- DimPlotMetacells(seurat_obj, group.by='Sample') + umap_theme() + ggtitle("Sample")
p1 | p2
3. 共表達(dá)網(wǎng)絡(luò)分析
3.1 Set up the expression matrix
這一步將指定將用于網(wǎng)絡(luò)分析的表達(dá)矩陣，使用SetDatExpr 函數(shù)，用于存儲(chǔ)將用于下游網(wǎng)絡(luò)分析的給定單元組的轉(zhuǎn)置表達(dá)式矩陣。這一步默認(rèn)情況下使用元單元表達(dá)式矩陣 (use_metacells=TRUE)。
由于我們只想分析抑制性神經(jīng)元 (INH) ，因此使用group_name進(jìn)行了指定。如果想包含多個(gè)細(xì)胞類型，可以用group_name = c("INH", "EX")
seurat_obj <- SetDatExpr(
seurat_obj,
group_name = "INH", # the name of the group of interest in the group.by column
group.by='cell_type', # the metadata column containing the cell type info. This same column should have also been used in MetacellsByGroups
assay = 'RNA', # using RNA assay
slot = 'data' # using normalized data
)
3.2 挑選軟閾值
接下來我們將選擇“軟閾值”。這是 hdWGNCA 分析（以及普通 WGCNA）中極其重要的一步。 hdWGCNA 構(gòu)建基因-基因相關(guān)鄰接矩陣來推斷基因之間的共表達(dá)關(guān)系。將相關(guān)性提高到一個(gè)冪以減少相關(guān)矩陣中存在的噪聲，從而保留強(qiáng)連接并去除弱連接。因此，確定軟功率閾值的適當(dāng)值至關(guān)重要。
我們使用 TestSoftPowers 函數(shù)來為不同的軟閾值進(jìn)行parameter sweep。此功能通過檢查不同power值的結(jié)果網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)，幫助我們指導(dǎo)我們選擇構(gòu)建共表達(dá)網(wǎng)絡(luò)的軟閾值。共表達(dá)網(wǎng)絡(luò)應(yīng)該具有無標(biāo)度拓?fù)洌虼?TestSoftPowers 函數(shù)模擬了共表達(dá)網(wǎng)絡(luò)在不同軟功率閾值下與無標(biāo)度圖的相似程度。此外，我們可以使用函數(shù) PlotSoftPowers 來可視化parameter sweep的結(jié)果。
# Test different soft powers:
seurat_obj <- TestSoftPowers(
seurat_obj,
networkType = 'signed' # you can also use "unsigned" or "signed hybrid"
)
# plot the results:
plot_list <- PlotSoftPowers(seurat_obj)
# assemble with patchwork
wrap_plots(plot_list, ncol=2)
WGCNA 和 hdWGCNA 的一般要求是選擇無標(biāo)度拓?fù)淠Ｐ蛿M合大于或等于 0.8 的最低軟閾值，因此在這種情況下，我們選擇9為軟閾值。如果不設(shè)置軟閾值，ConstructNetwork 將 自動(dòng)選擇軟閾值。
Parameter sweep的輸出結(jié)果存儲(chǔ)在 hdWGCNA experiment中，可以使用 GetPowerTable 函數(shù)訪問以進(jìn)行進(jìn)一步檢查：
power_table <- GetPowerTable(seurat_obj)
head(power_table)
# Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
# 1 1 0.192964167 10.1742363 0.9508966 5887.9290 5901.7852 6501.0283
# 2 2 0.001621749 0.4801145 0.9828730 3092.6379 3092.7878 3835.2202
# 3 3 0.150812707 -3.2345892 0.9918393 1657.6198 1646.6169 2349.6967
# 4 4 0.407513393 -4.5383131 0.9898440 905.8402 890.9431 1486.7695
# 5 5 0.585949838 -4.8303680 0.9917528 504.4288 490.2683 968.1157
# 6 6 0.677668481 -4.6749295 0.9935437 286.1585 274.3233 646.6741
3.3 構(gòu)建共表達(dá)網(wǎng)絡(luò)
隨后我們就可以使用ConstructNetwork函數(shù)構(gòu)建共表達(dá)網(wǎng)絡(luò)了。（The parameters for blockwiseConsensusModules can be passed directly to ConstructNetwork with the same parameter names.）
# construct co-expression network:
seurat_obj <- ConstructNetwork(
seurat_obj, soft_power=9,
setDatExpr=FALSE,
tom_name = 'INH' # name of the topoligical overlap matrix written to disk
)
共表達(dá)網(wǎng)絡(luò)的可視化
PlotDendrogram(seurat_obj, main='INH hdWGCNA Dendrogram')
模塊組成基因
seurat_obj@misc$tutorial$wgcna_modules %>% head
# gene_name module color kME_grey kME_black kME_green kME_grey60 kME_red kME_lightcyan kME_tan kME_royalblue kME_brown kME_turquoise kME_lightyellow
# AL627309.1 AL627309.1 grey grey 0.06467554 0.01169450 0.01346971 -0.03798701 0.06389461 -0.0234879179 0.024303002 0.038220179 0.04439502 0.06965388 0.032234528
# LINC01409 LINC01409 black black 0.09562911 0.14258271 0.04532450 -0.03838167 0.07259315 0.0088322467 0.089630925 0.090604216 0.05090351 0.08046802 0.020385146
# LINC01128 LINC01128 grey grey 0.15025185 0.01456062 0.08874896 0.05263830 0.10848233 0.0233505152 0.036381392 0.003267192 0.03539629 0.13056570 0.034278458
# NOC2L NOC2L grey grey 0.16386360 -0.03203850 0.03588840 0.02706705 0.09822920 0.0223949415 0.044098042 0.013722673 0.02077281 0.15718774 0.007470014
# AGRN AGRN grey grey 0.06157292 0.05102518 0.08148903 -0.01092934 0.07061436 -0.0000888185 -0.011453627 0.028133157 0.06998927 0.04722421 0.043001474
# C1orf159 C1orf159 grey grey 0.03297087 0.02592184 0.02742151 -0.01461799 0.01558220 -0.0035206063 0.001851281 0.010550009 0.01190812 0.02965008 0.025278788
# kME_lightgreen kME_yellow kME_salmon kME_pink kME_magenta kME_purple kME_cyan kME_greenyellow kME_midnightblue kME_blue kME_darkred kME_darkgreen
# AL627309.1 0.04998851 0.03999694 0.05183194 0.057363638 0.02197360 0.07381381 0.022132715 0.0309476545 0.020243486 0.039376677 0.04946240 -0.041431912
# LINC01409 0.06062908 0.12675277 0.06172935 0.086179445 0.03545920 0.10836689 -0.010795315 0.1022885505 -0.013422943 0.008179931 0.12100397 -0.005445180
# LINC01128 0.12536877 0.04878530 0.09510907 0.037100032 0.06866021 0.07657537 0.108189685 0.0433716590 0.045969133 0.112251246 0.08318723 -0.045492473
# NOC2L 0.07886350 0.01475673 0.13338651 0.049125991 0.07227435 0.07841803 0.157399769 0.0232836501 0.025693226 0.166539003 0.03711394 -0.050635030
# AGRN 0.05589207 0.11086059 0.02339966 0.032878117 -0.03134549 0.05606373 0.036350264 0.0422577304 0.009418002 0.034792567 0.01463038 0.012060664
# C1orf159 0.01596504 0.03185948 0.02135547 0.004147534 0.00511424 0.02773708 -0.003693451 0.0004414052 0.005552397 -0.002849447 0.03548955 0.002986991
table(seurat_obj@misc$tutorial$wgcna_modules$module)
# grey black green grey60 red lightcyan tan royalblue brown turquoise lightyellow lightgreen
# 9597 154 188 87 166 88 129 67 230 281 76 86
# yellow salmon pink magenta purple cyan greenyellow midnightblue blue darkred darkgreen
# 203 123 151 150 149 123 131 112 232 60 56
3.4 Optional: inspect the topoligcal overlap matrix (TOM)
hdWGCNA represents the co-expression network as a topoligcal overlap matrix (TOM). This is a square matrix of genes by genes, where each value is the topoligcal overlap between the genes. The TOM is written to the disk when running ConstructNetwork, and we can load it into R using the GetTOM function. Advanced users may wish to inspect the TOM for custom downstream analyses.
TOM <- GetTOM(seurat_obj)
4. Module Eigengenes and Connectivity
在這一部分，我們將介紹如何計(jì)算單個(gè)細(xì)胞中的模塊特征基因，以及如何計(jì)算每個(gè)基因的基于特征基因的連通性。
4.1 Compute harmonized module eigengenes
Module Eigengenes (MEs) 是一種常用的指標(biāo)，用于總結(jié)整個(gè)共表達(dá)模塊的基因表達(dá)譜。簡(jiǎn)而言之，模塊特征基因是通過對(duì)包含每個(gè)模塊的基因表達(dá)矩陣的子集執(zhí)行主成分分析 (PCA) 來計(jì)算的。這些 PCA 矩陣中的第一個(gè)PC 就是 ME。
hdWGCNA 包含一個(gè)函數(shù) ModuleEigengenes 來計(jì)算單個(gè)單元格中的模塊特征基因。此外，我們?cè)试S用戶對(duì) ME 應(yīng)用 Harmony 批量校正，從而產(chǎn)生協(xié)調(diào)的模塊特征基因 (hME)。以下代碼使用 group.by.vars 參數(shù)執(zhí)行由原始樣本協(xié)調(diào)的模塊特征基因計(jì)算。
# need to run ScaleData first or else harmony throws an error:
seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))
# compute all MEs in the full single-cell dataset
seurat_obj <- ModuleEigengenes(
seurat_obj,
group.by.vars="Sample"
) #時(shí)間較長(zhǎng)
ME矩陣存儲(chǔ)為一個(gè)矩陣，其中每一行是一個(gè)細(xì)胞，每一列是一個(gè)模塊。 可以使用 GetMEs 函數(shù)從 Seurat 對(duì)象中提取此矩陣，該函數(shù)默認(rèn)檢索 hME。
# harmonized module eigengenes:
hMEs <- GetMEs(seurat_obj) #每個(gè)細(xì)胞對(duì)于每個(gè)模塊的特征值
head(hMEs)
# lightcyan darkgreen red turquoise salmon purple black darkred lightgreen royalblue tan greenyellow
# AGTCTTTGTTGATTCG-11 -3.728266 -2.410602 3.8018569 5.540111 6.002367 4.829816 5.308432 4.066452 3.087318 5.677753 6.637234 7.003169
# AGCAGCCTCCAGATCA-11 -3.242143 -2.725894 3.1256861 5.243863 5.191766 4.053458 5.157324 3.919190 4.276658 5.821838 6.023351 7.905414
# CTTAGGATCTCATTCA-11 -3.697165 -1.910883 4.2795263 5.803626 5.108230 3.798885 5.144346 3.710556 3.451550 5.056269 6.701767 8.761117
# AACTCAGAGGAATGGA-11 -3.438431 -3.141482 0.9050427 2.689078 2.172585 2.790341 7.279000 4.412214 2.103700 4.536073 5.719890 6.730259
# ACGAGGAGTCTCCACT-11 -3.752210 -2.695012 3.8369253 1.600519 2.965712 3.075344 5.888933 2.554901 2.433528 5.539103 8.523871 7.847494
# ACTTACTTCGGAGCAA-11 -3.193025 -1.252224 3.3207978 4.391861 4.744844 2.290026 4.946544 2.965951 4.713082 4.360662 6.201866 6.620763
# midnightblue pink magenta brown lightyellow yellow grey green grey60 cyan blue
# AGTCTTTGTTGATTCG-11 6.576508 5.369716 6.578784 7.616044 4.476879 7.459703 37.92130 7.224577 4.975482 13.5867856 14.6903483
# AGCAGCCTCCAGATCA-11 5.517671 4.143108 5.390045 8.251945 4.413418 7.395514 38.43114 8.578032 6.748117 13.6829303 12.4873768
# CTTAGGATCTCATTCA-11 5.156724 5.510369 4.768236 8.243442 4.932943 8.103645 36.48388 7.867925 7.018552 12.0219104 10.7403398
# AACTCAGAGGAATGGA-11 5.603454 6.998764 6.598624 9.449461 6.600172 7.963021 29.11238 8.023318 1.401920 0.4015516 0.9184507
# ACGAGGAGTCTCCACT-11 7.024574 7.239748 7.934361 9.010840 5.841042 8.041680 32.85802 7.300763 4.790319 4.6499158 3.8687891
# ACTTACTTCGGAGCAA-11 5.045800 4.115106 5.537796 8.347864 6.416239 6.914440 30.72694 7.490030 5.289845 9.2345490 10.2069815
# module eigengenes:
MEs <- GetMEs(seurat_obj, harmonized=FALSE)
4.2 Compute module connectivity
在共表達(dá)網(wǎng)絡(luò)分析中，我們經(jīng)常希望關(guān)注“hub gene”，即在每個(gè)模塊內(nèi)高度連接的那些基因。 因此，我們希望確定每個(gè)基因的eigengene-based connectivity，也稱為 kME。 hdWGCNA 使用 ModuleConnectivity 函數(shù)來計(jì)算完整單細(xì)胞數(shù)據(jù)集中的 kME 值，而不是元細(xì)胞數(shù)據(jù)集。 該函數(shù)本質(zhì)上計(jì)算基因和模塊特征基因之間的成對(duì)相關(guān)性。 kME可被用于該數(shù)據(jù)集中的所有細(xì)胞，但我們建議在用于計(jì)算 ConstructNetwork 的細(xì)胞類型或分組中計(jì)算 kME。
# compute eigengene-based connectivity (kME):
seurat_obj <- ModuleConnectivity(
seurat_obj,
group.by = 'cell_type', group_name = 'INH'
)
We can visualize the genes in each module ranked by kME using the PlotKMEs function.
# plot genes ranked by kME for each module
p <- PlotKMEs(seurat_obj, ncol=5)
p
4.3 Getting the module assignment table
hdWGCNA 允許使用 GetModules 函數(shù)輕松訪問module assignment table。 該表由三列組成：gene_name: 存儲(chǔ)基因的符號(hào)或 ID，module: 存儲(chǔ)基因的模塊分配，color: 存儲(chǔ)每個(gè)模塊的顏色映射，用于許多下游繪圖步驟。 如果在此 hdWGCNA 實(shí)驗(yàn)中調(diào)用了 ModuleConnectivity，則此表將具有每個(gè)模塊的 kME 的附加列。
# get the module assignment table:
modules <- GetModules(seurat_obj)
# show the first 6 columns:
head(modules[,1:6])
# gene_name module color kME_INH-M1 kME_INH-M2 kME_INH-M3
# AL627309.1 AL627309.1 grey grey -0.032349090 0.029917426 0.0379323320
# LINC01409 LINC01409 INH-M19 brown -0.045140924 -0.013473381 0.0102194168
# LINC01128 LINC01128 grey grey 0.050793988 0.109578251 0.1153173093
# NOC2L NOC2L grey grey 0.032490535 0.164524557 0.1699131451
# AGRN AGRN INH-M19 brown -0.008488577 0.035558532 0.0379326966
# C1orf159 C1orf159 grey grey -0.015618737 0.002235229 -0.0003824554
鑒定模塊內(nèi)hub基因
# get hub genes
hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)
head(hub_df)
# gene_name module kME
# 1 SEPTIN2 INH-M1 0.3019739
# 2 MICU2 INH-M1 0.3140922
# 3 WDR37 INH-M1 0.3204428
# 4 UBE2Q2L INH-M1 0.3297472
# 5 HSD17B4 INH-M1 0.3399253
# 6 ASB3 INH-M1 0.3417400
This wraps up the critical analysis steps for hdWGCNA, so remember to save your output.
saveRDS(seurat_obj, file='hdWGCNA_object.rds')
- 計(jì)算每個(gè)細(xì)胞對(duì)于每個(gè)模塊hub基因的表達(dá)活性(module score)
可使用seurat的AddModuleScore 函數(shù)或者Ucell算法
# compute gene scoring for the top 25 hub genes by kME for each module
# with Seurat method
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='Seurat'
)
# compute gene scoring for the top 25 hub genes by kME for each module
# with UCell method
library(UCell)
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='UCell'
)
5. Basic Visualization
5.1 Module Feature Plots
hdWGCNA includes the ModuleFeaturePlot function to consruct FeaturePlots for each co-expression module colored by each module’s uniquely assigned color.
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=6)
We can also plot the hub gene signature score using the same function:
# make a featureplot of hub scores for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='scores', # plot the hub gene scores
order='shuffle', # order so cells are shuffled
ucell = TRUE # depending on Seurat vs UCell for gene scoring
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=6)
5.2 Module Correlations
hdWGCNA includes the ModuleCorrelogram function to visualize the correlation between each module based on their hMEs, MEs, or hub gene scores using the R package corrplot.
ModuleCorrelogram(seurat_obj)
5.3 Seurat plotting functions
The base Seurat plotting functions are also great for visualizing hdWGCNA outputs. Here we demonstrate plotting hMEs using DotPlot and VlnPlot. The key to using Seurat’s plotting functions to visualize the hdWGCNA data is to add it into the Seurat object’s @meta.data slot:
# get hMEs from seurat object
MEs <- GetMEs(seurat_obj, harmonized=TRUE)
mods <- colnames(MEs); mods <- mods[mods != 'grey']
# add hMEs to Seurat meta-data:
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs)
Now we can easily use Seurat’s DotPlot function:
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'cell_type')
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
coord_flip() +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue')
# plot output
p
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'cell_type')
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
coord_flip() +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue')
# plot output
p
6. 相關(guān)文獻(xiàn)
- 內(nèi)毒素血癥中肝細(xì)胞轉(zhuǎn)錄改變引起巨噬細(xì)胞招募和T細(xì)胞抑制
- 心臟成纖維細(xì)胞通過Htra3-TGF-β-IGFBP7軸調(diào)節(jié)心衰的發(fā)展
參考：https://smorabit.github.io/hdWGCNA/articles/basic_tutorial.html
