轉(zhuǎn)載來自：劉小澤 鏈接：http://www.itdecent.cn/p/e659a2f164f7
劉小澤寫于19.7.5-第二單元第七講：每個細(xì)胞檢測到的基因數(shù)量差異
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：https://www.bilibili.com/video/BV1dt411Y7nn?p=10
前言
原文中根據(jù)5個指標(biāo)對細(xì)胞進(jìn)行過濾，其中第四個是利用有表達(dá)量的基因數(shù)量進(jìn)行過濾
(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ǔ)，也就是有表達(dá)量的基因數(shù)量
之前在單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-5：http://www.itdecent.cn/p/33a7eb26bd31中提到過
# 這里檢測每個樣本中有多少基因是表達(dá)的，count值以1為標(biāo)準(zhǔn)，rpkm值可以用0為標(biāo)準(zhǔn)
n_g = apply(a,2,function(x) sum(x>1))
這里主要是重復(fù)文章的一個小提琴圖，目的是檢測細(xì)胞中可以表達(dá)的基因數(shù)量：
先分析一下：橫坐標(biāo)沒有說明，圖中也沒有分組，因此原文是將全部的基因都畫在了一起，于是之前構(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()
當(dāng)然，還可以實現(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)閉的可能性不是很大，可以換一種聚類算法試一試
目前單細(xì)胞也有很多采用dbscan算法進(jìn)行的聚類分析，后續(xù)會介紹
作者：劉小澤
鏈接：http://www.itdecent.cn/p/37d4dc8ca391
來源：簡書
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請注明出處。
