feature selection

比较不同feature selection方法

单细胞转录组测序的确可以一次性对所有细胞都检测到上千个基因的表达，但是大多数情况下，只有其中的少部分基因是有生物学意义。而且大多数基因之所以在不同的细胞里面表达有差异，其实是技术限制，背景噪音。这些技术限制，包括批次效应，都会阻碍我们发现那些真正的有生物学意义的基因。

feature selection分析去除那些技术噪音相关基因，可以显著的提高信噪比，降低后续分析的复杂度。

library(limma)
library(M3Drop)

## Warning: package 'M3Drop' was built under R version 3.5.2

library(scran)
usoskin1 <- readRDS("D:/paopaoR/testt/usoskin1.rds")
dim(usoskin1)

## [1] 25334   622

Normalize & filter expression matrix

uso_list <- M3Drop::M3DropCleanData(
    usoskin1,
    labels = colnames(usoskin1),
    min_detected_genes = 2000,
    is.counts = TRUE
)
expr_matrix <- uso_list$data 
dim(expr_matrix)

## [1] 15708   532

Brennecke

Brennecke_HVG <- M3Drop::BrenneckeGetVariableGenes(
    expr_matrix,
    fdr = 0.01,
    minBiolDisp = 0.5
)

image.png

High Dropout Genes

fit 米氏方程（按教程应该是用m3dropfeatureselection那个函数，我当时安装的版本可能有问题，一直没有这个函数，我就想既然原理是假设dropout是因为逆转录失败，这个酶促反应，那就挑偏离米氏方程model的吧，但是后来我删除了原来的包从github重新加载后，函数可以用了，心累。。。代码也放着吧）

m3drop <- M3DropDropoutModels(expr_matrix)

image.png

m3drop <- sort(m3drop$MMFit$predictions,decreasing = FALSE)
m3drop_gene <- names(m3drop)[1:1500]

m3dGenes <- as.character(
    M3DropFeatureSelection(deng)$Gene
)

fit model

拟合线性趋势

fit <- trendVar(expr_matrix, loess.args=list(span=0.05))
dec <- decomposeVar(expr_matrix, fit)
top.dec <- dec[order(dec$bio, decreasing=TRUE),]
top_genes <- rownames(top.dec)[1:1500]

correlation genes

这个很费时间，还很占内存。

pca first 2 components

用前两个pc

pca <- prcomp(log(expr_matrix+1)/log(2))
score <- rowSums(abs(pca$x[,c(1,2)]))
names(score) <- rownames(expr_matrix)
score <- score[order(-score)]
PCA_genes <- names(score[1:1500])

看看几种方法选择的基因overlap情况

venn.diag <- vennCounts(
 cbind(rownames(expr_matrix) %in% Brennecke_HVG,
 rownames(expr_matrix) %in% PCA_genes,rownames(expr_matrix) %in% m3drop_gene,rownames(expr_matrix) %in% top_genes))
vennDiagram(
 venn.diag,
 names = c("bre", "pca","m3d","fit"),
 circle.col = c("blue", "green","yellow","red"))

image.png

DEseq 看交集

DESeq_table <- readRDS("D:/paopaoR/testt/DESeq_table.rds")
DE_genes = unique(DESeq_table$Gene)

sum(Brennecke_HVG %in% DE_genes)/length(Brennecke_HVG)

## [1] 0.2275862

sum(PCA_genes %in% DE_genes)/length(PCA_genes)

## [1] 0.382

sum(m3drop_gene %in% DE_genes)/length(m3drop_gene)

## [1] 0.374

sum(top_genes %in% DE_genes)/length(top_genes)

## [1] 0.3646667