作者 | 單細(xì)胞天地小編 劉小澤
課程鏈接在：http://jm./index/mulitcourse/detail.html?cid=55
這次會介紹如何針對表達(dá)矩陣進(jìn)行分群,，不一定需要包裝好的函數(shù),。對應(yīng)視頻第三單元5-6講
前言目的是得到這個圖 
將會用到來自作者的包裝好的analysis_functions.R代碼： https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/scripts/analysis_functions.R 這個代碼有1800多行，將會貫穿整個分析,，正是這些DIY的代碼,，才讓文章的圖顯得與眾不同 1 首先創(chuàng)造表達(dá)矩陣先下載作者上游定量處理好的數(shù)據(jù)：female_rpkm.Robj https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/raw/master/data/female_rpkm.Robj 一會要用的基因列表：https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/data/prot_coding.csv
1load(file="female_rpkm.Robj")
2
3## 去掉重復(fù)細(xì)胞
4#（例如：同一個細(xì)胞建庫兩次，這里作者用“rep”進(jìn)行了標(biāo)記）
5grep("rep",colnames(female_rpkm))
6colnames(female_rpkm)[256:257]
7female_rpkm <- female_rpkm[,!colnames(female_rpkm) %in% grep("rep",colnames(female_rpkm), value=TRUE)]
8
9## 只保留編碼基因(去掉類似：X5430419D17Rik,、BC003331等)
10prot_coding_genes <- read.csv(file="prot_coding.csv", row.names=1)
11females <- female_rpkm[rownames(female_rpkm) %in% as.vector(prot_coding_genes$x),]
12save(females,file = 'female_rpkm.Rdata')
2 然后使用包裝好的代碼進(jìn)行tSNE2.1 對細(xì)胞操作=》細(xì)胞發(fā)育時期的獲取細(xì)胞是從6個時間點取出的,，于是先找到這6個時間點 1load('../female_rpkm.Rdata')
2> dim(females)
3[1] 21083 563
4
5> head(colnames(females))
6[1] "E10.5_XX_20140505_C01_150331_1" "E10.5_XX_20140505_C02_150331_1"
7[3] "E10.5_XX_20140505_C03_150331_1" "E10.5_XX_20140505_C04_150331_2"
8[5] "E10.5_XX_20140505_C06_150331_2" "E10.5_XX_20140505_C07_150331_3"
9
10## 取下劃線分隔的第一部分
11female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
12# 或者
13female_stages <- sapply(strsplit(colnames(females), "_"),
14 function(x)x[1])
15# 再或者
16female_stages <- stringr::str_split(colnames(females),'_', simplify = T)[,1]
17
18names(female_stages) <- colnames(females)
19> table(female_stages)
20female_stages
21E10.5 E11.5 E12.5 E13.5 E16.5 P6
22 68 100 103 99 85 108
2.2 對基因操作=》基因過濾與統(tǒng)計去掉在所有細(xì)胞都不表達(dá)的基因1> (dim(females))
2[1] 21083 563
3> females <- females[rowSums(females)>0,]
4> (dim(females))
5[1] 16765 563
可以看到去掉了4000多個 計算各種統(tǒng)計指標(biāo) 1# 利用apply函數(shù)對每行（每個基因）進(jìn)行統(tǒng)計
2mean_per_gene <- apply(females, 1, mean, na.rm = TRUE)
3sd_per_gene <- apply(females, 1, sd, na.rm = TRUE)
4mad_per_gene <- apply(females, 1, mad, na.rm = TRUE)
5cv = sd_per_gene/mean_per_gene
6library(matrixStats)
7var_per_gene <- rowVars(as.matrix(females))
8cv2=var_per_gene/mean_per_gene^2
9# 存儲統(tǒng)計結(jié)果
10cv_per_gene <- data.frame(mean = mean_per_gene,
11 sd = sd_per_gene,
12 mad=mad_per_gene,
13 var=var_per_gene,
14 cv=cv,
15 cv2=cv2)
16rownames(cv_per_gene) <- rownames(females)
17head(cv_per_gene)
18# 根據(jù)表達(dá)量過濾統(tǒng)計結(jié)果
19cv_per_gene=cv_per_gene[cv_per_gene$mean>1,]
20# 簡易的可視化
21with(cv_per_gene,plot(log10(mean),log10(cv2)))
CV值，它表示變異系數(shù)(coefficient of variation),。變異系數(shù)又稱離散系數(shù)或相對偏差 ,，我們肯定都知道標(biāo)準(zhǔn)偏差，也就是sd值,，sd描述了數(shù)據(jù)值偏離算術(shù)平均值的程度,。這個相對偏差CV描述的是標(biāo)準(zhǔn)偏差與平均值之比。 sd值,，它和均值mean,、方差var一樣，都是對一維數(shù)據(jù)進(jìn)行的分析,，如果出現(xiàn)兩組數(shù)據(jù)測量尺度差別太大或數(shù)據(jù)量綱存在差異的話,，直接用標(biāo)準(zhǔn)差就不合適了 CV變異系數(shù)就可以解決這個問題，它利用原始數(shù)據(jù)標(biāo)準(zhǔn)差和原始數(shù)據(jù)平均值的比值來各自消除尺度與量綱的差異,。
復(fù)雜一點的統(tǒng)計可視化：其實就是求每列之間的相關(guān)性 1library(psych)
2pairs.panels(cv_per_gene,
3 method = "pearson", # correlation method
4 hist.col = "#00AFBB",
5 density = TRUE, # show density plots
6 ellipses = TRUE # show correlation ellipses
7)
可以得到不同統(tǒng)計指標(biāo)的關(guān)系 
再用作者包裝的函數(shù)：getMostVarGenes()1females_data <- getMostVarGenes(females, fitThr=2)
2> dim(females_data)
3[1] 822 563
這個函數(shù)也找了822個變化比較大的基因,，用于下游分析，這其實也很像Seurat的FindVariableFeatures()做的事情 
1females_data <- log(females_data+1)
2> females_data[1:4,1:4]
3 E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1
4Ngfr 0 0
5Slc22a18 0 0
6Tspan32 0 0
7Gmpr 0 0
8 E10.5_XX_20140505_C03_150331_1 E10.5_XX_20140505_C04_150331_2
9Ngfr 0.4204863 3.619946
10Slc22a18 0.0000000 0.000000
11Tspan32 0.0000000 0.000000
12Gmpr 0.0000000 0.000000
13
14save(females_data,file = 'females_hvg_matrix.Rdata')
2.3 6個發(fā)育時期RtSNE分析先是PCA針對上面的822個HVGs進(jìn)行操作 1female_sub_pca <- FactoMineR::PCA(
2 t(females_data),
3 ncp = ncol(females_data),
4 graph=FALSE
5)
然后挑選最顯著的主成分,，作為tSNE的輸入 記得在Seurat中是使用ElbowPlot() 關(guān)注肘部的PC,，這里不需要觀察，直接返回最優(yōu)解
1significant_pcs <- jackstraw::permutationPA(
2 female_sub_pca$ind$coord,
3 B = 100,
4 threshold = 0.05,
5 verbose = TRUE,
6 seed = NULL
7)$r
8> significant_pcs
9[1] 9
然后使用上面`jackstraw`挑出的顯著主成分進(jìn)行tSNE 1# 6個時期給定6個顏色
2female_stagePalette <- c(
3 "#2754b5",
4 "#8a00b0",
5 "#d20e0f",
6 "#f77f05",
7 "#f9db21",
8 "#43f14b"
9)
10female_t_sne <- run_plot_tSNE(
11 pca=female_sub_pca,
12 pc=significant_pcs,
13 iter=5000,
14 conditions=female_stages,
15 colours=female_stagePalette
16)
2.4 根據(jù)PCA結(jié)果進(jìn)行層次聚類采用的方法是：Hierarchical Clustering On Principle Components (HCPC) 1# 使用9個顯著主成分重新跑PCA
2res.pca <- FactoMineR::PCA(
3 t(females_data),
4 ncp = significant_pcs,
5 graph=FALSE
6)
7# 作者根據(jù)經(jīng)驗認(rèn)為分成4群比較好解釋,，于是設(shè)置4
8res.hcpc <- FactoMineR::HCPC(
9 res.pca,
10 graph = FALSE,
11 min=4
12)
13# 得到分群結(jié)果
14female_clustering <- res.hcpc$data.clust$clust
15> table(female_clustering)
16female_clustering
17 1 2 3 4
18 90 240 190 43
19# 重新命名
20female_clustering <- paste("C", female_clustering, sep="")
21names(female_clustering) <- rownames(res.hcpc$data.clust)
22# 將C1和C2調(diào)換位置
23female_clustering[female_clustering=="C1"] <- "C11"
24female_clustering[female_clustering=="C2"] <- "C22"
25female_clustering[female_clustering=="C22"] <- "C1"
26female_clustering[female_clustering=="C11"] <- "C2"
27> table(female_clustering)
28female_clustering
29 C1 C2 C3 C4
30240 90 190 43
31
32write.csv(female_clustering, file="female_clustering.csv")
還是基于之前tSNE坐標(biāo),，對聚類得到的4個cluster可視化 1# 為4種cluster設(shè)置顏色
2female_clusterPalette <- c(
3 "#560047",
4 "#a53bad",
5 "#eb6bac",
6 "#ffa8a0"
7)
8
9> head(female_t_sne)
10 tSNE_1 tSNE_2 cond
11E10.5_XX_20140505_C01_150331_1 -2.714291 -24.47912 E10.5
12E10.5_XX_20140505_C02_150331_1 -1.580757 -26.45072 E10.5
13E10.5_XX_20140505_C03_150331_1 -1.577123 -25.36753 E10.5
14E10.5_XX_20140505_C04_150331_2 -6.677577 -20.00208 E10.5
15E10.5_XX_20140505_C06_150331_2 3.442235 -23.32570 E10.5
16E10.5_XX_20140505_C07_150331_3 3.793953 -23.33955 E10.5
17
18# 作者包裝的函數(shù)
19female_t_sne_new_clusters <- plot_tSNE(
20 tsne=female_t_sne,
21 conditions=female_clustering,
22 colours= female_clusterPalette
23)
24ggsave('tSNE_cluster.pdf')
|
