Differentially-Expressed Gene Analysis

clp

12 June, 2020

准备工作

加载前面学习程中的R环境变量和必要的R包。

library(data.table)
library(ggplot2)
library(cowplot)
library(Seurat)

load('out/02_immune_cons.rd') #immune.combined

鉴定不同条件下差异表达的基因

现在我们已经对齐了刺激细胞和对照细胞，我们可以开始进行比较分析，看看刺激引起的差异。广泛观察这些变化的一种方法是绘制受刺激细胞和对照细胞的平均表达图，并在散点图上寻找可见的异常值的基因。在这里，我们采用刺激和对照幼稚T细胞和CD14单核细胞群体的平均表达，并生成散点图，突出显示对干扰素刺激表现出戏剧性反应的基因。

DefaultAssay(immune.combined) <- "RNA"

theme_set(theme_cowplot())
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
Idents(t.cells) <- "stim"
avg.t.cells <- log1p(AverageExpression(t.cells, verbose = FALSE)$RNA) #the default slot is .data which was normalized in the prev step!
avg.t.cells$gene <- rownames(avg.t.cells)

cd14.mono <- subset(immune.combined, idents = "CD14 Mono")
Idents(cd14.mono) <- "stim"
avg.cd14.mono <- log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA)
avg.cd14.mono$gene <- rownames(avg.cd14.mono)

genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = FALSE, colour="red")
#> Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
#> Please use `as_label()` or `as_name()` instead.
#> This warning is displayed once per session.
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")
p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = FALSE, colour="red")
plot_grid(p1, p2)

image.png

两个样本间的差异表达基因分析

正如你所看到的，许多相同的基因在这两种类型的细胞中都上调，很可能代表了一种保守的干素扰反应途径。

因为我们有信心在整个条件下识别出共同的细胞类型，所以我们可以提出问题，对于相同类型的细胞，在不同的条件下，哪些基因会发生变化。首先，我们在元数据槽(meta.data slot)中创建一列来保存细胞类型和刺激信息，并将当前标识切换到该列。然后，我们使用FindMarkers来寻找刺激和对照B细胞之间的差异基因。请注意，这里显示的许多排名靠前的基因与我们早先绘制的核心干扰素反应基因图相同。此外，像CXCL10这样的基因，我们看到的是单核细胞和B细胞干扰素反应的特异性基因，在这个列表中也显示出非常重要的意义。

immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")
Idents(immune.combined) <- "celltype.stim"

# If you want more rigorous DGE analysis, you can use RNA.count and ZINB (zero-inflation negative binomial) as a test model.
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)
head(b.interferon.response, n = 15)
#>                 p_val avg_logFC pct.1 pct.2     p_val_adj
#> ISG15   5.387767e-159 3.2053742 0.998 0.233 7.571429e-155
#> IFIT3   1.945114e-154 3.1250966 0.965 0.052 2.733468e-150
#> IFI6    2.503565e-152 2.9325755 0.965 0.076 3.518260e-148
#> ISG20   6.492570e-150 2.0373399 1.000 0.668 9.124009e-146
#> IFIT1   1.951022e-139 2.8645627 0.907 0.029 2.741772e-135
#> MX1     6.897626e-123 2.2571348 0.905 0.115 9.693234e-119
#> LY6E    2.825649e-120 2.1597818 0.898 0.150 3.970885e-116
#> TNFSF10 4.007285e-112 2.6378243 0.786 0.020 5.631437e-108
#> IFIT2   2.672552e-108 2.5226736 0.786 0.037 3.755738e-104
#> B2M      5.283684e-98 0.4219851 1.000 1.000  7.425161e-94
#> PLSCR1   4.634658e-96 1.9730793 0.793 0.113  6.513085e-92
#> IRF7     2.411149e-94 1.8012358 0.835 0.187  3.388388e-90
#> CXCL10   3.708508e-86 3.6838772 0.651 0.010  5.211566e-82
#> UBE2L6   5.547472e-83 1.4893574 0.851 0.297  7.795863e-79
#> PSMB9    1.716262e-77 1.1333816 0.937 0.568  2.411863e-73

Gene expression level in UMAP or Violin Plot

可视化这些基因表达变化的另一个有用的方法是使用选项split.by，即FeaturePlot或VlnPlot函数。这将显示给定基因列表的FeaturePlots，按分组变量(此处为刺激条件)拆分。CD3D和GNLY等基因是典型的细胞类型标记物(T细胞和NK/CD8T细胞)，几乎不受干扰素刺激的影响，在对照组和刺激组中表现出相似的基因表达模式。另一方面，IFI6和ISG15是核心干扰素反应基因，在所有细胞类型中都相应上调。最后，CD14和CXCL10是显示细胞类型特异性干扰素反应的基因。CD14单核细胞在刺激后表达下降，这可能导致监督分析框架中的错误分类，突显了综合分析的价值。CXCL10在干扰素刺激后，在单核细胞和B细胞中明显上调，但在其他类型的细胞中没有明显上调。

FeaturePlot(immune.combined, 
            features = c("CD3D", "GNLY", "IFI6"), 
            split.by = "stim", 
            max.cutoff = 3, 
            cols = c("grey", "red"))

image.png

plots <- VlnPlot(immune.combined, 
                 features = c("LYZ", "ISG15", "CXCL10"), 
                 split.by = "stim", 
                 group.by = "seurat_annotations", 
                 pt.size = 0, 
                 combine = FALSE)

CombinePlots(plots = plots, ncol = 1)
#> Warning: CombinePlots is being deprecated. Plots should now be combined
#> using the patchwork system.

image.png

还可以进一步探索

• Pseudo time analysis (Slingshot)拟时序分析
• ZINB DGE analysis (DEsingle)
• Functional enrichment analysis (GSEA)功能富集分析