细胞聚类分群后检测marker基因
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆发送于19.10.1目标：得到了不同的cluster，就可以利用
findMarkers()去为每个cluster寻找marker基因。来自Bioconductor放出的单细胞数据分析教程
https://bioconductor.org/packages/release/workflows/vignettes/simpleSingleCell/inst/doc/reads.html
它的方法是：两两cluster之间对每个基因的log表达量进行 Welch t-tests (Soneson and Robinson 2018)。这样就能找到每个cluster中的差异表达基因（前提是掩盖掉非生物因素，例如批次信息）。得到的top 差异基因能更好地作为marker基因候选，原因就是：既然它们是最具有差异代表性的，那么利用它们就能区分不同cluster的细胞，也就实现了marker的目的。【当然，还有许多细胞类型自己特有的基因，也是作为marker的候选者】
markers <- findMarkers(sce, my.clusters, block=sce$Plate)
每个cluster与其他cluster的两两比较结果会整合到每个cluster的差异基因结果数据框中
例如：看一下第1个cluster的top10差异基因
marker.set <- markers[["1"]]
head(marker.set, 10)
## DataFrame with 10 rows and 7 columns
## Top p.value FDR logFC.2
## <integer> <numeric> <numeric> <numeric>
## Aurkb 1 6.65863261824492e-75 1.58362259559721e-70 -7.37163350569914
## Tk1 1 6.41442387689833e-64 3.81385607660686e-60 -4.92754579848056
## Myh11 1 4.95086465211539e-49 9.81220116843835e-46 4.42159318376489
## Cdca8 1 2.22334940769686e-46 3.5251945975503e-43 -6.84273526334783
## Ccna2 2 1.16841222330534e-68 1.38941739534356e-64 -7.3079756458424
## Rrm2 2 1.48873842020081e-56 5.05809512109088e-53 -5.52120322191947
## Cks1b 2 3.83636139889977e-39 2.40105745131667e-36 -6.6792118199827
## Pirb 2 1.83893803490065e-34 6.15992440620318e-32 5.25803749166673
## Pimreg 3 7.41004737119548e-68 5.87443855430467e-64 -7.30454210816126
## Pclaf 3 8.9651722100101e-51 2.13218690670669e-47 -5.60087985621813
## logFC.3 logFC.4 logFC.5
## <numeric> <numeric> <numeric>
## Aurkb -6.72799345321135 -1.95039440976238 -6.91128802096519
## Tk1 -7.74065113926215 -3.53749565362853 -4.63516273649457
## Myh11 4.30812918035861 4.45235717737968 1.0413149433198
## Cdca8 -4.88595092732349 -2.43821402084038 -7.12791471326961
## Ccna2 -6.9676852052539 -2.46589325823089 -7.12692721843331
## Rrm2 -7.94685699818751 -3.19173143688883 -5.42878091964762
## Cks1b -5.92137181826573 -4.37146346518812 -6.214592473138
## Pirb 5.18195596569259 5.87631057587633 0.0704964218555272
## Pimreg -5.91099762335684 -0.874660676141792 -7.01798853404623
## Pclaf -7.56997893312346 -2.36631043435985 -5.16956927698937
可以将全部cluster的marker基因保存起来
write.table(marker.set, file="416B_marker_1.tsv", sep="\t",
quote=FALSE, col.names=NA)
最后热图对cluster1的marker基因可视化
top.markers <- rownames(marker.set)[marker.set$Top <= 10]
plotHeatmap(sce, features=top.markers, columns=order(sce$cluster),
colour_columns_by=c("cluster", "Plate", "Oncogene"),
cluster_cols=FALSE, center=TRUE, symmetric=TRUE,
zlim=c(-5, 5),
show_colnames = F)
从图中粗略地看一看：cluster1包含了处理组（oncogene-induced）细胞，另外第5群也是处理组细胞。如果生物背景足够的话，也能推断出这些细胞与某些DE基因之间的联系。
「需要注意的地方」
当然，不是所有的marker基因都在所选cluster都是上调或下调的，它不仅仅是差异基因这么简单。如果只选差异基因，这样会过滤掉很多重要的marker基因。
例如，在CD4+-only, CD8+-only, double-positive and double-negative T cells混合的细胞群体中，重要的细胞marker基因Cd4或Cd8 都不会被检测到，因为这两个基因会在其中两个类型的细胞群体中都存在（比如Cd4 会存在于CD4+-only和 double-positive两种细胞类型中）
这里强烈建议挑选一些已经验证过的细胞marker基因，比如
CellMarker数据库（http://biocc.hrbmu.edu.cn/CellMarker/）：整理了100,000+发表的文献，包含人的158个组织 (亚组织)的467个细胞类型的13,605个Marker基因，和鼠的81个组织 (亚组织)的389个细胞类型的9, 148个Marker基因
PanglaoDB（https://panglaodb.se/）：小鼠的170种组织954个样本近400W细胞和来自人的68种组织279个样本100w+细胞
小鼠Mouse Cell Atlas（http://bis.zju.edu.cn/MCA/）
signatureDB（https://pldaniels.com/signaturedb/）
最好再加上一些定量实验，比如fluorescent in situ hybridisation or quantitative PCR，确定细胞亚群真实存在，而不仅仅是因为我们分析得到是这样
作者建议
如果只想提取cluster中的上调表达基因，就可以在
findMarkers中指定参数direction="up"。应用在高度异质性群体中，可以更快速判断不同的clusterfindMarkers除了将一个cluster与其他各个cluster两两比较之外，还能对某个cluster和其他全部clusters进行差异分析，设置参数pval.type="all"；如果再和direction="up"连用，就可以用来鉴定每个cluster中特异的marker基因。但是如果存在”过度聚类“的情况，原本一个cluster分成了两个，那么这个特异的marker基因就不会存在。另外，这里的p值不能简单理解为显著性。具体原因在:
https://bioconductor.org/packages/3.9/simpleSingleCell/vignettes/de.html#misinterpretation-of-de-$p$-values
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步。
我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到jieandze1314@gmail.com，每一条都会看到的哦~