劉小澤寫于19.10.17
筆記目的：根據(jù)生信技能樹(shù)的單細(xì)胞轉(zhuǎn)錄組課程探索10X Genomics技術(shù)相關(guān)的分析
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二單元第9講：細(xì)胞亞群之間的差異分析
前言
這次的任務(wù)是模仿原文的：
第五張：比較兩種CD8+細(xì)胞差異
在分群結(jié)果可以看到，CD8+主要分成了兩群烹植，一個(gè)是紅色的（170個(gè)CD8+ cytotoxic T cells看疗，即細(xì)胞毒性T細(xì)胞）女揭，一個(gè)是淺藍(lán)色的（429個(gè)CD8+ effector T cells永罚，即效應(yīng)T細(xì)胞）
【圖片注釋：Heat map of selected significantly differentially expressed genes comparing CD8+ T cells in the red activated cluster (n = 170) to those in the blue effector/EM cluster (n = 429) at response (day + 376) 】
首先進(jìn)行差異分析
準(zhǔn)備數(shù)據(jù)
=> SubsetData()取子集
既然是分析的day +376數(shù)據(jù)悯森，那么就先把這一小部分?jǐn)?shù)據(jù)提取出來(lái)：
rm(list = ls())
options(warn=-1)
suppressMessages(library(Seurat))
load('./patient1.PBMC.output.Rdata')
PBMC_RespD376 = SubsetData(PBMC,TimePoints =='PBMC_RespD376')
> table(PBMC_RespD376@ident)
0 1 2 3 4 5 6 7 8 9 10 11 12
800 433 555 677 636 516 119 324 204 200 170 11 39
然后這個(gè)圖是分析了紅色和淺藍(lán)色的兩群荐糜，結(jié)合之前得到的分群結(jié)果巷怜，紅色是第10群，淺藍(lán)色是第4群
于是再用這個(gè)函數(shù)提取出來(lái)第4狞尔、10群
PBMC_RespD376_for_DEG = SubsetData(PBMC_RespD376,
PBMC_RespD376@ident %in% c(4,10))
利用monocle V2構(gòu)建對(duì)象
=> newCellDataSet()
我們需要三樣?xùn)|西：表達(dá)矩陣丛版、細(xì)胞信息、基因信息
首先來(lái)看表達(dá)矩陣：上面取到的PBMC_RespD376_for_DEG中包含了多個(gè)數(shù)據(jù)接口偏序，單是表達(dá)矩陣相關(guān)就有三個(gè)
關(guān)于這三者的不同：
We view
object@raw.dataas a storage slot for the non-normalized data, upon which we perform normalization and returnobject@data.（來(lái)自：https://github.com/satijalab/seurat/issues/351）也就是說(shuō)页畦，raw.data是最原始的矩陣，data是歸一化之后的研儒，scale.data是標(biāo)準(zhǔn)化后的（一般是z-score處理）
data：The normalized expression matrix (log-scale)scale.data：scaled (default is z-scoring each gene) expression matrix; used for dimmensional reduction and heatmap visualization（來(lái)自：seurat的各個(gè)接口含義）
我們使用log處理過(guò)的歸一化表達(dá)矩陣
count_matrix=PBMC_RespD376_for_DEG@data
> dim(count_matrix)
[1] 17712 806
然后豫缨，看細(xì)胞信息
# 細(xì)胞分群信息
cluster=PBMC_RespD376_for_DEG@ident
> table(cluster)
cluster
4 10
636 170
# 細(xì)胞名稱（barcode名稱）
names(count_matrix)
最后是基因信息
gene_annotation <- as.data.frame(rownames(count_matrix))
萬(wàn)事俱備，開(kāi)始monocle
library(monocle)
> packageVersion('monocle')
[1] ‘2.12.0’
# 1.表達(dá)矩陣
expr_matrix <- as.matrix(count_matrix)
# 2.細(xì)胞信息
sample_ann <- data.frame(cells=names(count_matrix),
cellType=cluster)
rownames(sample_ann)<- names(count_matrix)
# 3.基因信息
gene_ann <- as.data.frame(rownames(count_matrix))
rownames(gene_ann)<- rownames(count_matrix)
colnames(gene_ann)<- "genes"
# 然后轉(zhuǎn)換為AnnotatedDataFrame對(duì)象
pd <- new("AnnotatedDataFrame",
data=sample_ann)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
# 最后構(gòu)建CDS對(duì)象
sc_cds <- newCellDataSet(
expr_matrix,
phenoData = pd,
featureData =fd,
expressionFamily = negbinomial.size(),
lowerDetectionLimit=1)
關(guān)于其中的參數(shù)expressionFamily ：負(fù)二項(xiàng)分布有兩種方法端朵，這里選用了negbinomial.size 好芭； 另外一種negbinomial稍微更準(zhǔn)確一點(diǎn)，但速度大打折扣冲呢，它主要針對(duì)非常小的數(shù)據(jù)集
monocle V2質(zhì)控過(guò)濾
=> detectGenes() + subset()
cds=sc_cds
cds <- detectGenes(cds, min_expr = 0.1)
# 結(jié)果保存在cds@featureData@data
> print(head(cds@featureData@data))
genes num_cells_expressed
RP11-34P13.7 RP11-34P13.7 0
FO538757.2 FO538757.2 20
AP006222.2 AP006222.2 11
RP4-669L17.10 RP4-669L17.10 0
RP11-206L10.9 RP11-206L10.9 5
LINC00115 LINC00115 0
在monocle版本2.12.0中舍败，取消了fData函數(shù)（此前在2.10版本中還存在）。如果遇到不能使用fData的情況敬拓，就可以采用備選方案：cds@featureData@data
然后進(jìn)行基因過(guò)濾 =>subset()
expressed_genes <- row.names(subset(cds@featureData@data,
num_cells_expressed >= 5))
> length(expressed_genes)
[1] 12273
cds <- cds[expressed_genes,]
monocle V2 聚類
step1：dispersionTable() 目的是判斷使用哪些基因進(jìn)行細(xì)胞分群
當(dāng)然可以使用全部基因邻薯，但這會(huì)摻雜很多表達(dá)量不高而檢測(cè)不出來(lái)的基因，反而會(huì)增加噪音乘凸。最好是挑有差異的厕诡，挑表達(dá)量不太低的
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
disp_table <- dispersionTable(cds) # 挑有差異的
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1) # 挑表達(dá)量不太低的
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id) # 準(zhǔn)備聚類基因名單
plot_ordering_genes(cds)
# 圖中黑色的點(diǎn)就是被標(biāo)記出來(lái)一會(huì)要進(jìn)行聚類的基因
[可省略]step2：plot_pc_variance_explained() 選一下主成分
plot_pc_variance_explained(cds, return_all = F) # norm_method='log'
step3: 聚類
# 進(jìn)行降維
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
reduction_method = 'tSNE', verbose = T)
# 進(jìn)行聚類
cds <- clusterCells(cds, num_clusters = 4)
# Distance cutoff calculated to 1.812595
plot_cell_clusters(cds, 1, 2, color = "cellType")
monocle V2 差異分析
=> differentialGeneTest()
# 這個(gè)過(guò)程比較慢！
start=Sys.time()
diff_test_res <- differentialGeneTest(cds,
fullModelFormulaStr = "~cellType")
end=Sys.time()
end-start
# Time difference of 5.729718 mins
然后得到差異基因
sig_genes <- subset(diff_test_res, qval < 0.1)
> nrow(sig_genes)
[1] 635
# 最后會(huì)得到635個(gè)差異基因
> head(sig_genes[,c("genes", "pval", "qval")] )
genes pval qval
ISG15 ISG15 1.493486e-03 0.040823046
CCNL2 CCNL2 2.228521e-03 0.055590716
MIB2 MIB2 4.954659e-05 0.002523176
MMP23B MMP23B 1.009464e-04 0.004657577
TNFRSF25 TNFRSF25 1.608900e-04 0.006785583
CAMTA1 CAMTA1 4.339489e-03 0.090162380
進(jìn)行熱圖可視化
文章使用的基因如下（也就是第一張圖中顯示的）
htmapGenes=c(
'GAPDH','CD52','TRAC','IL32','ACTB','ACTG1','COTL1',
'GZMA','GZMB','GZMH','GNLY'
)
# 看到作者挑選的這些基因也都出現(xiàn)在我們自己分析的結(jié)果中
> htmapGenes %in% rownames(sig_genes)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
下面開(kāi)始畫圖
先畫一個(gè)最最最原始的
library(pheatmap)
dat=count_matrix[htmapGenes,]
pheatmap(dat)
需要修改的地方：列名（barcode）不需要营勤、聚類不需要灵嫌、數(shù)據(jù)需要z-score
再來(lái)一版
n=t(scale(t(dat)))
# 規(guī)定上下限（和原文保持一致）
n[n>2]=2
n[n< -1]= -1
> n[1:4,1:4]
AAACCTGCAACGATGG.3 AAACCTGGTCTCCATC.3 AAACGGGAGCTCCTTC.3
GAPDH -1.0000000 -0.06057709 -0.37324949
CD52 0.6676167 0.31653833 0.04253204
TRAC 0.7272507 0.63277670 -0.51581312
IL32 -1.0000000 0.42064114 -0.16143114
AAAGCAATCATATCGG.3
GAPDH -1.0000000
CD52 -1.0000000
TRAC -0.5158131
IL32 0.1548693
# 加上細(xì)胞歸屬注釋（ac = annotation column）壹罚，去掉列名和聚類
ac=data.frame(group=cluster)
rownames(ac)=colnames(n)
pheatmap(n,annotation_col = ac,
show_colnames =F,
show_rownames = T,
cluster_cols = F,
cluster_rows = F)
