Filter genes

Low QC sample

We performed quality control to identify low quality single cells. Only keep high-quality single cells.

quality <- read.table("../data/quality-single-cells.txt", stringsAsFactors = FALSE)
colnames(quality) <- c("sample", "quality")
eset <- eset[, quality$quality]
dim(eset)

Features  Samples 
   54616     5597 

Over-expressed gene

We would like to remove genes which have greater than 4096 (44444*4) molecules in any given single cell. These are above our theoretical maximum number of 6-bp UMIs (it can happen when a highly expressed gene as multiple start sites). We did not find any genes with molecules greater than 4096. We found r length(overexpressed_genes_1024) mitochondrial genes which have molecules larger than 1024

## look at human genes
eset_hs <- eset[fData(eset)$source == "H. sapiens", ]
dim(eset_hs)

Features  Samples 
   20151     5597 

head(featureNames(eset_hs))

[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
[5] "ENSG00000000460" "ENSG00000000938"

## identify gene with molecule > 4096
overexpressed_rows <- apply(exprs(eset_hs), 1, function(x) any(x >= 4096))
overexpressed_genes <- rownames(exprs(eset_hs))[overexpressed_rows]
overexpressed_genes

character(0)

## identify gene with molecule > 1024
overexpressed_rows_1024 <- apply(exprs(eset_hs), 1, function(x) any(x >= 1024))
overexpressed_genes_1024 <- rownames(exprs(eset_hs))[overexpressed_rows]
overexpressed_genes_1024

character(0)

Lowly expressed

We identify the lower cutoff using the mean log2 molecule counts per million (cpm) in the high quality single cells.

## convert to log2 cpm
mol_hs_cpm <- cpm(exprs(eset_hs), log = TRUE)
mol_hs_cpm_means <- rowMeans(mol_hs_cpm)
summary(mol_hs_cpm_means)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.306   2.354   2.974   3.746   4.657  12.940 

hist(mol_hs_cpm_means, xlab = "Mean log2 molecule cpm in single cell samples",
     ylab = "Number of genes", main = "Identifying expression cutoff")

## set a cutoff
lower_exp_cutoff <- 3
abline(v = lower_exp_cutoff, col = "red")

genes_pass_filter <- rownames(mol_hs_cpm)[mol_hs_cpm_means > lower_exp_cutoff]

9957 genes have a mean log2 molecule cpm greater than 3

mol_hs_cpm <- data.frame(mol_hs_cpm)
mol_hs_cpm$mean <- apply(mol_hs_cpm, 1, mean)
mol_hs_cpm$filter <- mol_hs_cpm$mean > lower_exp_cutoff
write.table(data.frame(row.names(mol_hs_cpm), mol_hs_cpm[,"filter"]),
            file = "../data/genes-pass-filter.txt", quote = FALSE,
            sep = "\t", row.names = FALSE, col.names = FALSE)