說明：僅根據(jù)官網(wǎng)指南加個人理解消返，目前官網(wǎng)上發(fā)布了3.0版本的指南载弄，下面內(nèi)容是2.4版本耘拇，內(nèi)容上有些許出入。

整合分析的目的：

1. 鑒定在兩種數(shù)據(jù)集中都存在的細(xì)胞類型
2. 得到在實驗組和對照組中都存在細(xì)胞類型markers
3. 將兩種數(shù)據(jù)集進行對比宇攻，找到對刺激特異性反應(yīng)的細(xì)胞類型

1.構(gòu)建Seurat對象

library(Seurat)
#導(dǎo)入數(shù)據(jù)
ctrl.data <- read.table("~/Downloads/immune_control_expression_matrix.txt.gz", sep = "\t")
stim.data <- read.table("~/Downloads/immune_stimulated_expression_matrix.txt.gz", sep = "\t")

#對control對象進行處理
#創(chuàng)建Seurat對象惫叛，min.cells參數(shù)控制一個基因至少在5個細(xì)胞中被檢測到
ctrl <- CreateSeuratObject(raw.data = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)

#對ctrl對象設(shè)定標(biāo)簽，后面作圖的時候會用到
ctrl@meta.data$stim <- "CTRL"

#根據(jù)基因的表達過濾細(xì)胞逞刷，篩選出至少有00個基因的細(xì)胞
ctrl <- FilterCells(ctrl, subset.names = "nGene", low.thresholds = 500, high.thresholds = Inf)

#數(shù)據(jù)標(biāo)準(zhǔn)化
ctrl <- NormalizeData(ctrl)
ctrl <- ScaleData(ctrl, display.progress = F)

#對stim組進行相同的處理
stim <- CreateSeuratObject(raw.data = stim.data, project = "IMMUNE_STIM", min.cells = 5)
stim@meta.data$stim <- "STIM"
stim <- FilterCells(stim, subset.names = "nGene", low.thresholds = 500, high.thresholds = Inf)
stim <- NormalizeData(stim)
stim <- ScaleData(stim, display.progress = F)

#挑選差異基因進行后續(xù)的CCA分析 
ctrl <- FindVariableGenes(ctrl, do.plot = F) 
stim <- FindVariableGenes(stim, do.plot = F) 
g.1 <- head(rownames(ctrl@hvg.info), 1000) 
g.2 <- head(rownames(stim@hvg.info), 1000) 
genes.use <- unique(c(g.1, g.2)) 
genes.use <- intersect(genes.use, rownames(ctrl@scale.data)) 
genes.use <- intersect(genes.use, rownames(stim@scale.data))

2.進行CCA分析

#CCA類似于PCA嘉涌，降維的一種方法
#RunCCA這一個函數(shù)將兩個對象整合成一個對象
#num.cc指的是Number of canonical vectors to calculate
immune.combined <- RunCCA(ctrl, stim, genes.use = genes.use, num.cc = 30)
#可視化CCA結(jié)果
#reduction.use參數(shù)表示選用哪一種降維的方法，默認(rèn)是pca亲桥，還可以選擇tsne洛心、ica固耘、cca
#features.plot參數(shù)表示用來展示的特征（這里是CC1）题篷，可以是基因表達、pc分?jǐn)?shù)等
p1 <- DimPlot(object = immune.combined, reduction.use = "cca", group.by = "stim", pt.size = 0.5, do.return = TRUE)
p2 <- VlnPlot(object = immune.combined, features.plot = "CC1", group.by = "stim", do.return = TRUE)
plot_grid(p1, p2)

#這一步可以查看定義特定CC的基因
PrintDim(object = immune.combined, reduction.type = "cca", dims.print = 1:2, genes.print = 10)

#根據(jù)曲線的走勢確定
#根據(jù)下方的圖厅目，選定CC1-20進行分析
p3 <- MetageneBicorPlot(immune.combined, grouping.var = "stim", dims.eval = 1:30, display.progress = FALSE)

#用DimHeatmap函數(shù)查看每一個CC中的top gene對細(xì)胞分群的影響
#這里查看的是9個CC
DimHeatmap(object = immune.combined, reduction.type = "cca", cells.use = 500, dim.use = 1:9, do.balanced = TRUE)

4.整合分析

#reduction.use參數(shù)指的是選擇降維方法番枚，默認(rèn)是pca
#dims.use參數(shù)指的是上述選擇的CC個數(shù)
#resolution參數(shù)指的是分辨率，越大群分的越細(xì)
immune.combined <- RunTSNE(immune.combined, reduction.use = "cca.aligned", dims.use = 1:20, do.fast = T)
immune.combined <- FindClusters(immune.combined, reduction.type = "cca.aligned", resolution = 0.6, dims.use = 1:20)

# 可視化
p1 <- TSNEPlot(immune.combined, do.return = T, pt.size = 0.5, group.by = "stim")
p2 <- TSNEPlot(immune.combined, do.label = T, do.return = T, pt.size = 0.5)
plot_grid(p1, p2)

5.鑒定保守的細(xì)胞類型marker

#這一步是鑒定在不同條件下都保守的細(xì)胞類型marker
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 7, grouping.var = "stim", print.bar = FALSE)
head(nk.markers)

#查看特定基因在細(xì)胞群的表達情況
FeaturePlot(object = immune.combined, features.plot = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A", "CCL2", "PPBP"), min.cutoff = "q9", cols.use = c("lightgrey", "blue"), pt.size = 0.5)

#在所有群當(dāng)中查看特定基因的表達和占比情況
#確定細(xì)胞群類型
immune.combined@ident <- factor(immune.combined@ident, levels = (c("pDC", "Eryth", "Mk", "DC", "CD14 Mono", "CD16 Mono", "B activated", "B", "CD8 T", "NK", "T activated", "CD4 Naive T", "CD4 Memory T")))
#確定基因
markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5", "CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1", "GPR183", "PPBP", "GNG11", "HBA2", "HBB", "TSPAN13", "IL3RA", "IGJ")
#出圖
sdp <- SplitDotPlotGG(immune.combined, genes.plot = rev(markers.to.plot), cols.use = c("blue", "red"), x.lab.rot = T, plot.legend = T, dot.scale = 8, do.return = T, grouping.var = "stim")

6.鑒定在不同條件下差異表達的基因

方法一：

#一種方法就是分別查看實驗組和對照組的基因表達情況的平均值损敷，以此為依據(jù)葫笼，鑒定出離群的基因
#這一部分主要是構(gòu)建函數(shù)用于后續(xù)分析的使用
LabelPoint <- function(plot, genes, exp.mat, adj.x.t = 0, adj.y.t = 0, adj.x.s = 0, 
    adj.y.s = 0, text.size = 2.5, segment.size = 0.1) {
    for (i in genes) {
        x1 <- exp.mat[i, 1]
        y1 <- exp.mat[i, 2]
        plot <- plot + annotate("text", x = x1 + adj.x.t, y = y1 + adj.y.t, 
            label = i, size = text.size)
        plot <- plot + annotate("segment", x = x1 + adj.x.s, xend = x1, y = y1 + 
            adj.y.s, yend = y1, size = segment.size)
    }
    return(plot)
}

LabelUR <- function(plot, genes, exp.mat, adj.u.t = 0.1, adj.r.t = 0.15, adj.u.s = 0.05, 
    adj.r.s = 0.05, ...) {
    return(LabelPoint(plot, genes, exp.mat, adj.y.t = adj.u.t, adj.x.t = adj.r.t, 
        adj.y.s = adj.u.s, adj.x.s = adj.r.s, ...))
}

LabelUL <- function(plot, genes, exp.mat, adj.u.t = 0.1, adj.l.t = 0.15, adj.u.s = 0.05, 
    adj.l.s = 0.05, ...) {
    return(LabelPoint(plot, genes, exp.mat, adj.y.t = adj.u.t, adj.x.t = -adj.l.t, 
        adj.y.s = adj.u.s, adj.x.s = -adj.l.s, ...))
}

#挑選出實驗組和對照組中特定的細(xì)胞亞群(這里是CD4 Naive T和CD14 Mono)
#算出這一亞群基因表達的平均值，標(biāo)記出已知的對實驗敏感的基因
t.cells <- SubsetData(immune.combined, ident.use = "CD4 Naive T", subset.raw = T)
t.cells <- SetAllIdent(t.cells, id = "stim")
avg.t.cells <- log1p(AverageExpression(t.cells, show.progress = FALSE))
avg.t.cells$gene <- rownames(avg.t.cells)

cd14.mono <- SubsetData(immune.combined, ident.use = "CD14 Mono", subset.raw = T)
cd14.mono <- SetAllIdent(cd14.mono, id = "stim")
avg.cd14.mono <- log1p(AverageExpression(cd14.mono, show.progress = FALSE))
avg.cd14.mono$gene <- rownames(avg.cd14.mono)

genes.to.label1 = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1")
genes.to.label2 = c("IFIT2", "IFIT1")
genes.to.label3 = c("CXCL10", "CCL8")
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")
p1 <- LabelUR(p1, genes = c(genes.to.label1, genes.to.label2), avg.t.cells, 
    adj.u.t = 0.3, adj.u.s = 0.23)
p1 <- LabelUL(p1, genes = genes.to.label3, avg.t.cells, adj.u.t = 0.5, adj.u.s = 0.4, 
    adj.l.t = 0.25, adj.l.s = 0.25)
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")
p2 <- LabelUR(p2, genes = c(genes.to.label1, genes.to.label3), avg.cd14.mono, 
    adj.u.t = 0.3, adj.u.s = 0.23)
p2 <- LabelUL(p2, genes = genes.to.label2, avg.cd14.mono, adj.u.t = 0.5, adj.u.s = 0.4, 
    adj.l.t = 0.25, adj.l.s = 0.25)
plot_grid(p1, p2)

#這一步就是用于查找在不同處理條件下差異表達的基因
immune.combined@meta.data$celltype.stim <- paste0(immune.combined@ident, "_", immune.combined@meta.data$stim)
immune.combined <- StashIdent(immune.combined, save.name = "celltype")
immune.combined <- SetAllIdent(immune.combined, id = "celltype.stim")
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", print.bar = FALSE)
head(b.interferon.response, 15)

方法二：

#將會展示給定基因在不同處理條件下的表達情況
FeatureHeatmap(immune.combined, features.plot = c("CD3D", "GNLY", "IFI6", "ISG15", "CD14", "CXCL10"), group.by = "stim", pt.size = 0.25, key.position = "top", max.exp = 3)