比較5種scRNA鑑定HVGs方法

2021-02-15 單細胞天地

單細胞轉錄組雖然說不太可能跟傳統的bulk轉錄組那樣對每個樣本都測定到好幾萬基因的表達量,如果是10x這樣的技術,每個細胞也就幾百個基因是有表達量的,但是架不住細胞數量多,對大量細胞來說,這個表達矩陣的基因數量就很可觀了。一般來說,gencode資料庫的gtf文件注釋到的五萬多個基因裡面,包括蛋白編碼基因和非編碼的,可以在單細胞表達矩陣裡面過濾低表達量後,還剩下2萬多個左右。

2萬個基因,近萬個細胞的表達矩陣,降維起來是一個浩大的計算量,所以有一個步驟,就是從2萬個基因裡面挑選一下 highly variable genes (HVG) 進行下遊分析,從此就假裝自己的單細胞轉錄組就僅僅是測到了HVGs,數量嘛,兩千多吧!

那麼,問題就來了, 什麼樣的標準,算法來挑選 highly variable genes (HVG) 進行下遊分析呢?

搜索最新文獻

簡單谷歌搜索highly variable genes (HVG) 加上單細胞,這樣的關鍵詞就可以看到很多手把手教程:

2019六月的:Current best practices in single‐cell RNA‐seq analysis: a tutorial 地址:https://www.embopress.org/doi/10.15252/msb.20188746

bioconductor教程:Using scran to analyze single-cell RNA-seq data https://bioconductor.org/packages/devel/bioc/vignettes/scran/inst/doc/scran.html

一步步手把手流程:https://dockflow.org/workflow/simple-single-cell/

基本上經過四五個小時的模式,你就會總結到一些常見的挑選 highly variable genes (HVG) 的算法和R包,我簡單羅列5個,會對比一下:

seurat包的FindVariableFeatures函數,默認挑選2000個

monocle包的monocle_hvg函數

scran包的decomposeVar函數

M3Drop的BrenneckeGetVariableGenes

以及文獻參考::(Brennecke et al 2013 method) Extract genes with a squared coefficient of variation >2 times the fit regression

在scRNAseq數據集比較這5個方法

接下來我就沒有空繼續做這些瑣碎的比較啦,老規矩,跟我們的單細胞一期和二期學習視頻一下,學習筆記我們有專員操刀,下面看公司學習者的表演:

目的:利用scRNA-seq包的表達矩陣測試幾個R包尋找HVGs,然後畫upsetR圖看看不同方法的HVG的重合情況。

rm(list = ls())  
options(stringsAsFactors = F)

關於測試數據scRNAseq

包中內置了Pollen et al. 2014 的數據集(https://www.nature.com/articles/nbt.2967),到19年8月為止,已經有446引用量了。只不過原文完整的數據是 23730 features, 301 samples,這個包中只選取了4種細胞類型:pluripotent stem cells 分化而成的 neural progenitor cells (NPC,神經前體細胞) ,還有 GW16(radial glia,放射狀膠質細胞) 、GW21(newborn neuron,新生兒神經元) 、GW21+3(maturing neuron,成熟神經元)

library(scRNAseq)
#
data(fluidigm)

# 得到RSEM矩陣
assay(fluidigm)  <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]
#
#
#
#
#
# 樣本注釋信息
pheno_data <- as.data.frame(colData(fluidigm))
table(pheno_data$Coverage_Type)
#
#
#
table(pheno_data$Biological_Condition)
#
#
#

利用Seurat V3

suppressMessages(library(Seurat))
packageVersion("Seurat") 

seurat_sce <- CreateSeuratObject(counts = ct,
                          meta.data = pheno_data,
                          min.cells = 5,
                          min.features = 2000,
                          project = "seurat_sce")
seurat_sce <- NormalizeData(seurat_sce)

seurat_sce <- FindVariableFeatures(seurat_sce, selection.method = "vst", nfeatures=2000)
VariableFeaturePlot(seurat_sce)

image-20191127150351131

seurat_hvg <- VariableFeatures(seurat_sce)
length(seurat_hvg)
## [1] 2000
head(seurat_hvg)
## [1] "FOS"     "ERBB4"   "SCD"     "SGPL1"   "ARID5B"  "IGFBPL1"

利用Monocle

library(monocle)
gene_ann <- data.frame(
  gene_short_name = row.names(ct),
  row.names = row.names(ct)
)
sample_ann <- pheno_data
# 然後轉換為AnnotatedDataFrame對象
pd <- new("AnnotatedDataFrame",
          data=sample_ann)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)
monocle_cds <- newCellDataSet(
  ct,
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)
monocle_cds
#
#
#
#
#
#
#
#
#
#
#
#
#
#
monocle_cds <- estimateSizeFactors(monocle_cds)
monocle_cds <- estimateDispersions(monocle_cds)
disp_table <- dispersionTable(monocle_cds)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
monocle_cds <- setOrderingFilter(monocle_cds, unsup_clustering_genes$gene_id)
plot_ordering_genes(monocle_cds)

image-20191127150429288

monocle_hvg <- unsup_clustering_genes[order(unsup_clustering_genes$mean_expression, decreasing=TRUE),][,1]
length(monocle_hvg)
## [1] 13986
# 也取前2000個
monocle_hvg <- monocle_hvg[1:2000]
head(monocle_hvg)
## [1] "MALAT1" "TUBA1A" "MAP1B"  "SOX4"   "EEF1A1" "SOX11"

利用scran

library(SingleCellExperiment)
library(scran)
scran_sce <- SingleCellExperiment(list(counts=ct))
scran_sce <- computeSumFactors(scran_sce)
scran_sce <- normalize(scran_sce)
fit <- trendVar(scran_sce, parametric=TRUE,use.spikes=FALSE)
dec <- decomposeVar(scran_sce, fit)

plot(dec$mean, dec$total, xlab="Mean log-expression",
     ylab="Variance of log-expression", pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)

image-20191127150436995

scran_df <- dec
scran_df <- scran_df[order(scran_df$bio, decreasing=TRUE),]
scran_hvg <- rownames(scran_df)[scran_df$FDR<0.1]
length(scran_hvg)
## [1] 875
head(scran_hvg)
## [1] "IGFBPL1" "STMN2"   "ANP32E"  "EGR1"    "DCX"     "XIST"

利用M3Drop

library(M3DExampleData)
# 需要提供表達矩陣(expr_mat)=》normalized or raw (not log-transformed)
HVG <-M3Drop::BrenneckeGetVariableGenes(expr_mat=ct, spikes=NA, suppress.plot=FALSE, fdr=0.1, minBiolDisp=0.5, fitMeanQuantile=0.8)

image-20191127150443846

M3Drop_hvg <- rownames(HVG)
length(M3Drop_hvg)
## [1] 917
head(M3Drop_hvg)
## [1] "CCL2"   "EMP1"   "ZNF630" "NRDE2"  "PLCL1"  "KCTD21"

自定義函數

來自:(Brennecke et al 2013 method) Extract genes with a squared coefficient of variation >2 times the fit regression

實現了:Select the highly variable genes based on the squared coefficient of variation and the mean gene expression and return the RPKM matrix the the HVG

if(T){
  getMostVarGenes <- function(
    data=data,              
    fitThr=1.5,             
    minMeanForFit=1         
  ){
    
    
    data_no0 <- as.matrix(
      data[rowSums(data)>0,]
    )
    
    meanGeneExp <- rowMeans(data_no0)
    names(meanGeneExp)<- rownames(data_no0)

    
    varGenes <- rowVars(data_no0)
    cv2 <- varGenes / meanGeneExp^2

    
    useForFit <- meanGeneExp >= minMeanForFit

    
    fit <- glmgam.fit( cbind( a0 = 1,
                              a1tilde = 1/meanGeneExp[useForFit] ),
                       cv2[useForFit] )
    a0 <- unname( fit$coefficients["a0"] )
    a1 <- unname( fit$coefficients["a1tilde"])

    
    fit_genes <- names(meanGeneExp[useForFit])
    cv2_fit_genes <- cv2[useForFit]
    fitModel <- fit$fitted.values
    names(fitModel) <- fit_genes
    HVGenes <- fitModel[cv2_fit_genes>fitModel*fitThr]
    print(length(HVGenes))

    
    plot_meanGeneExp <- log10(meanGeneExp)
    plot_cv2 <- log10(cv2)
    plotData <-  data.frame(
      x=plot_meanGeneExp[useForFit],
      y=plot_cv2[useForFit],
      fit=log10(fit$fitted.values),
      HVGenes=log10((fit$fitted.values*fitThr))
    )
    p <- ggplot(plotData, aes(x,y)) +
      geom_point(size=0.1) +
      geom_line(aes(y=fit), color="red") +
      geom_line(aes(y=HVGenes), color="blue") +
      theme_bw() +
      labs(x = "Mean expression (log10)", y="CV2 (log10)")+
      ggtitle(paste(length(HVGenes), " selected genes", sep="")) +
      theme(
        axis.text=element_text(size=16),
        axis.title=element_text(size=16),
        legend.text = element_text(size =16),
        legend.title = element_text(size =16 ,face="bold"),
        legend.position= "none",
        plot.title = element_text(size=18, face="bold", hjust = 0.5),
        aspect.ratio=1
      )+
      scale_color_manual(
        values=c("#595959","#5a9ca9")
      )
    print(p)

    
    HVG <- data_no0[rownames(data_no0) %in% names(HVGenes),]
    return(HVG)
  }
}
library(statmod)
diy_hvg <- rownames(getMostVarGenes(ct))

看起來代碼比較長,主要是繪圖的時候拖拉了。

image-20191127150504326

length(diy_hvg)
## [1] 2567
head(diy_hvg)
## [1] "A2M"     "A2MP1"   "AAMP"    "ABCA11P" "ABCA3"   "ABCB10"

upsetR作圖

require(UpSetR)
input <- fromList(list(seurat=seurat_hvg, monocle=monocle_hvg,
                       scran=scran_hvg, M3Drop=M3Drop_hvg, diy=diy_hvg))
upset(input)

很多圖都是可以美化的,不過並不是我們的重點

image-20191127150523192

可以看到,不同的R包和方法,差異不是一般的大啊!!!

更多方法

比如基尼係數:findHVG: Finding highly variable genes (HVG) based on Gini 見:https://rdrr.io/github/jingshuw/descend/man/findHVG.html

趕快試一下吧,探索的過程中,你就會對單細胞數據處理的認知加強了。

什麼才是你應該學的單細胞

如果是10X儀器的單細胞轉錄組數據走cellranger流程,我們在單細胞天地多次分享過流程筆記:

如果是smart-seq2技術,首先走單細胞下遊分析標準流程啊,就是那些R包的認知,包括 scater,monocle,Seurat,scran,M3Drop 需要熟練掌握它們的對象,:一些單細胞轉錄組R包的對象 ,分析流程也大同小異:

step1: 創建對象

step2: 質量控制

step3: 表達量的標準化和歸一化

step4: 去除幹擾因素(多個樣本整合)

step5: 判斷重要的基因

step6: 多種降維算法

step7: 可視化降維結果

step8: 多種聚類算法

step9: 聚類後找每個細胞亞群的標誌基因

step10: 繼續分類

最後說幾句話

相關焦點

  • 胎兒鑑定性別方法 這5種方法準確又科學
    如今網絡上流傳的鑑別胎兒性別的方法層出不窮,數不勝數,很多父母又不敢輕易嘗試,害怕會傷害到孩子,那麼,除了網上的鑑別胎兒性別的方法之外,還有沒有其他的簡單有效的鑑別胎兒性別的方法?今天我們就來看看胎兒鑑定性別的方法。
  • 種與品種——物種鑑定
    ):通常指人工選育出的具有一定優良農藝、園藝性狀的種下單位,一個品種對應一個實體;如:   桐鄉菊花栽培品種——菊科植物菊花Dendranthema morifolium (Ramat.)Tzvel.渭黨一號   拉丁學名—雙名法:《植物命名法規》規定每一個物種均只對應一個拉丁正名,即:種名=屬名+種加詞+定名人(可縮寫)。 黨參:Codonopsispilosula(Franch.)Nannf   菊花:Dendranthema morifolium(Ramat.) Tzvel.
  • 比較二次根式大小的8種方法,第5種最常用!
    別擔心,本節小隴整理的8種比較大小的方法,如果你能全掌握,那就可以對比較大小的題目「通吃」了,這8種方法不僅適用於二次根式大小的比較,對於其他數的大小比較也適用,當然,本節是結合二次根式比較大小的題型來講述這8種方法,既學會了二次根式大小的比較,又掌握了8種比較大小的方法,可謂收穫良多。
  • 和田玉真假的鑑定方法,如何鑑定和田玉的真假?
    正因為和田玉的市場需求,市場上也充斥著不少假和田,而接下來,玉雕名家要在這裡給大家介紹幾種和田玉鑑定方法,希望喜愛玉石的朋友都能夠結緣到適合自己的天然玉石。和田玉真假的鑑定方法--燈照法可以用強光手電或者手機閃光燈側面照和田玉,如果發現玉石內部有空心的小氣泡,則合成玉石;如果看到棉絮狀結構,則為和田玉。如果看到顆粒裝結構則為東陵石或石英石。
  • 死活細胞的鑑定方法
    活細胞的鑑定在生物學和醫學上具有很重要的意義。細胞培養過程中要隨時記錄細胞的生長情況,需要經常測定細胞的存活率;在腫瘤細胞的研究中,為了檢驗各種藥物對腫瘤細胞的殺傷力,也需要測定腫瘤細胞的存活率。在臨床醫學中死活細胞的鑑定也有很大的應用,例如為了檢測某一男子的生育能力,測定精子細胞的存活力是比較常用的辦法。
  • 鑽石鑑定知識科普 鑽戒鑑定方法有哪幾種
    鑽戒之於大多數消費者來說,屬於價格比較昂貴的物品,每個人都希望物有所值,甚至物超所值。然而珠寶市場充斥著以次充好、以假亂真的諸多亂象,給消費者選購鑽戒造成了很大的困擾。因此,適當掌握一些鑽石鑑定知識成為了一項婚前必備技能。接下來,我們一起學習幾個鑽戒鑑定的小方法吧。
  • 5種方法比較分數大小,將分母進行通分只是其中最普通的一種
    上一篇文章我們說了整數以及小數,怎麼比較大小?今天我們繼續小學階段的分數如何比較大小?簡單介紹5種常用的方法。1.分母通分。最常見的分母通分比如說3/4和4/5,將這兩個分數分母進行通分,分母的最小公倍數是20,進行通分後,3/4=15/20,4/5=16/20,所以4/5大於3/4。2.分子「通分」。
  • 科學家提出鑑定罕見病基因新方法
    科學家提出鑑定罕見病基因新方法 作者:小柯機器人 發布時間:2019/7/9 13:30:34 史丹福大學Stephen B.Montgomery研究團隊的一項最新研究,提出了利用血液轉錄組測序和大型對照組鑑定罕見病基因的方法。 2019年6月出版的《Nature Medicine》發表了這項成果。 該課題組試圖評估血液中的RNA-seq作為診斷不同病理生理學罕見疾病工具的效用。研究人員從94名患有16種不同疾病類別的未確診罕見疾病的患者身上提取了全血進行RNA-seq。
  • Method | 利用Proximity-CLIP的方法鑑定亞細胞水平的RNA及其結合蛋白
    本文基於此,發展了新的proximity-CLIP方法,將細胞內區域特異性的RNA結合蛋白(RBP)生物素化後利用UV交聯RNA,作者認為這樣的方法是可以鑑定到局部轉錄產物的。可以用此來檢測局部RBPome、定量局部轉錄組、鑑定RBP結合部分的RNA等。作者將此方法用於研究細胞核、細胞質、細胞間界面的RNA。
  • 鑑定青銅器,科學得方法原來是這樣
    如果您恰巧有一件青銅器,小編給大家介紹一下科學鑑定的主要方法,與其相比,人工鑑定都是浮雲。第一種方法:金相分析如果您是材料專業畢業的,應該知道什麼是金屬的金相分析,通過測金屬的相圖,能夠掌握合金中幾種金屬的成分及分散方式,現代仿品外表可以做得很像,結構還會那麼像嗎?這種方法有一個缺點,就是過程需要時間較長,而且費用較高。
  • 論文|3種鶇科鳥類的性別鑑定研究
    自1995年CHD基因被應用於鳥類性別鑑定以來,CHD基因已成為突胸鳥類性別鑑定最重要的分子標記,也是目前鳥類性別鑑定應用最廣泛的方法之一。但實踐中發現該方法在很多情況下難以鑑定或鑑定準確性難以保證。目 的總結現有文獻中3種鶇科鳥類性別鑑定的形態特徵,討論3種鶇科鳥類形態學方法和微衛星性別鑑定的準確率和可行性,旨在為這些鳥類的野外研究和種群性別結構分析提供借鑑。
  • 【珠寶知識】 教你肉眼鑑定11種常見寶石
    原標題:【珠寶知識】 教你肉眼鑑定11種常見寶石 珠寶贗品主要是在材料方面作假,以偽劣材料冒充真品,以躋身於高檔珠寶之列。下面來教你怎麼用肉眼去區分它們。
  • 一種快速鑑定RNA 質量的方法
    方法 在常規瓊脂糖凝膠電泳的基礎上,結合甲醛變性瓊脂糖凝膠電泳加以改進而成。結果 提取的總RNA 用此方法電泳後可得到28S、18S、5. 8S(5S) 三條清晰的rRNA 條帶。結論 此方法可以達到對RNA 完整性進行快速鑑定的目的。 隨著分子生物學技術的迅猛發展及在醫學領域中的普及,以RNA 為研究對象的實驗與日俱增。RNA 質量的好壞直接關係到後續實驗的成敗。
  • Biochemistry ▏穩定α/β-水解酶的五種蛋白質工程策略的比較
    蛋白質工程策略的複雜性和先決條件知識各不相同,但很少有對不同策略進行並排比較的例子。該研究旨在通過比較用於穩定蛋白質的各種策略來彌補文獻中的這種疏忽。首先是基於結構的理性設計和半理性設計。其次,是基於序列的的設計,一種方法是將靶蛋白中的胺基酸替換為更穩定的同源物中的胺基酸。另一種方法是用同源物中保守的胺基酸替換靶蛋白中的胺基酸,而不管同源物是否更穩定(共識突變)。
  • 13種單細胞RNA擴增測序方法的比較
    然而,各種scRNA-seq技術之間存在著重要的差異,目前並不清楚哪種方法最適合繪製組織、器官和生物體的細胞圖譜。Nature biotechnology最近發表的一篇文章通過多個實驗室合作,對同樣的細胞樣本進行了13種方法的單細胞測序,系統地評估了這幾種技術在全面描述細胞類型和狀態方面的能力。
  • 蔣玉晨:明清古玉的痕跡學與鑑定方法
    今天,我們請到了鬼市平臺玉器、錢幣、銅器鑑寶師蔣玉晨老師為大家講解明清古玉的鑑定方法。 01玉器收藏的趨勢:精品導向型 在我國,古玩收藏的歷史可以追溯到宋代。宋代盛行金石學。這裡的「金」就是古青銅器,「石」指石刻碑碣。這個定義還可以延伸到竹、骨、玉、石等常見的文玩藏品。
  • 黃金密度法鑑定及貴金屬鑑定方法
    密度法(比重法):密度法又稱靜水力學法,根據黃金標準密度而估算出被測黃金飾品的成色含量,它是根據物理學裡的阿基米德定律而來,即物體的浮力等於排開這種液體的重量來測算的
  • 張承青電鏡實驗室環境約稿[5]:幾種改善電磁環境方法比較
    (本文經授權發布,分享內容為作者個人觀點, 僅供讀者學習參考,不代表本網觀點)系列之五 幾種改善電磁環境方法比較被動式低頻電磁屏蔽根據屏蔽材料不同主要分為兩種:一種是使用高導磁材料(如鋼、矽鋼、玻莫合金等),另一種是使用高導電材料(
  • 「蛋白晶片法or質譜法」蛋白質組學研究方法的選擇及比較
    經常會有人將兩種方法進行比較,科學上很難有一勞永逸的方法,每一種方法都有其優勢和局限,有各自的用武之地。那麼,研究蛋白質組學應該如何選擇更合適的研究方法呢?以下將從六個主要方面進行比較和推薦。近年來,高通量技術平臺的開發大大推動了蛋白質組學研究工作。這些適用於蛋白質分析和標誌物開發的平臺讓研究人員能夠在單次分析中檢測數千種蛋白質。
  • 質量鑑定工作中常用的物相分析方法(物理法)
    物相分析(phase analysis)是指用化學或物理方法測定材料礦物組成及其存在裝填的分析方法。材料的性質一般取決於其化學礦物的組成和結構狀況。一個複雜產品也是由各種基本材料所構成的,在鑑定一個產品的失效原因時,其往往取決於其中的物相組成、分布及各相的特性,包括礦物種類、數量、晶型、晶粒大小、分布狀況、結合方式、形成固溶體及玻璃相等。天縱鑑定(SKYLABS)在日常產品質量鑑定工作中,常用的物相分析方法總結起來有X射線衍射分析、雷射拉曼分析、傅立葉紅外分析以及微區電子衍射分析等。