劉小澤寫于19.8.8-第三單元前四講：學(xué)習(xí)scRNAseq這個(gè)R包
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

前言

內(nèi)容在：https://github.com/jmzeng1314/scRNA_smart_seq2/blob/master/scRNA/study_scRNAseq.html

要使用scRNAseq這個(gè)R包悦穿，首先要對(duì)它進(jìn)行了解，包中內(nèi)置了Pollen et al. 2014 的數(shù)據(jù)集（https://www.nature.com/articles/nbt.2967）待笑，到19年8月為止滥崩，已經(jīng)有446引用量了珊拼。只不過(guò)原文完整的數(shù)據(jù)是 23730 features扣泊， 301 samples，這個(gè)包中只選取了4種細(xì)胞類型：pluripotent stem cells 分化而成的 neural progenitor cells (NPC催束，神經(jīng)前體細(xì)胞) 集峦，還有 GW16（radial glia，放射狀膠質(zhì)細(xì)胞） 抠刺、GW21（newborn neuron塔淤，新生兒神經(jīng)元） 、GW21+3(maturing neuron矫付，成熟神經(jīng)元) 凯沪，它們的關(guān)系如下圖（NPC和其他三類存在較大差別）：

要想知道數(shù)據(jù)怎么處理的，可以看：https://hemberg-lab.github.io/scRNA.seq.datasets/human/tissues/

簡(jiǎn)單看看文章怎么說(shuō)

文章題目是：Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex

粗略看一下买优，不要求全文通讀妨马。作者利用低覆蓋度單細(xì)胞轉(zhuǎn)錄組測(cè)序，揭示了大腦皮層發(fā)育過(guò)程的細(xì)胞異質(zhì)性和激活的信號(hào)通路

研究背景

大規(guī)模的單細(xì)胞表達(dá)譜測(cè)序具有鑒定罕見細(xì)胞類型和發(fā)育關(guān)系的潛力杀赢，但需要有效地細(xì)胞捕獲和mRNA測(cè)序方法烘跺。目前使用cell barcoding技術(shù)可以實(shí)現(xiàn)極低深度的并行測(cè)序，但是這種低深度測(cè)序有沒有什么弊端還不知道脂崔。

文章使用了11個(gè)細(xì)胞群體的301個(gè)細(xì)胞進(jìn)行了低深度測(cè)序(大約每個(gè)細(xì)胞測(cè)50000條reads)滤淳，發(fā)現(xiàn)這種方法也可以和高深度一樣進(jìn)行細(xì)胞類型鑒定和biomarker鑒定。

材料方法

• 分為高深度測(cè)序和低深度測(cè)序砌左，總共301個(gè)細(xì)胞
  高深度：~8.9 × 10⁶ reads per cell
  低深度：~2.7 x 10⁵ reads per cell

• 數(shù)據(jù)在： http://www.ncbi.nlm.nih.gov/sra/?term=SRP041736

結(jié)論一：低覆蓋度和高覆蓋度的測(cè)序結(jié)果沒太大區(qū)別

論點(diǎn)一：利用低覆蓋度得到的spike-in含量和它們已知的濃度相關(guān)性很高(r = 0脖咐，968)铺敌，當(dāng)每次反應(yīng)有32拷貝以上的spike-in時(shí)，所有的spike-in在所有樣本中都能檢測(cè)出屁擅，并且差異不大

論點(diǎn)二：利用高覆蓋度檢測(cè)到的大部分基因偿凭，在低覆蓋度方法中也能檢測(cè)到

論點(diǎn)三：只在高覆蓋度方法中檢測(cè)到的基因，98%不是高表達(dá)(TPM>100)派歌，大多數(shù)(63%)是低表達(dá)(1<TPM<10)

論點(diǎn)四：對(duì)于不同來(lái)源的的301個(gè)細(xì)胞弯囊，低覆蓋度和高覆蓋度得到的不同基因表達(dá)量估值的平均相關(guān)性為0.91

疑問(wèn)一：但是呢，利用低覆蓋度方法測(cè)序得到的低表達(dá)量基因(1<TPM<10)胶果，它和高覆蓋度的相關(guān)性掉到了0.25匾嘱，也就是說(shuō)，在檢測(cè)低表達(dá)基因方面早抠，低覆蓋度方法還是有局限性的

論點(diǎn)五： 雖然有局限霎烙，但是就取樣技術(shù)來(lái)說(shuō)，使用微流控(Microfluidic)獲得的10個(gè)K562細(xì)胞和使用流式細(xì)胞術(shù)(flow cytometry)得到的一大群細(xì)胞贝或，它們得到的表達(dá)量相關(guān)性很強(qiáng)(r = 0.955)

結(jié)論二: 低覆蓋度測(cè)序也能區(qū)分不同細(xì)胞類型

論點(diǎn)一： 先利用PCA看看細(xì)胞分群的效果

論點(diǎn)二：看看不同分組的基因表達(dá)情況

論點(diǎn)三： 低覆蓋度和高覆蓋度得到的細(xì)胞分布很相似吼过，而且PCA得到的貢獻(xiàn)最大的500個(gè)基因中，有78%是在低覆蓋度和高覆蓋度中都存在的

然后主要看這個(gè)包中的數(shù)據(jù)

第一次見不會(huì)怎么辦咪奖？看幫助文檔和Bioconductor相關(guān)包的教程

# 先看幫助文檔
library(scRNAseq)
?scRNAseq
# 其中包含了三種數(shù)據(jù)集：fluidigm、th2酱床、allen羊赵，我們用到的是第一個(gè)fluidigm
# The dataset fluidigm contains 65 cells from Pollen et al. (2014), each sequenced at high and low coverage (SRA: SRP041736).

也就是說(shuō)，雖然原文總共做了301個(gè)樣本扇谣，但這里一共有130個(gè)樣本文庫(kù)(高覆蓋度昧捷、低覆蓋度各65個(gè))

具體數(shù)據(jù)的處理可以看這個(gè)R包的Bioconductor詳細(xì)文檔：https://bioconductor.org/packages/release/data/experiment/vignettes/scRNAseq/inst/doc/scRNAseq.html

其中介紹了另外兩個(gè)數(shù)據(jù)集的來(lái)歷，瀏覽一下：

• The dataset th2 contains 96 T helper cells from (Mahata et al. 2014) (ArrayExpress: E-MTAB-2512).
• The dataset allen contains 379 cells from the mouse visual cortex. This is a subset of the data published in (Tasic et al. 2016) (SRA: SRP061902).

還介紹了數(shù)據(jù)的預(yù)處理：

直接看看數(shù)據(jù)應(yīng)該怎么獲裙拚靡挥？

這個(gè)包中的數(shù)據(jù)都是以SummarizedExperiment對(duì)象形式存放的，那么什么是SummarizedExperiment對(duì)象鸯绿？

使用?assay得到的幫助結(jié)果：

The SummarizedExperiment class is a matrix-like container where rows represent features of interest (e.g. genes, transcripts, exons, etc...) and columns represent samples (with sample data summarized as a DataFrame).

以第一個(gè)數(shù)據(jù)fluidigm為例進(jìn)行讀劝掀啤：

library(scRNAseq)
data(fluidigm)

> fluidigm
class: SummarizedExperiment 
dim: 26255 130 
metadata(3): sample_info clusters which_qc
assays(4): tophat_counts cufflinks_fpkm rsem_counts rsem_tpm
rownames(26255): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
rowData names(0):
colnames(130): SRR1275356 SRR1274090 ... SRR1275366 SRR1275261
colData names(28): NREADS NALIGNED ... Cluster1 Cluster2

下面就是對(duì)這個(gè)對(duì)象的探索了：

# 例如要提取基因表達(dá)量的信息，就用assay函數(shù)(注意上面??assays那一行瓶蝴，其中包含了4個(gè)結(jié)果：tophat_counts毒返、 cufflinks_fpkm、rsem_counts舷手、rsem_tpm)
names(assays(fluidigm))
## [1] "tophat_counts"  "cufflinks_fpkm" "rsem_counts"    "rsem_tpm"

# 默認(rèn)訪問(wèn)第一個(gè)拧簸，也就是原始的表達(dá)量tophat_counts
head(assay(fluidigm)[,1:3])
##          SRR1275356 SRR1274090 SRR1275251
## A1BG              0          0          0
## A1BG-AS1          0          0          0
## A1CF              0          0          0
## A2M               0          0          0
## A2M-AS1           0          0          0
## A2ML1             0          0          0

# 如果要得到RPKM值，可以使用assay：
head(assay(fluidigm, 2)[,1:3])
##          SRR1275356 SRR1274090 SRR1275251
## A1BG              0  0.0000000          0
## A1BG-AS1          0  0.3256690          0
## A1CF              0  0.0687904          0
## A2M               0  0.0000000          0
## A2M-AS1           0  0.0000000          0
## A2ML1             0  1.3115300          0

# 或者使用assays
head(assays(fluidigm)$cufflinks_fpkm)

看完表達(dá)矩陣男窟，少不了的是樣本的注釋信息盆赤，這些就存放在了：colData中

# 包含了太多的信息贾富，如果你直接使用colData(fluidigm)，會(huì)得到眼花繚亂的結(jié)果
# 于是可以先大體看看有哪些類
names(metadata(fluidigm))
## [1] "sample_info" "clusters"    "which_qc"

# 然后假如我們想看QC相關(guān)的信息(也是最常用的)
metadata(fluidigm)$which_qc
##  [1] "NREADS"                       "NALIGNED"                    
##  [3] "RALIGN"                       "TOTAL_DUP"                   
##  [5] "PRIMER"                       "INSERT_SZ"                   
##  [7] "INSERT_SZ_STD"                "COMPLEXITY"                  
##  [9] "NDUPR"                        "PCT_RIBOSOMAL_BASES"         
## [11] "PCT_CODING_BASES"             "PCT_UTR_BASES"               
## [13] "PCT_INTRONIC_BASES"           "PCT_INTERGENIC_BASES"        
## [15] "PCT_MRNA_BASES"               "MEDIAN_CV_COVERAGE"          
## [17] "MEDIAN_5PRIME_BIAS"           "MEDIAN_3PRIME_BIAS"          
## [19] "MEDIAN_5PRIME_TO_3PRIME_BIAS"

# 因此我們想獲得樣本QC信息牺六，就可以
sample_qc <- as.data.frame(colData(fluidigm)[metadata(fluidigm)$which_qc])

探索完祷安，開始基本操作

# 我們要對(duì)RSEM得到的count值進(jìn)行操作，之所以使用floor函數(shù)兔乞，是因?yàn)檫@個(gè)RSEM矩陣存在小數(shù)點(diǎn)汇鞭。猜測(cè)：因?yàn)镽SEM計(jì)算表達(dá)量是考慮了reads比對(duì)到不同基因的情況，這樣的話就不能直接判斷這個(gè)reads到底屬于哪個(gè)基因庸追，于是就用帶小數(shù)的expected count(也就是真實(shí)值)表示霍骄。其實(shí)我們使用的時(shí)候，是需要變成整數(shù)（raw count）的淡溯，于是簡(jiǎn)單使用了floor向下取整
mtx <- floor(assay(fluidigm,3))

> mtx[1:3,1:3]
         SRR1275356 SRR1274090 SRR1275251
A1BG              0          0          0
A1BG-AS1          0          0          0
A1CF              0          0          0
> dim(mtx)
[1] 26255   130

看表型信息并過(guò)濾

想法是：對(duì)每個(gè)QC指標(biāo)都做個(gè)箱線圖读整，這些指標(biāo)會(huì)以向量的形式保存，然后一個(gè)一個(gè)循環(huán)操作咱娶，那么會(huì)向量循環(huán)就用lapply

# 目前QC指標(biāo)都存在：colnames(sample_qc)中米间，一共19個(gè)；使用更容易調(diào)參數(shù)的ggboxplot
library(ggpubr)
box <- lapply(colnames(sample_qc),function(i) {
  dat <-  sample_qc[,i,drop=F] 
  dat$all_cells="all_cells"
  ggboxplot(dat,x=dat[,2],y=i,
            xlab=F,add = "jitter")
})
plot_grid(plotlist=box, ncol=5 )

然后利用表型信息對(duì)樣本進(jìn)行過(guò)濾：

# 從QC數(shù)據(jù)中挑選一些指標(biāo)膘侮，作為過(guò)濾條件
choose_anno <-  colnames(sample_qc[,c(1:9,11:16,18,19)])

# 下面就是將一個(gè)個(gè)的QC條件進(jìn)行細(xì)胞的過(guò)濾屈糊，如果細(xì)胞滿足設(shè)定的QC過(guò)濾條件，就為1琼了；否則為0逻锐；并且用cbind按列組合在一起
filter <- lapply(choose_anno,function(i) {
  # 寫循環(huán)時(shí)可以先用一個(gè)值作為測(cè)試：例如 i=choose_anno[1]
  dat <-  sample_qc[,i]  
  dat <- abs(log10(dat))
  fivenum(dat)
  (up <- mean(dat)+2*sd(dat))
  (down <- mean(dat)- 2*sd(dat) ) 
  valid <- ifelse(dat > down & dat < up, 1,0 ) 
})
filter <- do.call(cbind,filter)

# 得到了列為QC，行為細(xì)胞的過(guò)濾結(jié)果雕薪，那么就將QC條件全部為1的細(xì)胞挑出來(lái)(也就是找全部為1的行)
choosed_cells <- apply(filter,1,function(x) all(x==1))

# 進(jìn)行對(duì)比：原來(lái)的細(xì)胞
> table(colData(fluidigm)$Biological_Condition)

  GW16   GW21 GW21+3    NPC 
    52     16     32     30 
# 過(guò)濾后的細(xì)胞
> table(colData(fluidigm)[choosed_cells,]$Biological_Condition)

  GW16   GW21 GW21+3    NPC 
    36     11     23     29

# 將表達(dá)矩陣進(jìn)行過(guò)濾
mtx <- mtx[,choosed_cells]
> dim(mtx)
[1] 26255    99

看基因表達(dá)信息

> fivenum(apply(mtx,1,function(x) sum(x>0) ))
     A1CF     OR8G1 LINC01003    MRPS36     YWHAZ 
        0         0         4        26        99 
# 看到至少有25%的基因表達(dá)量為0.那么具體有多少個(gè)呢昧诱？可以看看：
choosed_genes=apply(mtx,1,function(x) sum(x>0) )>0
> table(choosed_genes)

FALSE  TRUE 
 9496 16759 
# 看到有9000多個(gè)基因在所有細(xì)胞中都沒有表達(dá)量(可能原因：原文分析的確實(shí)是2w多個(gè)基因，但他是在300多個(gè)細(xì)胞中分析的所袁，我們這個(gè)包里過(guò)濾后只剩下99個(gè)細(xì)胞盏档，所以存在很多基因不在這部分細(xì)胞中表達(dá)，因此需要去掉)
boxplot(apply(mtx,1,function(x) sum(x>0) ))
              
# 最后根據(jù)基因表達(dá)量對(duì)矩陣進(jìn)行過(guò)濾
mtx <- mtx[choosed_genes,]

從下面這個(gè)箱線圖中也可以看到燥爷，很少有基因在99個(gè)過(guò)濾后的細(xì)胞中都有表達(dá)蜈亩，大部分基因還是在部分細(xì)胞中表達(dá)量為0 (這也是單細(xì)胞一個(gè)很特殊的現(xiàn)象：dropout情況，意思就是真實(shí)情況下基因是有表達(dá)量的局劲，但技術(shù)問(wèn)題沒有檢測(cè)到)

看細(xì)胞間基因表達(dá)量相關(guān)性

# 都要基于CPM標(biāo)準(zhǔn)化數(shù)值勺拣，并做一個(gè)備份
mtx <- log2(edgeR::cpm(mtx) + 1)
mtx[1:4, 1:4]
mtx_back <- mtx

# 對(duì)相關(guān)性進(jìn)行初步的可視化
exprSet <- mtx_back
> dim(exprSet)
[1] 16759    99

pheatmap::pheatmap(cor(exprSet))
# 注意：cor函數(shù)計(jì)算的是列與列間的相關(guān)系數(shù)

需要加上分組的信息：

# 使用細(xì)胞過(guò)濾后的分組信息(GW16：36  GW21：11  GW21+3：23  NPC：29）         
group_list <- colData(fluidigm)[choosed_cells,]$Biological_Condition
tmp <- data.frame(g = group_list)
rownames(tmp) <-  colnames(exprSet)
# 組內(nèi)相似性高于組間，并且看到NPC組和其他組差異更大
pheatmap::pheatmap(cor(exprSet), annotation_col = tmp)

好鱼填，接著設(shè)置閾值對(duì)表達(dá)矩陣過(guò)濾

# 這一次設(shè)置閾值為5药有，表示至少要滿足基因在5個(gè)細(xì)胞中的表達(dá)量都大于1
exprSet = exprSet[apply(exprSet, 1, function(x) sum(x > 1) > 5), ]
> dim(exprSet)
[1] 11337    99

過(guò)濾完，按照mad統(tǒng)計(jì)方法取前500個(gè)表達(dá)量變化最大的基因

# 絕對(duì)中位差來(lái)估計(jì)方差,先計(jì)算出數(shù)據(jù)與它們的中位數(shù)之間的偏差，然后這些偏差的絕對(duì)值的中位數(shù)就是mad
exprSet <-  exprSet[names(sort(apply(exprSet, 1, mad), decreasing = T)[1:500]), ]
> dim(exprSet)
[1] 500  99

# 對(duì)組間差異最大的基因再進(jìn)行相關(guān)性分析
M <-cor(log2(exprSet + 1))
tmp <- data.frame(g = group_list)
rownames(tmp) <-  colnames(M)
pheatmap::pheatmap(M, annotation_col = tmp)

小結(jié)：NPC跟另外的GW細(xì)胞群可以區(qū)分的很好愤惰，但是GW本身的3個(gè)小群體并沒有那么好的區(qū)分度苇经。后來(lái)簡(jiǎn)單選取mad前500的基因重新計(jì)算，也沒有改善宦言；另外可以看到每個(gè)細(xì)胞測(cè)了兩次(圖中對(duì)角線中有紅色和橙色扇单，表示兩次不同深度)

表達(dá)矩陣簡(jiǎn)單的層次聚類

mtx <- mtx_back
hc <- hclust(dist(t(mtx))) # dist以行為輸入
plot(hc,labels = FALSE)

clus <-  cutree(hc, 4) #對(duì)hclust()函數(shù)的聚類結(jié)果進(jìn)行剪枝，即選擇輸出指定類別數(shù)的系譜聚類結(jié)果奠旺。
group_list <-  as.factor(clus) ##轉(zhuǎn)換為因子屬性
> table(group_list) ##統(tǒng)計(jì)頻數(shù)
group_list
 1  2  3  4 
29 25 39  6 

filtered_anno <- colData(fluidigm)[choosed_cells,]$Biological_Condition

> table(group_list,filtered_anno)   
          filtered_anno
group_list GW16 GW21 GW21+3 NPC
         1    0    0      0  29
         2   20    3      2   0
         3   15    8     16   0
         4    1    0      5   0
# 結(jié)果看到：NPC蜘澜、GW21+3這兩個(gè)利用普通的層次聚類還是可以區(qū)分開，但是GW16响疚、GW21就不太能區(qū)分了

最常規(guī)的PCA降維結(jié)果

算法很多鄙信，比如：主成分分析PCA、多維縮放(MDS)忿晕、線性判別分析(LDA)装诡、等度量映射(Isomap)、局部線性嵌入(LLE)践盼、t-SNE鸦采、Deep Autoencoder Networks
以下會(huì)采用PCA 和 t-SNE

mtx <- mtx_back
mtx <- t(mtx) # PCA也是對(duì)行操作：需要先轉(zhuǎn)置一下，讓行為樣本
mtx <- as.data.frame(mtx)
plate <- filtered_anno # 這里定義分組信息
mtx <-  cbind(mtx, plate) # 添加分組信息

> mtx[1:4, 1:4]
           A1BG A1BG-AS1      A2M A2M-AS1
SRR1274090    0        0 0.000000       0
SRR1275287    0        0 4.216768       0
SRR1275364    0        0 0.000000       0
SRR1275269    0        0 3.552694       0
> table(mtx$plate)

  GW16   GW21 GW21+3    NPC 
    36     11     23     29 

# 進(jìn)行PCA操作(實(shí)際是降維咕幻，映射到二維坐標(biāo))
mtx.pca <- PCA(mtx[, -ncol(mtx)], graph = FALSE)

> head(mtx.pca$var$coord) ## 每個(gè)主成分的基因重要性占比
               Dim.1       Dim.2       Dim.3        Dim.4         Dim.5
A1BG      0.19046450  0.09601240 -0.17840553 -0.001507970 -0.0006057691
A1BG-AS1 -0.02510451  0.29821319  0.03571804  0.020001929 -0.0105727109
A2M       0.03403042  0.25458727  0.24264958  0.228512329  0.5414019044
A2M-AS1   0.23140893  0.02900348 -0.07952678  0.356461354  0.1283450099
A2ML1    -0.15776536  0.13831288  0.10065788  0.004060288 -0.0353422367
A2MP1    -0.04068586 -0.05584736 -0.02857416  0.018287992  0.0069603680
> head(mtx.pca$ind$coord) ## 每個(gè)細(xì)胞的前5個(gè)主成分取值渔伯。
                Dim.1      Dim.2      Dim.3      Dim.4      Dim.5
SRR1274090  40.251912 -13.231641 -12.358891 -20.038100 -12.704947
SRR1275287   1.196637  15.386256  30.566235  14.262858  -4.852418
SRR1275364 -34.731051 -14.782146  -7.716928   7.046918   1.951473
SRR1275269  21.760471   3.307309  17.985263 -18.382512   9.270646
SRR1275263  -3.313968 -15.856721   8.929275 -36.358830  20.275875
SRR1274117  59.378486  16.453551  -5.098901  56.245455  19.257598

PCA函數(shù)返回的結(jié)果如下：現(xiàn)在細(xì)胞在行，基因在列谅河，所以細(xì)胞是ind咱旱，基因是var

# 進(jìn)行PCA可視化
fviz_pca_ind(
  mtx.pca,
  #repel =T,
  geom.ind = "point",
  # show points only (nbut not "text")
  col.ind = mtx$plate,
  # color by groups
  #palette = c("#00AFBB", "#E7B800"),
  addEllipses = TRUE,
  # Concentration ellipses
  legend.title = "Groups"
) 

小結(jié)：NPC跟另外的GW細(xì)胞群可以區(qū)分的很好，但是GW本身的3個(gè)小群體并沒有那么好的區(qū)分度

進(jìn)階一點(diǎn)的tSNE降維

這里先選取PCA后的主成分绷耍，然進(jìn)行tSNE；其實(shí)還可以選取變化高的基因鲜侥，顯著差異基因等等

# 選取前面PCA分析的5個(gè)主成分褂始。
tsne_mtx <- mtx.pca$ind$coord
# Set a seed if you want reproducible results
set.seed(42)
library(Rtsne) 
# 如果使用原始表達(dá)矩陣進(jìn)行 tSNE耗時(shí)會(huì)很久
# 如果出現(xiàn)Remove duplicates before running TSNE 則check_duplicated = FALSE
# tsne_out <- Rtsne(dat_matrix,pca=FALSE,perplexity=30,theta=0.0, check_duplicates = FALSE) 

# Run TSNE
tsne_out <- Rtsne(tsne_mtx,perplexity=10)
plate <- filtered_anno # 這里定義分組信息
plot(tsne_out$Y,col= rainbow(4)[as.numeric(as.factor(plate))], pch=19) 

降維后呢？

降維和聚類不是一回事描函，各有各的算法崎苗、參數(shù)，比如降維我們常用PCA舀寓、tsne胆数，聚類就有kmeans、dbscan

> # 前面我們的層次聚類是針對(duì)全部表達(dá)矩陣互墓，tsne選取前面PCA分析的5個(gè)主成分必尼，因此為了節(jié)省計(jì)算量，選取tsne_out$Y這個(gè)結(jié)果
> head(tsne_out$Y)
          [,1]        [,2]
[1,]  4.855236 -26.9704714
[2,]  0.179925  -0.5475169
[3,]  6.256713  24.9241040
[4,]  2.471635 -20.2250523
[5,]  2.615960 -12.6056267
[6,] -2.384375 -22.3821087

opt_tsne=tsne_out$Y

> table(kmeans(opt_tsne,centers = 4)$clust)

 1  2  3  4 
31 24 24 20 

plot(opt_tsne,  col=kmeans(opt_tsne,centers = 4)$clust, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")

# 換一種dbscan
library(dbscan)
plot(opt_tsne,  col=dbscan(opt_tsne,eps=3.1)$cluster, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")

> table(dbscan(opt_tsne,eps=3.1)$cluster)

 0  1  2  3  4 
 2 22 38 31  6 

進(jìn)行一個(gè)比較：

# 比較兩個(gè)聚類算法區(qū)別
> table(kmeans(opt_tsne,centers = 4)$clust,dbscan(opt_tsne,eps=3.1)$cluster)
   
     0  1  2  3  4
  1  0  0 38  0  6
  2  0  0  0 16  0
  3  0  0  0 15  0
  4  2 22  0  0  0

下面使用M3Drop包處理單細(xì)胞數(shù)據(jù)

第一步 構(gòu)建對(duì)象

## 重新加載數(shù)據(jù)
rm(list=ls())
data(fluidigm)
# names(assays(fluidigm))
counts <- floor(assay(fluidigm, 3))
dim(counts)

## 過(guò)濾
sample_qc <- as.data.frame(colData(fluidigm)[metadata(fluidigm)$which_qc])
choose_anno <-  colnames(sample_qc[,c(1:9,11:16,18,19)])
filter <- lapply(choose_anno,function(i) {
  # i=choose_anno[1]
  dat <-  sample_qc[,i]  
  dat <- abs(log10(dat))
  fivenum(dat)
  (up <- mean(dat)+2*sd(dat))
  (down <- mean(dat)- 2*sd(dat) ) 
  valid <- ifelse(dat > down & dat < up, 1,0 ) 
})

filter <- do.call(cbind,filter)
choosed_cells <- apply(filter,1,function(x) all(x==1))

counts <- counts[,choosed_cells]

## 開始M3Drop分析
library(M3Drop) 
Normalized_data <- M3DropCleanData(counts, 
                                   labels = colData(fluidigm)[choosed_cells,]$Biological_Condition, 
                                   is.counts=TRUE, min_detected_genes=2000)

> dim(Normalized_data$data)
[1] 13405    97
                       
## 檢查
> str(Normalized_data) #只返回了一個(gè)list，而不是S4對(duì)象
List of 2
 $ data  : num [1:13405, 1:97] 0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:13405] "A1BG" "A2M" "A2ML1" "AAAS" ...
  .. ..$ : chr [1:97] "SRR1274090" "SRR1275287" "SRR1275364" "SRR1275269" ...
 $ labels: chr [1:97] "NPC" "GW21+3" "GW16" "GW21" ...

第二步 表達(dá)矩陣檢驗(yàn)--Michaelis-Menten算法

具體的算法可以不用了解太深

fits <- M3DropDropoutModels(Normalized_data$data)

## 查看檢驗(yàn)的結(jié)果
## Sum absolute residuals
data.frame(MM=fits$MMFit$SAr, Logistic=fits$LogiFit$SAr,
           DoubleExpo=fits$ExpoFit$SAr) 
#    MM Logistic DoubleExpo
1 1651     1646       4033


## Sum absolute residuals
data.frame(MM=fits$MMFit$SAr, Logistic=fits$LogiFit$SAr,
           DoubleExpo=fits$ExpoFit$SAr) 
#   MM Logistic DoubleExpo
1 403      345       1962

第三步 找差異基因

這里需要注意：如果提示找不到M3DropDifferentialExpression這個(gè)函數(shù)判莉，那么可能由于安裝的M3Drop包版本是舊版豆挽，新版對(duì)應(yīng)的函數(shù)是：M3DropFeatureSelection

DE_genes <-  M3DropFeatureSelection(Normalized_data$data, 
                                         mt_method="fdr", mt_threshold=0.01)
> dim(DE_genes)
[1] 182   4
> head(DE_genes)
           Gene effect.size      p.value      q.value
IGFBPL1 IGFBPL1    13.94805 1.790836e-21 4.801233e-18
DLK1       DLK1    11.66363 1.104816e-22 4.936685e-19
BCAT1     BCAT1    10.97600 6.335191e-31 8.492323e-27
CA12       CA12    10.94265 1.020863e-08 3.110152e-06
ZNF30     ZNF30    10.78263 2.503772e-06 3.813985e-04
FOS         FOS    10.01977 1.389763e-12 9.314889e-10

第四步 差異基因畫熱圖

# 就是看一下差異基因在不同細(xì)胞類型的表達(dá)分布
par(mar=c(1,1,1,1)) 
heat_out <- M3DropExpressionHeatmap(DE_genes$Gene, Normalized_data$data, 
                                    cell_labels = Normalized_data$labels)

第五步 重新聚類，找marker

同樣的券盅，如果函數(shù)M3DropGetHeatmapCellClusters找不到帮哈，就替換成M3DropGetHeatmapClusters，因?yàn)樽髡邥r(shí)刻在改函數(shù)名字锰镀，有困難找：https://bioconductor.org/packages/release/bioc/vignettes/M3Drop/inst/doc/M3Drop_Vignette.R

# 重新聚類
cell_populations <- M3DropGetHeatmapClusters(heat_out, k=4)
# 重聚類得到的和自帶的表型數(shù)據(jù)比較
> table(cell_populations,Normalized_data$labels)
                
cell_populations GW16 GW21 GW21+3 NPC
               1    0    0      0  29
               2   14    8     19   0
               3    4    1      2   0
               4   16    2      2   0

# 找marker
library("ROCR") 
marker_genes <- M3DropGetMarkers(Normalized_data$data, cell_populations)

第六步 可以看每個(gè)分類的marker基因

# 比如想看新得到的第4組的marker基因
> head(marker_genes[marker_genes$Group==4,],10) 
             AUC Group         pval
ADGRV1 0.9707792     4 1.831217e-11
TFAP2C 0.9451299     4 1.885637e-11
EGR1   0.9409091     4 1.159852e-09
PLCE1  0.9233766     4 7.676760e-11
FOS    0.9058442     4 1.806229e-08
SLC1A3 0.9048701     4 1.542660e-09
AASS   0.9032468     4 3.136961e-08
ITGB8  0.8886364     4 8.823909e-08
BCAN   0.8844156     4 1.513264e-12
NFIA   0.8831169     4 7.724889e-08

# 或者想知道某個(gè)marker基因的分配位置
> marker_genes[rownames(marker_genes)=="FOS",] 
          AUC Group         pval
FOS 0.9058442     4 1.806229e-08

當(dāng)然娘侍，可以挑選一些marker基因進(jìn)行作圖

# 挑選的代碼
choosed_marker_genes=as.character(unlist(lapply(split(marker_genes,marker_genes$Group), function(x) (rownames(head(x,20))))))

# 看不懂？沒關(guān)系泳炉，可以拆解~
# 其中最核心的是split(marker_genes,marker_genes$Group)憾筏，它返回什么東西，用str()查看就對(duì)了
> str(split(marker_genes,marker_genes$Group))
List of 4
 $ 1:'data.frame':  5370 obs. of  3 variables:
  ..$ AUC  : num [1:5370] 1 0.998 0.998 0.997 0.993 ...
  ..$ Group: chr [1:5370] "1" "1" "1" "1" ...
  ..$ pval : num [1:5370] 8.68e-22 1.20e-17 1.04e-14 1.18e-14 1.81e-14 ...
 $ 2:'data.frame':  1844 obs. of  3 variables:
  ..$ AUC  : num [1:1844] 0.99 0.936 0.933 0.927 0.927 ...
  ..$ Group: chr [1:1844] "2" "2" "2" "2" ...
# 還有兩段省略............

# 不難看出胡桃，它的作用是根據(jù)分組信息（marker_genes$Group）踩叭，將marker_genes分成了4個(gè)部分，共同組成一個(gè)列表〈湟龋現(xiàn)在我們可以最粗略地選出每個(gè)分組中的前20基因(當(dāng)然容贝，后續(xù)可以自定義，根據(jù)AUC和P值來(lái)挑選)
# 對(duì)列表循環(huán)用lapply之景，要做的事情就是挑出前20=》rownames(head(x,20))
                                                
# 接下來(lái)作圖就容易了
par(mar=c(1,1,1,1)) 
heat_out <- M3DropExpressionHeatmap(choosed_marker_genes, Normalized_data$data, cell_labels =  cell_populations)