【單細(xì)胞轉(zhuǎn)錄組 實戰(zhàn)】七、scRNAseq、M3Drop Package

這里是佳奧！我們補充了相關(guān)知識，現(xiàn)在進一步學(xué)習(xí)相關(guān)R包的使用方法。

文章數(shù)據(jù)：Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex

簡要的流程在Git Hub下載的壓縮包中，section02-scRNA文件夾的study_scRNAseq文件。

QQ截圖20220901134733.png

打開R Project文件，我們正式開始吧。

注：從這里開始，難度會明顯提高，先從運行代碼開始，再通過幫助文檔和參考代碼逐漸理解，最后理解下列代碼。這需要時間。

1 R包安裝與鏡像設(shè)置

rm(list = ls())
Sys.setenv(R_MAX_NUM_DLLS=999) ##Sys.setenv修改環(huán)境設(shè)置，R的namespace是有上限的，如果導(dǎo)入包時超過這個上次就會報錯,R_MAX_NUM_DLLS可以修改這個上限
options(stringsAsFactors = F) ##options:允許用戶對工作空間進行全局設(shè)置，stringsAsFactors防止R自動把字符串string的列辨認(rèn)成factor

options()$repos  ## 查看使用install.packages安裝時的默認(rèn)鏡像
options()$BioC_mirror ##查看使用bioconductor的默認(rèn)鏡像
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") ##指定鏡像，這個是中國科技大學(xué)鏡像
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) ##指定install.packages安裝鏡像，這個是清華鏡像
options()$repos 
options()$BioC_mirror

# http://www.bio-info-trainee.com/3727.html 周末生物信息學(xué)培訓(xùn)班準(zhǔn)備工作
 
if (!requireNamespace("BiocManager", quietly = TRUE)) 
  install.packages("BiocManager") ##判斷是否存在BiocManager包，不存在的話安裝

if(!require('Seurat')){
  BiocManager::install('Seurat',ask = F,update = F)
}
if(!require('scran')){
  BiocManager::install(c( 'scran'),ask = F,update = F)
}
if(!require('monocle')){
  BiocManager::install(c( 'monocle'),ask = F,update = F)
}

## http://cole-trapnell-lab.github.io/monocle-release/monocle3/ 
BiocManager::install('destiny')
BiocManager::install(c( 'flexclust','mcclust'),ask = F,update = F)

BiocManager::install(c( 'scRNAseq'),ask = F,update = F)
BiocManager::install(c( 'dbscan'),ask = F,update = F)
BiocManager::install(c( 'M3Drop'),ask = F,update = F)

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(reticulate))
suppressMessages(library(devtools))
suppressMessages(library(monocle))
suppressMessages(library(flexclust))
suppressMessages(library(mcclust))

全部不報錯后，我們開始第一個包的學(xué)習(xí)。

2 scRNAseq Package

需要自行下載安裝一些必要的R包，主要是scRNAseq，然后是其它一些輔助包用來探索這個數(shù)據(jù)集。

if (!requireNamespace("Rtsne"))
    install.packages("Rtsne")
if (!requireNamespace("FactoMineR"))
    install.packages("FactoMineR")
if (!requireNamespace("factoextra"))
    install.packages("factoextra")
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
if (!requireNamespace("scater"))
    BiocManager::install("scater")
if (!requireNamespace("scRNAseq"))
    BiocManager::install("scRNAseq") 
if (!requireNamespace("M3Drop"))
    BiocManager::install("M3Drop") 
if (!requireNamespace("ROCR"))
    BiocManager::install("ROCR") 

##加載R包，前提是已經(jīng)成功安裝了
rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(scater))
suppressMessages(library(scRNAseq))
library(ggplot2)
library(tidyr)
library(cowplot)
library("FactoMineR")
library("factoextra")
library("ROCR")

scRNAseq R包中的數(shù)據(jù)集

這個包內(nèi)置的是 Pollen et al. 2014 數(shù)據(jù)集，人類單細(xì)胞細(xì)胞，分成4類，分別是 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，還有 “GW16” and “GW21” ，“GW21+3” 這種孕期細(xì)胞，理解這些需要一定的生物學(xué)背景知識，如果不感興趣，可以略過。

不過本例子只使用了數(shù)據(jù)集的4種細(xì)胞類型而已，因為 scRNAseq 這個R包就提供了這些，完整的數(shù)據(jù)是 23730 features，301 samples, 地址為：

https://hemberg-lab.github.io/scRNA.seq.datasets/human/tissues/

這里面的表達(dá)矩陣是由 RSEM (Li and Dewey 2011) 軟件根據(jù) hg38 RefSeq transcriptome 得到的，總是130個文庫，每個細(xì)胞測了兩次，測序深度不一樣。

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm) 
fluidigm<-ReprocessedFluidigmData()##作者修改了提取數(shù)據(jù)的方法，需要多運行這一步
ct <- floor(assays(fluidigm)$rsem_counts)##counts含有小數(shù)，因此需要floor()
ct[1:4,1:4] 
sample_ann <- as.data.frame(colData(fluidigm))
DT::datatable(sample_ann)

先探索表型信息

前面說到，這個數(shù)據(jù)集是130個文庫，每個細(xì)胞測了兩次，測序深度不一樣，這130個細(xì)胞，分成4類，分別是： pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，還有 “GW16” and “GW21” ，“GW21+3” 這種孕期細(xì)胞。

批量，粗略的看一看各個細(xì)胞的一些統(tǒng)計學(xué)指標(biāo)的分布情況：

```{r fig.width=10, fig.height=15}##rmarkdown畫圖技巧
box <- lapply(colnames(sample_ann[,1:19]),function(i) {
    dat <-  sample_ann[,i,drop=F] 
    dat$sample=rownames(dat)
    ## 畫boxplot 
   ggplot(dat, aes('all cells', get(i))) +
          geom_boxplot() +
          xlab(NULL)+ylab(i)
})
plot_grid(plotlist=box, ncol=5 )
# ggsave(file="stat_all_cells.pdf")
QQ截圖20220901145444.png

因為進行了簡單探索，對表型數(shù)據(jù)就有了把握，接下來可以進行一定程度的過濾，因為細(xì)節(jié)太多，這里重點是批量，代碼技巧非常值得學(xué)習(xí)。教程略過了這一部分。

pa <- colnames(sample_ann[,c(1:9,11:16,18,19)])
tf <- lapply(pa,function(i) {
 # i=pa[1]
  dat <-  sample_ann[,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 ) 
})

tf <- do.call(cbind,tf)
choosed_cells <- apply(tf,1,function(x) all(x==1))
table(sample_ann$Biological_Condition)
sample_ann=sample_ann[choosed_cells,]
table(sample_ann$Biological_Condition)
ct <- ct[,choosed_cells]

再探索基因表達(dá)情況

ct[1:4,1:4] 
counts <- ct
fivenum(apply(counts,1,function(x) sum(x>0) ))
boxplot(apply(counts,1,function(x) sum(x>0) ))
fivenum(apply(counts,2,function(x) sum(x>0) ))
hist(apply(counts,2,function(x) sum(x>0) ))
choosed_genes=apply(counts,1,function(x) sum(x>0) )>0
table(choosed_genes)
counts <- counts[choosed_genes,]

##查看一下表達(dá)不為0的基因個數(shù)
> table((apply(counts,1,function(x) sum(x>0)>0 )))

 TRUE 
17096

幾乎沒有基因在所有細(xì)胞都表達(dá)，中位線為35，如果基因在35個細(xì)胞內(nèi)不表達(dá)即可過濾掉（閾值）。

QQ截圖20220901150747.png

接下來要利用自己的常規(guī)轉(zhuǎn)錄組數(shù)據(jù)分析知識。

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

下面的計算，都是基于log后的表達(dá)矩陣。

dat <- log2(edgeR::cpm(counts) + 1)
dat[1:4, 1:4]
dat_back <- dat

先備份這個表達(dá)矩陣，后面的分析都用得上。

exprSet <- dat_back
colnames(exprSet)
pheatmap::pheatmap(cor(exprSet))
group_list <- sample_ann$Biological_Condition
tmp <- data.frame(g = group_list)
rownames(tmp) <- colnames(exprSet)

##組內(nèi)的樣本的相似性應(yīng)該是要高于組間的！
pheatmap::pheatmap(cor(exprSet), annotation_col = tmp)
dim(exprSet)
exprSet = exprSet[apply(exprSet, 1, function(x)
sum(x > 1) > 5), ]
dim(exprSet)
 
dim(exprSet)
exprSet <-  exprSet[names(sort(apply(exprSet, 1, mad), decreasing = T)[1:500]), ]
dim(exprSet)
M <-cor(log2(exprSet + 1))
tmp <- data.frame(g = group_list)
rownames(tmp) <-  colnames(M)
pheatmap::pheatmap(M, annotation_col = tmp)

table(sample_ann$LibraryName)
QQ截圖20220901152639.png

可以看到,從細(xì)胞的相關(guān)性角度來看，到NPC跟另外的GW細(xì)胞群可以區(qū)分的很好，但是GW本身的3個小群體并沒有那么好的區(qū)分度。

而且簡單選取top的sd的基因來計算相關(guān)性，并沒有很明顯的改善。

但是可以看到每個細(xì)胞測了兩次，所以它們的相關(guān)性要好于其它同類型的細(xì)胞。

對表達(dá)矩陣進行簡單的層次聚類

如果計算機資源不夠，這里可以先對基因進行一定程度的挑選，最簡單的就是選取top的sd的基因，這里略。

dat <- dat_back
hc <- hclust(dist(t(dat))) 
plot(hc,labels = FALSE)
clus <-  cutree(hc, 4) #對hclust()函數(shù)的聚類結(jié)果進行剪枝，即選擇輸出指定類別數(shù)的系譜聚類結(jié)果。
group_list <-  as.factor(clus) ##轉(zhuǎn)換為因子屬性
table(group_list) ##統(tǒng)計頻數(shù)
table(group_list,sample_ann$Biological_Condition)
QQ截圖20220901152800.png

可以看到GW16和GW21是很難區(qū)分開來的，如果是普通的層次聚類的話。

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

降維算法很多，詳情可以去自行搜索學(xué)習(xí)，比如：

  1. 主成分分析PCA
  2. 多維縮放(MDS)
  3. 線性判別分析(LDA)
  4. 等度量映射(Isomap)
  5. 局部線性嵌入(LLE)
  6. t-SNE
  7. Deep Autoencoder Networks

這里只介紹 PCA 和 t-SNE

dat <- dat_back
dat <- t(dat)
dat <- as.data.frame(dat)
plate <- sample_ann$Biological_Condition # 這里定義分組信息
dat <-  cbind(dat, plate) #cbind根據(jù)列進行合并，即疊加所有列 #矩陣添加批次信息
dat[1:4, 1:4]
table(dat$plate)

# The variable plate (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[, -ncol(dat)], graph = FALSE)
head(dat.pca$var$coord) ## 每個主成分的基因重要性占比
head(dat.pca$ind$coord) ## 每個細(xì)胞的前5個主成分取值。
fviz_pca_ind(
      dat.pca,
      #repel =T,
      geom.ind = "point",
      # show points only (nbut not "text")
      col.ind = dat$plate,
      # color by groups
      #palette = c("#00AFBB", "#E7B800"),
      addEllipses = TRUE,
      # Concentration ellipses
      legend.title = "Groups"
)

同樣的，很明顯可以看到NPC跟另外的GW細(xì)胞群可以區(qū)分的很好，但是GW本身的3個小群體并沒有那么好的區(qū)分度。

接著是稍微高大上的tSNE降維

因為計算量的問題，這里先選取PCA后的主成分，然進行tSNE，當(dāng)然，也有其它做法，比如選取變化高的基因，顯著差異基因等等。

# 選取前面PCA分析的5個主成分。
dat_matrix <- dat.pca$ind$coord
# Set a seed if you want reproducible results
set.seed(42)
library(Rtsne) 
# 如果使用原始表達(dá)矩陣進行 tSNE耗時很可怕，dat_matrix = dat_back
# 出現(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(dat_matrix,perplexity=10)
plate <- sample_ann$Biological_Condition # 這里定義分組信息
plot(tsne_out$Y,col= rainbow(4)[as.numeric(as.factor(plate))], pch=19)

對PCA或者tSNE結(jié)果進行kmeans或者dbscan算法聚類

降維是降維，聚類是聚類，需要理解其中的區(qū)別。

降維與否，不同的降維算法選擇，不同參數(shù)的選擇得到的結(jié)果都不一樣。

聚類也是一樣，不同的算法，不同的參數(shù)。

# 前面我們的層次聚類是針對全部表達(dá)矩陣，這里我們?yōu)榱斯?jié)省計算量，可以使用tsne_out$Y這個結(jié)果
head(tsne_out$Y)
opt_tsne=tsne_out$Y
table(kmeans(opt_tsne,centers = 4)$clust)
plot(opt_tsne,  col=kmeans(opt_tsne,centers = 4)$clust, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
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)
# 比較兩個聚類算法區(qū)別
table(kmeans(opt_tsne,centers = 4)$clust,dbscan(opt_tsne,eps=3.1)$cluster)

3 M3Drop Package分析scRNA數(shù)據(jù)

首先構(gòu)建M3Drop需要的對象

library(M3Drop) 
Normalized_data <- M3DropCleanData(counts, 
                                   labels = sample_ann$Biological_Condition , 
                                   is.counts=TRUE, min_detected_genes=2000)

##查看對象
dim(Normalized_data$data)
length(Normalized_data$labels)
class(Normalized_data)
str(Normalized_data)

這個包設(shè)計比較簡單，并沒有構(gòu)建S4對象，只是一個簡單的list而已。

統(tǒng)計學(xué)算法 Michaelis-Menten

需要深入讀該文章，了解其算法，這里略過，總之它對單細(xì)胞轉(zhuǎn)錄組的表達(dá)矩陣進行了一系列的統(tǒng)計檢驗。

fits <- M3DropDropoutModels(Normalized_data$data)

# Sum absolute residuals
data.frame(MM=fits$MMFit$SAr, Logistic=fits$LogiFit$SAr,
           DoubleExpo=fits$ExpoFit$SAr) 
# Sum squared residuals
data.frame(MM=fits$MMFit$SSr, Logistic=fits$LogiFit$SSr,
           DoubleExpo=fits$ExpoFit$SSr)
QQ截圖20220901154758.png

找差異基因

DE_genes <- M3DropFeatureSelection(Normalized_data$data, 
                                   mt_method="fdr", mt_threshold=0.01)
dim(DE_genes)
head(DE_genes)

這里是針對上面的統(tǒng)計結(jié)果來的。


QQ截圖20220901155120.png

針對差異基因畫熱圖

par(mar=c(1,1,1,1)) 
heat_out <- M3DropExpressionHeatmap(DE_genes$Gene, Normalized_data$data, 
                                    cell_labels = Normalized_data$labels)

可視化了解一下找到的差異基因在不同的細(xì)胞類型的表達(dá)分布情況。


QQ截圖20220901155236.png

聚類

這里可以重新聚類后，針對自己找到的類別來分別找marker基因，不需要使用測試數(shù)據(jù)自帶的表型信息。

cell_populations <- M3DropGetHeatmapCellClusters(heat_out, k=4)
library("ROCR") 
marker_genes <- M3DropGetMarkers(Normalized_data$data, cell_populations)
table(cell_populations,Normalized_data$labels)

每個類別的marker genes

head(marker_genes[marker_genes$Group==4,],20) 
marker_genes[rownames(marker_genes)=="FOS",]

也可以針對這些 marker genes去畫熱圖，當(dāng)然，得根據(jù)AUC和P值來挑選符合要求的差異基因去繪圖。

par(mar=c(1,1,1,1)) 
choosed_marker_genes=as.character(unlist(lapply(split(marker_genes,marker_genes$Group), function(x) (rownames(head(x,20))))))
heat_out <- M3DropExpressionHeatmap(choosed_marker_genes, Normalized_data$data, cell_labels =  cell_populations)

如果遇到Error in plot.new() : figure margins too large報錯，則單獨將heat_out這行命令復(fù)制出來運行

對感興趣基因集進行注釋

通常是GO/KEGG等數(shù)據(jù)庫，通常是超幾何分布，GSEA，GSVA等算法。

這里就略過。

顯示運行環(huán)境

sessionInfo()

下一篇我們繼續(xù)用單細(xì)胞的R包來處理數(shù)據(jù)集。

加油！我們下一篇再見！

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容