這次跟著課程（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-1-前言及上游介紹）要練習(xí)的文章是：Dissecting Cell Lineage Specification and Sex Fate Determination in Gonadal Somatic Cells Using Single-Cell Transcriptomics不撑。
課程里是從下載sra文件開始的照筑，但是由于這篇文章的數(shù)據(jù)實(shí)在是太大了肠阱。立倍。刑棵。
都下載下來的話毁兆，按照我這網(wǎng)速屑迂，可能要一年以后才能下載完了瑞驱。所以我直接從課程的第二講開始練習(xí)（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-2-根據(jù)表達(dá)矩陣進(jìn)行分群）。關(guān)于單細(xì)胞測(cè)序的數(shù)據(jù)如何批量下載和屎，如何比對(duì)拴驮，還有featureCount定量，我在上一個(gè)練習(xí)里有寫到（單細(xì)胞測(cè)序?qū)崙?zhàn)（第一部分））柴信。這次主要練習(xí)下游的分析部分套啤，這部分是最重要的，也是核心部分随常。還有個(gè)原因就像課程說的（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-1-前言及上游介紹）：
smartseq2得到的兩個(gè)R1潜沦、R2都是測(cè)序信息，于是它的操作和我們常規(guī)bulk轉(zhuǎn)錄組是類似的绪氛∷艏Γ可以用hisat2+featureCounts進(jìn)行操作
需要明確的是：文章表明了6個(gè)發(fā)育時(shí)期、4群細(xì)胞枣察、2個(gè)發(fā)育軌跡這三種細(xì)胞屬性
首先先下載原始count矩陣https://raw.githubusercontent.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/master/data/female_count.Robj
作者rpkm標(biāo)準(zhǔn)化后的矩陣https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/raw/master/data/female_rpkm.Robj
（一）過濾重復(fù)的細(xì)胞争占，構(gòu)建表達(dá)矩陣（使用原始count矩陣)
#導(dǎo)入作者上傳的原始矩陣
> load(file="female_count.Robj")
> load(file="female_rpkm.Rdata")
# 去掉重復(fù)細(xì)胞（作者用“rep”進(jìn)行了標(biāo)記），先看一下哪些列是重復(fù)的
> grep("rep",colnames(female_count))
[1] 257 #第257列是重復(fù)的
# 直接對(duì)細(xì)胞和基因過濾
female_count <- female_count[rownames(female_count) %in% rownames(females),!colnames(female_count) %in% grep("rep",colnames(female_count), value=TRUE)]
# %in%相當(dāng)于match()函數(shù)的一個(gè)縮寫序目。用來判斷一個(gè)數(shù)組或矩陣（前面的）是否包含在另一個(gè)數(shù)組或矩陣（后面的）里臂痕。
#看一下過濾之后的
> female_count[1:3,1:3]
E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1 E10.5_XX_20140505_C03_150331_1
eGFP 19582 526 4786
Gnai3 2218 122 4
Pbsn 0 0 0
#保存一下
> save(female_count,file = 'female_count.Rdata')
（二）tSNE分析（Seurat方法）
在分析之前，先簡(jiǎn)單的了解一下什么是Seurat猿涨，這是這個(gè)包的官方網(wǎng)站（SATIJA LAB）握童，據(jù)說這個(gè)包今年剛升過級(jí)，這個(gè)包還可以很友好的直接讀取10x Genomics的輸出結(jié)果（單細(xì)胞轉(zhuǎn)錄組測(cè)序分析--初探Seurat）叛赚。有幾個(gè)比較好的文章對(duì)于這個(gè)包進(jìn)行介紹：
1.Seurat使用教程（v3.0）
2.單細(xì)胞Seurat包升級(jí)澡绩，換湯不換藥
3.單細(xì)胞轉(zhuǎn)錄組測(cè)序分析--初探Seurat
下面進(jìn)行分析，首先要對(duì)細(xì)胞分群俺附，按照不同的發(fā)育時(shí)期來獲取細(xì)胞群：
> load(file="female_count.Rdata")
> dim(females) #看一下多少行英古，多少列
[1] 21083 563
> head(colnames(females)) #看一下列名
[1] "E10.5_XX_20140505_C01_150331_1" "E10.5_XX_20140505_C02_150331_1" "E10.5_XX_20140505_C03_150331_1"
[4] "E10.5_XX_20140505_C04_150331_2" "E10.5_XX_20140505_C06_150331_2" "E10.5_XX_20140505_C07_150331_3"
#按照胚胎的天數(shù)來分細(xì)胞群，所以要按照E10.5這部分來劃分
> female_stages <- stringr::str_split(colnames(females),'_', simplify = T)[,1]
#str_split意思是按照模式分割字符串昙读，這里是按照行名里的下劃線分割召调，后面的1是第一個(gè)下劃線的意思
> head(female_stages)
[1] "E10.5" "E10.5" "E10.5" "E10.5" "E10.5" "E10.5"
#表達(dá)矩陣?yán)镉卸嗌倭校ǘ嗌偌?xì)胞），female_stages里就有多少元素蛮浑，把列名賦值給它
> names(female_stages) <- colnames(females)
> table(female_stages)
female_stages
E10.5 E11.5 E12.5 E13.5 E16.5 P6
68 100 103 99 85 108
現(xiàn)在就得到了6個(gè)不同的群，每個(gè)群的細(xì)胞數(shù)不同
下面使用seurat進(jìn)行tSNE分析沮稚，分幾步進(jìn)行：
（1）構(gòu)建對(duì)象
#安裝seurat包
> install.packages('Seurat')
> library(Seurat)
#構(gòu)建對(duì)象
> sce_female <- CreateSeuratObject(counts = female_count,
project = "sce_female",
min.cells = 1, min.features = 0)
# 這一步還有過濾功能艺沼，基因至少在min.cells個(gè)細(xì)胞中表達(dá)，每個(gè)細(xì)胞中至少表達(dá)min.genes個(gè)基因/features蕴掏。這句的代碼是這個(gè)基因至少要在1個(gè)細(xì)胞里有表達(dá)障般，不然就被過濾掉调鲸。
> sce_female
An object of class Seurat
16765 features across 563 samples within 1 assay
Active assay: RNA (16765 features)
view一下這個(gè)sce_female:
點(diǎn)開里面的meta.data，可以看到UMI的情況:
我搜了一下網(wǎng)上挽荡，想知道這些都是些啥藐石，在這篇文章里得知（單細(xì)胞轉(zhuǎn)錄組測(cè)序分析--初探Seurat）：創(chuàng)建完Seurat對(duì)象后，Seurat將數(shù)據(jù)保存在不同的slot中定拟，如sce_female@raw.data, sce_female@data, sce_female@meta.data, sce_female@ident于微，其中meta.data存放的是每個(gè)細(xì)胞的統(tǒng)計(jì)數(shù)據(jù)如UMI數(shù)目，gene數(shù)目等青自，ident此時(shí)存放的是project信息株依。
（2）添加樣品注釋
#使用AddMetaData添加meta信息
> sce_female <- AddMetaData(object = sce_female,
metadata = apply(female_count, 2, sum),
col.name = 'nUMI_raw')
> sce_female <- AddMetaData(object = sce_female,
metadata = female_stages,
col.name = 'female_stages') #把上面我們分的細(xì)胞的不同時(shí)期加到這個(gè)對(duì)象里
（3）標(biāo)準(zhǔn)化
Seurat官網(wǎng)上對(duì)于標(biāo)準(zhǔn)化的這部分是這樣寫的：
After removing unwanted cells from the dataset, the next step is to normalize the data. By default, we employ a global-scaling normalization method “LogNormalize” that normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. Normalized values are stored in @data.
去除掉不想要的細(xì)胞以后，下一步就是標(biāo)準(zhǔn)化延窜。默認(rèn)值下恋腕，我們使用global-scaling標(biāo)準(zhǔn)化方法，稱為“LogNormalize”逆瑞。對(duì)每個(gè)基因進(jìn)行標(biāo)準(zhǔn)化：該基因的UMI除以該細(xì)胞內(nèi)轉(zhuǎn)錄物的總UMI并乘以10000（默認(rèn)值）荠藤，再取log值作為結(jié)果。
> sce_female <- NormalizeData(sce_female) #這是簡(jiǎn)化版的代碼呆万，也是官方推薦使用的
#這句代碼的完整格式應(yīng)該是這樣：
#sce_female <- NormalizeData(object =sce_female商源，normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> sce_female[["RNA"]]@data[1:3,1:3] #看一眼標(biāo)準(zhǔn)化之后的
3 x 3 sparse Matrix of class "dgCMatrix"
E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1 E10.5_XX_20140505_C03_150331_1
eGFP 3.753731 0.7024916 2.8485868
Gnai3 1.744141 0.2121183 0.0135009
Cdc45 2.001816 1.6581178 0.6510173
（4）找差異基因Identification of highly variable features (feature selection)
> sce_female <- FindVariableFeatures(sce_female,
selection.method = "vst",
nfeatures = 2000)
# nfeatures = 2000是默認(rèn)值
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
# HVGs可視化
> VariableFeaturePlot(sce_female)
如果想看上面這個(gè)圖里偏上的紅點(diǎn)都是什么基因车份，還可以再用兩句代碼：
> top10 <- head(VariableFeatures(sce_female), 10)
> plot1<- VariableFeaturePlot(sce_female)
> plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
> plot2
（5）數(shù)據(jù)比例化（調(diào)整數(shù)據(jù)范圍）
在降維之前還可以有一步就是數(shù)據(jù)的比例化谋减，官方是這樣解釋為什么要進(jìn)行這一步的：
(1)Shifts the expression of each gene, so that the mean expression across cells is 0
(2)Scales the expression of each gene, so that the variance across cells is 1
(3)This step gives equal weight in downstream analyses, so that highly-expressed genes do not dominate
The results of this are stored in pbmc[["RNA"]]@scale.data
說的簡(jiǎn)單一點(diǎn)就是，有些基因的表達(dá)量非常高扫沼，會(huì)掩蓋掉低表達(dá)基因的變化出爹，這一步就是降低高表達(dá)基因的影響。
#默認(rèn)只對(duì)FindVariableFeatures得到的HVGs進(jìn)行操作
> sce_female <- ScaleData(object = sce_female,
vars.to.regress = c('nUMI_raw'),
model.use = 'linear',
use.umi = FALSE)
#vars.to.regress的意思就是需要去除的變異來源缎除。該函數(shù)會(huì)先進(jìn)行先行回歸分析严就，然后再進(jìn)行數(shù)據(jù)的比例化
Regressing out nUMI_raw
|=============================================================================================================| 100%
Centering and scaling data matrix
|=============================================================================================================| 100%
（6）PCA降維
對(duì)于降維：PCA和tSNE都是數(shù)據(jù)降維分析方法，PCA屬于線性降維器罐，tSNE屬于非線性降維（單細(xì)胞轉(zhuǎn)錄組測(cè)序分析--初探Seurat）梢为。
為什么要降維？單細(xì)胞測(cè)序數(shù)據(jù)是一個(gè)高維的復(fù)雜數(shù)據(jù)轰坊。降維铸董，就是把高維數(shù)據(jù)通過優(yōu)化保留原始數(shù)據(jù)中的關(guān)鍵特征后投射到低維空間，從而可以通過二維或三維的形式把數(shù)據(jù)展示出來肴沫。先執(zhí)行PCA分析粟害，使高維數(shù)據(jù)的信息最大程度保留在低維數(shù)據(jù)中。
> sce_female <- RunPCA(sce_female,
features = VariableFeatures(object = sce_female))
#然后會(huì)彈出下面一堆：
PC_ 1
Positive: Kctd14, Cgn, Itga6, Kitl, Gstm2, Bhlhe41, Lzts1, eGFP, Serpinb6a, Ttyh2
Soat1, Ifitm1, Rnf213, Gstm1, Idh1, Pank1, Anxa2, Gpr56, Klhl23, Sdc4
Hmgcs2, Laptm4b, Tmprss13, Tmem206, Clu, Rgs11, Serpinb6b, Gadd45g, Cpa2, Gatm
Negative: Col1a2, Col1a1, Sfrp1, Sept11, Ncam1, Ptn, Mmp2, Fn1, Col3a1, Hmga2
Nrep, Sparcl1, Tcf21, Mfap4, Sox11, Meis1, Cxcl12, Nfib, Bgn, Cnn2
Pdgfra, Efna5, Dcn, Cdh11, Fstl1, Fblim1, Vcan, Top2a, Tuba1a, Lpar1
PC_ 2
Positive: Itih5, Cfh, Dcn, Col6a2, Igf1, Akr1cl, Ogn, Col1a1, Lama2, Sparcl1
Lum, Zbtb20, Cd34, Col6a1, Ptrf, Rarres2, Col1a2, Bgn, Gas6, Lrrc17
Cped1, Col4a4, Srpx2, Col6a3, Prss35, Plxnc1, Serpine2, Slc26a7, Tcf21, Pdlim3
Negative: Hmga2, Peg3, Mest, Sox11, Cdkn1c, Rrm2, Rgs5, Ccnb1, Aldh1a2, Ccnd1
Oraov1, Asb4, Tpx2, Cenpf, Flrt3, Hapln1, Kif20a, Ect2, Top2a, Meis1
Cenpe, F11r, Kif22, Bub1, Sfrp2, Kif2c, Krt18, Ube2c, Ncaph, Sulf1
PC_ 3
Positive: Cdkn1c, Ifitm1, Peg3, Fstl1, Acaa2, Nfib, Lzts1, Cgn, Anxa6, Tnfrsf19
Lhx9, Laptm4b, Col3a1, Sema6d, Col1a1, Sulf1, Zfp57, Serpinb6b, Bgn, Timp3
eGFP, Anxa2, Lgr4, Col1a2, Colec12, Itm2a, Hmcn1, Mest, Serpini1, Bhlhe41
Negative: Serpine2, Socs2, Inha, Kazald1, Lect1, Ephx2, Ptges, Hsd3b1, Hsd17b1, Apoa1
Myo1e, Ubash3b, Esr2, Fam13a, Arhgap31, Gls, Gja1, Slc18a2, Ablim1, Tanc2
Slc26a7, Map1b, Sytl4, Myo5c, Tenm4, Ivns1abp, Thbs1, Epn2, Top2a, Akr1c14
PC_ 4
Positive: Top2a, Prc1, Ube2c, Ccnb1, Cenpf, Tpx2, Plk1, Kif22, Ckap2l, Bub1
Kif23, Mki67, Aspm, Cdk1, eGFP, Bub1b, Ccna2, Kif2c, Ncaph, Ect2
Cenpe, Cdc25c, Ccnb2, Kif20a, Ifitm1, Nr6a1, Zwilch, Sgol1, Rrm2, Gm3550
Negative: Alcam, Ptgs1, Fras1, Gata2, Upk1b, Ccnd2, Dpp4, Runx1t1, Pkhd1l1, Frem2
Vstm2a, Cpa6, Efna5, Nfib, Itm2a, Flrt1, Mn1, Sdc4, Plscr2, Gpm6a
Sema3c, Egr1, Mybpc1, Btg2, Krt18, Atp2b1, Zfp36, Sorbs2, Timp2, Ece1
PC_ 5
Positive: Nup62cl, Cyp11a1, Pnlip, Gadd45g, Sphkap, Asb4, Pnliprp2, Akr1cl, Inha, Plac1
Gfra3, Asns, Siglecg, Foxp2, Pls3, Nefm, St3gal5, Adam8, Mc2r, Vav3
Gk5, Cpa1, Sprr2d, Tenm4, Hmga2, Cxcl12, Gbp5, Trpc7, Ttc38, Inhbb
Negative: Sorbs2, Cav1, Btg2, Pdpn, Lgr5, Prc1, Egr1, Cst12, Mt1, Ccnb1
Ube2c, Fos, Zfp36, Aspm, Sdc4, Cenpf, Ckap2l, Top2a, Tpx2, Mki67
Egr2, Cdk1, Kif23, Plk1, Tnfrsf19, Kif20a, Pde8b, Bub1, Kif22, Enpp2
Seurat包提供了多種PCA可視化的方法颤芬，感興趣的同學(xué)可以參考：https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html官網(wǎng)學(xué)習(xí)悲幅。這里只舉例一種：
> DimPlot(sce_female, reduction = "pca")
還可以畫個(gè)熱圖套鹅，舉個(gè)例子，畫PC_1和PC_2兩個(gè)dim的熱圖：
> DimHeatmap(object = sce_female, dims = 1:2, cells = 500, balanced = TRUE)
（7）降維以后聚類
在進(jìn)行聚類分析之前汰具，需要選擇使用對(duì)少個(gè)主成分進(jìn)行計(jì)算卓鹿。每個(gè)主成分實(shí)際上代表的是相關(guān)基因集的信息，因此確定多少個(gè)主成分是一個(gè)重要的步驟郁副，我們可以根據(jù)PCElbowPlot函數(shù)來判斷减牺。
#Seurat采用的是基于圖形的聚類方法，即利用PCA空間中的歐幾里德距離構(gòu)造一個(gè)KNN圖存谎。官網(wǎng)上寫了一堆數(shù)學(xué)概念拔疚，我也看不懂。既荚。稚失。
> sce_female <- FindNeighbors(sce_female, dims = 1:20)
Computing nearest neighbor graph
Computing SNN
> sce_female <- FindClusters(sce_female, resolution = 0.3)
#官方解釋：將resolution參數(shù)設(shè)置在0.4-1.2之間，對(duì)于3K左右的單細(xì)胞數(shù)據(jù)集通常會(huì)得到良好的結(jié)果恰聘。對(duì)于較大的數(shù)據(jù)集句各，最佳分辨率通常會(huì)增加。課程里這里用了0.3晴叨，為了盡量還原課程里的圖凿宾，我也用了0.3。（resolution會(huì)影響聚類的個(gè)數(shù)）
(8)tSNE
#在進(jìn)行tSNE之前要先確定有多少個(gè)主成分兼蕊，可以用PCElbowPlot函數(shù)來看
> ElbowPlot(sce_female)
可以看到初厚，拐點(diǎn)出現(xiàn)在9左右，所以我們可以選擇9來進(jìn)行聚類分析孙技。
（9）tSNE可視化
Seurat提供了幾種非線性的降維技術(shù)产禾，如tSNE和UMAP。
> sce_female_tsne <- RunTSNE(sce_female, dims = 1:9)
# 6個(gè)發(fā)育時(shí)間
> DimPlot(object = sce_female_tsne, reduction = "tsne",
group.by = 'female_stages')
# 4個(gè)cluster
> DimPlot(sce_female_tsne, reduction = "tsne")
如果用UMAP來降維：
> sce_female_umap <- RunUMAP(sce_female, dims = 1:9)
> DimPlot(sce_female_umap, reduction = "umap",group.by = 'female_stages')
> DimPlot(sce_female_umap, reduction = "umap")
（三）tSNE分析（DBSCAN牵啦、kmeans方法）
上面用了Seurat包做了分群亚情，課程里還用了DBSCAN方法跑了一遍，所以我也練習(xí)了一下哈雏。
# 這個(gè)運(yùn)行會(huì)非常慢楞件！
> if(T){
library(Rtsne)
N_tsne <- 10
# 其實(shí)這里應(yīng)該多運(yùn)行一些，例如N_tsne <- 50裳瘪，但是會(huì)非常非常慢土浸。。盹愚。
tsne_out <- list(length = N_tsne)
KL <- vector(length = N_tsne)
set.seed(1234)
for(k in 1:N_tsne)
{
tsne_out[[k]]<-Rtsne(t(log2(females+1)),initial_dims=30,verbose=FALSE,check_duplicates=FALSE,
perplexity=27, dims=2,max_iter=5000)
KL[k]<-tail(tsne_out[[k]]$itercosts,1)
print(paste0("FINISHED ",k," TSNE ITERATION"))
}
names(KL) <- c(1:N_tsne)
opt_tsne <- tsne_out[[as.numeric(names(KL)[KL==min(KL)])]]$Y
#在Y里儲(chǔ)存了畫圖坐標(biāo)
}
這個(gè)代碼運(yùn)行的時(shí)候栅迄，會(huì)彈出來一堆一下的東西：
[1] "FINISHED 1 TSNE ITERATION"
[1] "FINISHED 2 TSNE ITERATION"
[1] "FINISHED 3 TSNE ITERATION"
[1] "FINISHED 4 TSNE ITERATION"
[1] "FINISHED 5 TSNE ITERATION"
[1] "FINISHED 6 TSNE ITERATION"
[1] "FINISHED 7 TSNE ITERATION"
[1] "FINISHED 8 TSNE ITERATION"
[1] "FINISHED 9 TSNE ITERATION"
[1] "FINISHED 10 TSNE ITERATION"
# DBSCAN結(jié)果
> library(dbscan)
> plot(opt_tsne, col=dbscan(opt_tsne,eps=3.1)$cluster,
pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
# kmeans方法
plot(opt_tsne, col=kmeans(opt_tsne,centers = 4)$clust,
pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
（四）標(biāo)記基因可視化
這一部分將跟著課程（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-3-標(biāo)記基因可視化）進(jìn)行練習(xí)。我們已經(jīng)把細(xì)胞分了群皆怕，接下來就要分析一下這些群間的差異性毅舆。Seurat包也可以做這個(gè)事情西篓，F(xiàn)indAllMarkers能夠同時(shí)計(jì)算所有亞群的差異性（分別計(jì)算每個(gè)亞群與剩下的所有細(xì)胞的差異性）。但課程里寫到原文作者根據(jù)胚胎細(xì)胞發(fā)育的不同時(shí)期憋活，選取了一些能代表不同時(shí)期的特異性marker進(jìn)行后續(xù)分析岂津。所以這里我還是跟隨教程走一遍。
目標(biāo)：根據(jù)之前分析的6個(gè)發(fā)育時(shí)期和4個(gè)cluster的tSNE結(jié)果悦即，畫一些marker基因的熱圖吮成、小提琴圖。
首先要明確作者選了哪些基因（這些基因在原文中作者有提到是分別代表了哪些細(xì)胞群的marker基因）：
# 作者選擇的14個(gè)marker基因
> markerGenes <- c(
"Nr2f1",
"Nr2f2",
"Maf",
"Foxl2",
"Rspo1",
"Lgr5",
"Bmp2",
"Runx1",
"Amhr2",
"Kitl",
"Fst",
"Esr2",
"Amh",
"Ptges"
)
小提琴圖辜梳，可以直接查看選取的這些marker基因在不同的cluster里的表達(dá)情況：
> VlnPlot(object = sce_female_tsne, features = markerGenes ,
pt.size = 0.2,ncol = 4)
熱圖：
> FeaturePlot(object = sce_female_tsne, features = markerGenes ,
pt.size = 0.2,ncol = 3)
如果是想單獨(dú)畫一個(gè)箱線圖（比如只畫Nrf2f這個(gè)基因）：
# 分類信息
> group <- Seurat::Idents(sce_female)
## 表達(dá)矩陣
> nr2f2 <- as.numeric(log(female_count['Nr2f2',]+1))
#出圖
> boxplot(nr2f2~group)
覺得丑？那就畫個(gè)看起來高大上的小提琴圖吧作瞄。用ggviolin畫
#nrf2f的小提琴圖
> df <- data.frame(expr=nr2f2,
group=group) #首先整一個(gè)data.frame茶宵，以為ggplot要用data.frame才行
看一下這個(gè)df長(zhǎng)什么樣
> view(df)
> female_clusterPalette <- c(
"#560047",
"#a53bad",
"#eb6bac",
"#ffa8a0"
) #這一步是給小提琴圖上顏色
> my_comparisons <- list( c("0", "1"), c("1", "2"), c("2", "3") ) #還想給這四組兩兩比較
> ggviolin(df, "group", "expr", fill = "group",
palette = female_clusterPalette,
add = "boxplot", add.params = list(fill = "white"))+
stat_compare_means(comparisons = my_comparisons)
#這里"group"和"expr"分別代表了x軸和y軸
（五）monocle V2進(jìn)行差異分析
這一部分按照教程（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-4-差異分析及功能注釋）來練習(xí)∽诨樱現(xiàn)在就要對(duì)之前分出來的4群細(xì)胞進(jìn)行差異分析乌庶。這部分主要是利用原文文獻(xiàn)作者上傳的代碼來進(jìn)行分析。另外契耿，ROTS包也可以用來進(jìn)行單細(xì)胞測(cè)序的差異分析瞒大，我這里沒有練習(xí)，之后有時(shí)間可以再練習(xí)一下搪桂。具體ROTS包流程代碼可見：Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-4-差異分析及功能注釋
（1）準(zhǔn)備表達(dá)矩陣和分群信息
#獲取6個(gè)發(fā)育時(shí)期
> head(colnames(female_count))
[1] "E10.5_XX_20140505_C01_150331_1" "E10.5_XX_20140505_C02_150331_1" "E10.5_XX_20140505_C03_150331_1" "E10.5_XX_20140505_C04_150331_2"
[5] "E10.5_XX_20140505_C06_150331_2" "E10.5_XX_20140505_C07_150331_3"
> female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
> names(female_stages) <- colnames(female_count)
> table(female_stages)
female_stages
E10.5 E11.5 E12.5 E13.5 E16.5 P6
68 100 103 99 85 108
#這部分跟前面是一樣的
接下來我用的是下載下來的clustering文件（https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/data/female_clustering.csv）：
> cluster <- read.csv('female_clustering.csv')
> female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
> table(female_clustering)
female_clustering
C1 C2 C3 C4
240 90 190 43
（2）構(gòu)建對(duì)象
需要：表達(dá)矩陣透敌、細(xì)胞信息、基因信息
# 表達(dá)矩陣
> dim(female_count)
# 細(xì)胞分群信息（包括6個(gè)stage和4個(gè)cluster）
> table(female_stages)
female_stages
E10.5 E11.5 E12.5 E13.5 E16.5 P6
68 100 103 99 85 108
> table(female_clustering)
female_clustering
C1 C2 C3 C4
240 90 190 43
# 基因信息
> gene_annotation <- as.data.frame(rownames(female_count))
#開始構(gòu)建
> library(monocle)
這里教程使用了文獻(xiàn)中的源代碼prepare_for_DE锅棕，作者把下面這些代碼包裝成一個(gè)function:
> prepare_for_DE <- function(count_matrix=count_matrix, clustering=clustering, stages=stages){
# Prepare tables for monocle object
expr_matrix <- as.matrix(count_matrix)
sample_sheet <- data.frame(cells=names(count_matrix), stages=stages, cellType=clustering)
rownames(sample_sheet)<- names(count_matrix)
gene_annotation <- as.data.frame(rownames(count_matrix))
rownames(gene_annotation)<- rownames(count_matrix)
colnames(gene_annotation)<- "genes"
pd <- new("AnnotatedDataFrame", data = sample_sheet)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
# Create a CellDataSet from the relative expression levels
HSMM <- newCellDataSet(
as(expr_matrix, "sparseMatrix"),
phenoData = pd,
featureData = fd,
lowerDetectionLimit=0.5,
expressionFamily=negbinomial.size()
)
HSMM <- detectGenes(HSMM, min_expr = 5)
# HSMM <- HSMM[fData(HSMM)$num_cells_expressed > 5, ]
HSMM <- HSMM[fData(HSMM)$num_cells_expressed > 10, ]
HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)
return(HSMM)
}
然后直接使用這個(gè)function:
#包裝好的函數(shù)只需要提供表達(dá)矩陣和細(xì)胞信息
> DE_female <- prepare_for_DE (
female_count,
female_clustering,
female_stages
)
下面這一步是篩選那些差異較大的基因（FDR<0.05拙泽，也就是q值小于0.05）淌山，這里教程還是使用了文獻(xiàn)里的源代碼進(jìn)行篩選裸燎。源代碼包成了一個(gè)包：
#第一次運(yùn)行的話，需要自己運(yùn)行一遍這個(gè)泼疑，這樣自己的R里就會(huì)儲(chǔ)存這個(gè)function德绿，之后就可以直接調(diào)用了
> findDEgenes <- function(HSMM=HSMM, qvalue=qvalue){
diff_test_res <- differentialGeneTest(
HSMM,
fullModelFormulaStr="~cellType",
cores = 3
)
sig_genes_0.05 <- subset(diff_test_res, qval < 0.05)
sig_genes_0.01 <- subset(diff_test_res, qval < 0.01)
print(paste(nrow(sig_genes_0.05), " significantly DE genes (FDR<0.05).", sep=""))
print(paste(nrow(sig_genes_0.01), " significantly DE genes (FDR<0.01).", sep=""))
diff_test_res <- subset(diff_test_res, qval< qvalue)
return(diff_test_res)
}
然后直接運(yùn)行這個(gè)function:
> start_time <- Sys.time()
> female_DE_genes <- findDEgenes(
DE_female,
qvalue=0.05
)
[1] "4435 significantly DE genes (FDR<0.05)." #也就是說最后我們篩選出來了4435個(gè)差異顯著的基因
[1] "3268 significantly DE genes (FDR<0.01)."
> end_time <- Sys.time()
> end_time - start_time
Time difference of 1.661446 mins
得到了顯著的差異基因，需要知道這些基因都是屬于哪一個(gè)cluster的退渗。這里文獻(xiàn)的作者仍然包了一個(gè)函數(shù)get_up_reg_clusters()移稳，這個(gè)函數(shù)的作用是：先得到了每個(gè)差異基因在不同cluster的平均表達(dá)量，然后找平均表達(dá)量最大的那個(gè)cluster会油，就認(rèn)為這個(gè)基因?qū)儆谶@個(gè)cluster个粱。get_up_reg_clusters源代碼是：
> get_up_reg_clusters <- function(count, clustering, DE_genes){
cluster_nb <- unique(clustering)
mean_per_cluster <- vector()
DE_genes <- DE_genes[order(rownames(DE_genes)),]
count <- count[order(rownames(count)),]
count_de_genes <- count[rownames(count) %in% DE_genes$genes,]
print(dim(count_de_genes))
for (clusters in cluster_nb) {
# print(head(count_de_genes[,
# colnames(count_de_genes) %in% names(clustering[clustering==clusters])
# ]))
mean <- rowMeans(
as.matrix(count_de_genes[,
colnames(count_de_genes) %in% names(clustering[clustering==clusters])
])
)
names(mean) <- clusters
mean_per_cluster <- cbind(
mean_per_cluster,
mean
)
}
colnames(mean_per_cluster) <- cluster_nb
up_reg_cluster <- colnames(mean_per_cluster)[apply(mean_per_cluster,1,which.max)]
de_genes_table <- data.frame(
DE_genes,
mean_per_cluster,
cluster=up_reg_cluster
)
return(de_genes_table)
}
然后用我們最開始下載的RPKM矩陣
> load(file="female_rpkm.Rdata")
> de_clusters <- get_up_reg_clusters(
females,
female_clustering,
female_DE_genes
)
（六）差異基因功能注釋
之前用了monocle分析得到了差異基因，接下來就是對(duì)這些差異基因進(jìn)行注釋了翻翩。首先我們要先找到這些差異基因的基因名：
> de_genes <- subset(de_clusters, qval<0.05)
> length(de_genes$genes)
[1] 4435
基因的ID進(jìn)行轉(zhuǎn)換：
基因ID轉(zhuǎn)換
> entrez_genes <- bitr(de_genes$genes, fromType="SYMBOL",
toType="ENTREZID",
OrgDb="org.Mm.eg.db")
NOTE:接下來一步都许，根據(jù)課程（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-4-差異分析及功能注釋）里寫的稻薇，作者剔除掉了一個(gè)基因名
> entrez_genes <- entrez_genes[!entrez_genes$ENTREZID %in% "101055843",]
> length(entrez_genes$SYMBOL)
[1] 4281 #有些基因SYMBOL沒有對(duì)應(yīng)的ENTREZID，所以少了100多個(gè)基因
#把那些存在ENTREZID的基因的基因名和cluster信息提取出來
> de_gene_clusters <- de_genes[de_genes$genes %in% entrez_genes$SYMBOL,
c("genes", "cluster")]
# 保持de_gene_clusters$genes的順序不變胶征，將symbol變成entrez ID
> de_gene_clusters <- data.frame(
ENTREZID=entrez_genes$ENTREZID[entrez_genes$SYMBOL %in% de_gene_clusters$genes],
cluster=de_gene_clusters$cluster
)
現(xiàn)在看一下de_gene_clusters塞椎，可以看出來這4000多個(gè)差異基因的ENTREZID名字都有一個(gè)對(duì)應(yīng)的cluster了：
下面就是把這些混在一起的cluster們都拆分出來，
> list_de_gene_clusters <- split(de_gene_clusters$ENTREZID,
de_gene_clusters$cluster)
然后再看看我們拆分的情況：
下面就是注釋的步驟了：
> formula_res <- compareCluster(
ENTREZID~cluster,
data=de_gene_clusters,
fun="enrichGO",
OrgDb="org.Mm.eg.db",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05
)
最后就是可視化：
在教程里睛低，大神用的是dotplot來畫圖案狠，可奈何我的電腦不知道為何一直報(bào)錯(cuò)。钱雷。骂铁。明明就是很簡(jiǎn)單的一句代碼：
> dotplot(formula_res, showCategory=5)
解決了兩天也沒弄好，后來尋思著那就換一個(gè)包畫圖吧罩抗。从铲。。在網(wǎng)上搜了好久澄暮，發(fā)現(xiàn)了這篇文章：聽說你也在畫dotplot名段，但是我不服! 文章里說“clusterProfiler的輸出你不用as.data.frame，就可以直接傳給ggplot2”,那我就用ggplot2試試吧
> ggplot(formula_res, aes(Cluster, Description), showCategory=5)
+ geom_point(aes(color=qvalue, size=GeneRatio)) + theme_bw()+scale_colour_gradient(low="#132B43",high="#56B1F7")
到此為止泣懊，smartseq2的單細(xì)胞測(cè)序基本的分析就做完了伸辟。因?yàn)槲也蛔雠咛ミ@塊，發(fā)育軌跡接觸的機(jī)會(huì)也幾乎沒有馍刮，所以沒有跟著教程做發(fā)育譜系可視化這一塊練習(xí)（Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-5-發(fā)育譜系推斷及可視化）,有興趣和有需要的同學(xué)可以繼續(xù)練習(xí)這一塊信夫。總之卡啰，需要深入學(xué)習(xí)的還有很多静稻，這次練習(xí)的收獲就是初次接觸了Seurat包，也是單細(xì)胞測(cè)序分析的主流包匈辱。對(duì)于這些經(jīng)常能用得到的包振湾，還是多研究一下比較好，最好是能去官網(wǎng)閱讀一下說明書亡脸，加深理解押搪。
