劉小澤寫于19.7.5-第二單元第七講：每個細胞檢測到的基因數(shù)量差異
筆記目的：根據(jù)生信技能樹的單細胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
前言
原文中根據(jù)5個指標對細胞進行過濾银择，其中第四個是利用有表達量的基因數(shù)量進行過濾
(g) number of exon mapping reads. Cutoff: 10000 (8 cells failed).
(h) Percentage of uniquely mapping reads. Cutoff: 26.11 % (56 cells failed).
(i) percentage of exon mapping reads. Cutoff: 31.15% (40 cells failed).
(j) Number of genes with reads per kilobase gene model and million mappable reads (RPKM)>1. Cutoff: 1112.76 and 7023 (56 cells failed).
(k) Maximum correlation of each cell to all other cells. Cutoff: 0.34 (54 cells failed).
但是要過濾就要有個基礎(chǔ)，也就是有表達量的基因數(shù)量
之前在單細胞轉(zhuǎn)錄組學(xué)習筆記-5：http://www.reibang.com/p/33a7eb26bd31中提到過
# 這里檢測每個樣本中有多少基因是表達的蚕冬，count值以1為標準傲霸，rpkm值可以用0為標準
n_g = apply(a,2,function(x) sum(x>1))
這里主要是重復(fù)文章的一個小提琴圖疆瑰，目的是檢測細胞中可以表達的基因數(shù)量：
先分析一下：橫坐標沒有說明，圖中也沒有分組昙啄，因此原文是將全部的基因都畫在了一起穆役，于是之前構(gòu)建的樣本meta信息中的all這一列就用上了
實際操作
原文使用的是RPKM值
rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')
# 以下是檢查數(shù)據(jù)
dat[1:4,1:4]
> head(metadata)
g plate n_g all
SS2_15_0048_A3 1 0048 3065 all
SS2_15_0048_A6 2 0048 3036 all
SS2_15_0048_A5 1 0048 3742 all
SS2_15_0048_A4 3 0048 5012 all
SS2_15_0048_A1 1 0048 5126 all
SS2_15_0048_A2 3 0048 5692 all
就根據(jù)這個metadata就可以作圖了，步驟就是先鎖定對象(這里的metadata數(shù)據(jù)框)梳凛，然后做映射(y軸用n_g耿币，x軸用all)
先畫第一張圖
library(ggpubr)
ggviolin(metadata, x = "all", y = "n_g", fill = "all",
add = "boxplot", add.params = list(fill = "white"))
然后可以考慮看下批次之間比較
ggviolin(metadata, x = "plate", y = "n_g", fill = "plate",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add = "boxplot", add.params = list(fill = "white"))
另外還有hclust分組間比較
在上圖的基礎(chǔ)上還可以加上p-value，參考STHDA網(wǎng)站韧拒，http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/76-add-p-values-and-significance-levels-to-ggplots/
ggviolin(metadata, x = "g", y = "n_g", fill = "g",
add = "boxplot", add.params = list(fill = "white")) + stat_compare_means()
可以看到差異極顯著
對比一下原始count矩陣得到的hclust分組結(jié)果：
rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input.Rdata')
ggviolin(df, x = "g", y = "n_g", fill = "g",
add = "boxplot", add.params = list(fill = "white")) + stat_compare_means()
當然淹接，還可以實現(xiàn)兩兩比較十性，想比較誰就比較誰：
# 只需要之前構(gòu)建一個比較組合
my_comparisons <- list( c("1", "2"), c("2", "3"), c("3", "4") )
ggviolin(df, x = "g", y = "n_g", fill = "g",
add = "boxplot", add.params = list(fill = "white")) + stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
stat_compare_means(label.y = 50) # Add global p-value
小tip：如果說可視化分群結(jié)果，發(fā)現(xiàn)群組間基因數(shù)量差異太大塑悼，就要考慮技術(shù)差異問題劲适，因為由于生物學(xué)導(dǎo)致幾千個基因關(guān)閉的可能性不是很大，可以換一種聚類算法試一試
目前單細胞也有很多采用dbscan算法進行的聚類分析厢蒜，后續(xù)會介紹
