比較不同單細胞轉錄組資料尋找features方法
挑選到的跟feature相關的基因集，有點類似於在某些組間差異表達的基因集，都需要後續功能註釋。
背景介紹
單細胞轉錄組測序的確可以一次性對所有細胞都檢測到上千個基因的表達，但是，大多數情況下，只有其中的少部分基因是有生物學意義的，比如可以區分不同的細胞型別，或者分化發育相關的基因，或者細胞應對外界刺激的。
而且大多數基因之所以在不同的細胞裡面表達有差異，其實是技術限制，背景噪音。這些技術限制，包括批次效應，都會阻礙我們發現那些真正的有生物學意義的基因。
所以做 feature selection 分析來去除那些技術噪音相關基因，可以顯著的提高信噪比，降低後續分析的複雜度。
包的安裝
所以需要安裝並且載入一些包,安裝程式碼如下；
install.packages('ROCR')
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("M3Drop")
install.packages("devtools")
library("devtools")
install_github("BPSC","nghiavtr")
install_github("hemberg-lab/scRNA.seq.funcs") 載入程式碼如下：
library(scRNA.seq.funcs) library(matrixStats) library(M3Drop) library(RColorBrewer) set.seed(1)
載入測試資料
這裡選取的是Usoskin et al 文章單細胞測序文章的資料，包含4種細胞型別：
- NP = non-peptidergic nociceptors
- PEP = peptidergic nociceptors
- NF = neurofilament containing
- TH = tyrosine hydroxylase containing neurons.
對應著25334個基因在 622 個細胞裡面的表達量
# http://www.biotrainee.com/jmzeng/scRNA/usoskin1.rds usoskin1 <- readRDS("../usoskin/usoskin1.rds") dim(usoskin1)
## [1] 25334 622table(colnames(usoskin1))##
## NF NP PEP TH
## 139 169 81 233用M3Drop對錶達矩陣進行一站式過濾
uso_list <- M3Drop::M3DropCleanData(
usoskin1,
labels = colnames(usoskin1),
min_detected_genes = 2000,
is.counts = TRUE
)
expr_matrix <- uso_list$data # Normalized & filtered expression matrix
dim(expr_matrix)## [1] 15708 532celltype_labs <- uso_list$labels # filtered cell-type labels
cell_colors <- brewer.pal(max(3,length(unique(celltype_labs))), "Set3")這個M3Drop的M3DropCleanData函式會自動過濾那些表達基因數量很少的細胞，過濾低表達基因，然後把reads counts的表達矩陣轉換為 counts per million (CPM) 。
過濾之後，只剩下 15708個基因在532個細胞 的表達了。
尋找highly variable genes (HVG)
那些在樣本群體裡面表達量變異比較大的基因可能是真正的生物學現象，也有可能是技術誤差，而且變異程度總是跟基因的表達量成正相關。如下圖所示：
plot(rowMeans(expr_matrix),rowVars(expr_matrix),log='xy')
Brennecke et al.提出了演算法來矯正這一影響，這個方法也被包裝成了Brennecke_getVariableGenes(counts, spikes) 函式，但是這個資料並沒有ERCC spike-in，所以直接對整個表達矩陣處理即可。
Brennecke_HVG <- M3Drop::BrenneckeGetVariableGenes(
expr_matrix,
fdr = 0.01,
minBiolDisp = 0.5
)
探究 High Dropout Genes
另外一個尋找HVGs是檢視它們是否有非常多的0表達量情況，這種多0表達的情況叫做dropout rate，通常單細胞轉錄組表達矩陣裡面過半數的基因都是0表達的。因為單細胞裡面的mRNA很多無法被反轉錄，這種情況可以用Michaelis-Menten等式來模擬，如下圖所示：
K = 49
S_sim = 10^seq(from=-3, to=4, by=0.05) # range of expression values
MM = 1-S_sim/(K+S_sim)
plot(S_sim, MM, type="l", lwd=3, xlab="Expression", ylab="Dropout Rate", xlim=c(1,1000))
S1 = 10; P1 = 1-S1/(K+S1) # Expression & dropouts for cells in condition 1
S2 = 750; P2 = 1-S2/(K+S2) # Expression & dropouts for cells in condition 2
points(c(S1,S2),c(P1,P2), pch=16, col="grey85", cex=3)
mix = 0.5; # proportion of cells in condition 1
points(S1*mix+S2*(1-mix), P1*mix+P2*(1-mix), pch=16, col="grey35", cex=3)
用來M3Drop包的M3DropFeatureSelection函式來挑選那些顯著偏離了Michaelis-Menten曲線的基因，這裡的閾值取1% FDR.
但是這個函式M3DropFeatureSelection依賴於正確的M3Drop包版本，下面就不運行了。
M3Drop_genes <- M3Drop::M3DropFeatureSelection(
expr_matrix,
mt_method = "fdr",
mt_threshold = 0.01
)
title(main = "Usoskin")
M3Drop_genes <- M3Drop_genes$GeneDepth-Adjusted Negative Binomial (DANB)
下面這個 NBumiConvertToInteger 也依賴於正確的M3Drop包版本，下面就不運行了。
usoskin_int <- NBumiConvertToInteger(usoskin1)
DANB_fit <- NBumiFitModel(usoskin_int) # DANB is fit to the raw count matrix
# Perform DANB feature selection
DropFS <- NBumiFeatureSelectionCombinedDrop(DANB_fit)
DANB_genes <- names(DropFS[1:1500])基因表達相關性
這個表達矩陣很大，所以計算所有基因之間的相關性耗時很長，為節約時間，也不運行了。
cor_mat <- cor(t(expr_matrix), method="spearman") #Gene-gene correlations
diag(cor_mat) <- rep(0, times=nrow(expr_matrix))
score <- apply(cor_mat, 1, function(x) {max(abs(x))}) #Correlation of highest magnitude
names(score) <- rownames(expr_matrix);
score <- score[order(-score)]
Cor_genes = names(score[1:1500])PCA挑選
PCA 速度還行，挑選第1，2主成分的前1500個基因
pca <- prcomp(log(expr_matrix+1)/log(2));
# PCA is typically performed on log-transformed expression data
plot(pca$rotation[,1], pca$rotation[,2], pch=16, col=cell_colors[as.factor(celltype_labs)]) # plot projection
score <- rowSums(abs(pca$x[,c(1,2)]))
# calculate loadings for components 1 and 2
names(score) <- rownames(expr_matrix)
score <- score[order(-score)]
PCA_genes = names(score[1:1500])檢查挑選的基因集的效果
熱圖+聚類可以看看基因是否在各個細胞型別差異表達，並且把細胞型別比較好的分開。
這個熱圖非常耗時，如無必要，請不要執行這個程式碼
M3Drop::M3DropExpressionHeatmap(
PCA_genes , ## 或者 M3Drop_genes 等其它方法挑到的基因
uso_list$data,
cell_labels = uso_list$labels
)挑選的基因集跟DEseq得到的差異基因列表看交集
載入用DEseq得到的差異基因列，跟前面得到的M3Drop_genes比較一下。
# Load DE genes
# http://www.biotrainee.com/jmzeng/scRNA/DESeq_table.rds
DESeq_table <- readRDS("../usoskin/DESeq_table.rds")
DE_genes = unique(DESeq_table$Gene)
# Calculate precision
# sum(M3Drop_genes %in% DE_genes)/length(M3Drop_genes)
sum(PCA_genes %in% DE_genes)/length(PCA_genes)
## [1] 0.382文中提到的測試資料：
http://www.biotrainee.com/jmzeng/scRNA/DESeq_table.rds
http://www.biotrainee.com/jmzeng/scRNA/pollen.rds
http://www.biotrainee.com/jmzeng/scRNA/tung_umi.rds
http://www.biotrainee.com/jmzeng/scRNA/usoskin1.rds