Filter genes and low-quality single cells

Filter genes

We filter the genes to exlude both those that are lowly expressed or over-expressed (>= 1024 molecules in a given cell).

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

molecules_cpm_mean <- rowMeans(cpm(molecules_filter, log = TRUE))
hist(molecules_cpm_mean, xlab = "Mean log2 molecule cpm in single cells",
     ylab = "Number of genes", main = "Identifying expression cutoff")
lower_exp_cutoff <- 2
abline(v = lower_exp_cutoff, col = "red")

genes_pass_filter <- rownames(molecules_filter)[molecules_cpm_mean > lower_exp_cutoff]

13112 genes have a mean log₂ molecule cpm greater than 2, including 48 ERCC genes.

Next we identify any genes which have greater than 1024 molecules in any given single cell. These are above our theoretical maximum number of UMIs (it can happen when a highly expressed gene as multiple start sites), and thus we cannot correct them for the collision probability.

overexpressed_rows <- apply(molecules_filter, 1, function(x) any(x >= 1024))
overexpressed_genes <- rownames(molecules_filter)[overexpressed_rows]
overexpressed_genes

[1] "ENSG00000198888" "ENSG00000198804" "ENSG00000198712" "ENSG00000198899"
[5] "ENSG00000198938" "ENSG00000198886"

ensembl <- useMart(host = "grch37.ensembl.org",
                   biomart = "ENSEMBL_MART_ENSEMBL",
                   dataset = "hsapiens_gene_ensembl")
overexpressed_genes_info <- getBM(attributes = c("ensembl_gene_id", "chromosome_name",
                                                 "external_gene_name", "transcript_count",
                                                 "description"),
                                  filters = "ensembl_gene_id",
                                  values = overexpressed_genes[grep("ENSG", overexpressed_genes)],
                                  mart = ensembl)
overexpressed_genes_info

  ensembl_gene_id chromosome_name external_gene_name transcript_count
1 ENSG00000198712              MT             MT-CO2                1
2 ENSG00000198804              MT             MT-CO1                1
3 ENSG00000198886              MT             MT-ND4                1
4 ENSG00000198888              MT             MT-ND1                1
5 ENSG00000198899              MT            MT-ATP6                1
6 ENSG00000198938              MT             MT-CO3                1
                                                                     description
1  mitochondrially encoded cytochrome c oxidase II [Source:HGNC Symbol;Acc:7421]
2   mitochondrially encoded cytochrome c oxidase I [Source:HGNC Symbol;Acc:7419]
3     mitochondrially encoded NADH dehydrogenase 4 [Source:HGNC Symbol;Acc:7459]
4     mitochondrially encoded NADH dehydrogenase 1 [Source:HGNC Symbol;Acc:7455]
5           mitochondrially encoded ATP synthase 6 [Source:HGNC Symbol;Acc:7414]
6 mitochondrially encoded cytochrome c oxidase III [Source:HGNC Symbol;Acc:7422]

6 genes have molecule counts greater than or equal to 1024 in at least one single cell, which includes 0 ERCC control genes.

Update the list of genes passing the filters.

genes_pass_filter <- setdiff(genes_pass_filter, overexpressed_genes)

Filter the data to only include the subset of 13106 genes which pass the lower and upper expression cutoffs. This subset includes 48 ERCC genes.

reads_filter <- reads_filter[rownames(reads_filter) %in% genes_pass_filter, ]
molecules_filter <- molecules_filter[rownames(molecules_filter) %in% genes_pass_filter, ]
reads_bulk_filter <- reads_bulk[rownames(reads_bulk) %in% genes_pass_filter, ]
stopifnot(nrow(reads_filter) == length(genes_pass_filter),
          dim(reads_filter) == dim(molecules_filter),
          nrow(reads_bulk_filter) == nrow(molecules_filter))

Output filtered data.

write.table(anno_filter, "../data/annotation-filter.txt", quote = FALSE,
            sep = "\t", row.names = FALSE)
write.table(reads_filter, "../data/reads-filter.txt", quote = FALSE,
            sep = "\t", col.names = NA)
write.table(molecules_filter, "../data/molecules-filter.txt", quote = FALSE,
            sep = "\t", col.names = NA)
write.table(reads_bulk_filter, "../data/reads-bulk-filter.txt", quote = FALSE,
            sep = "\t", col.names = NA)

Output list of genes that passed filters.

write.table(genes_pass_filter, "../data/genes-pass-filter.txt", quote = FALSE,
            sep = "\t", row.names = FALSE, col.names = FALSE)

PCA of filtered data

pca_reads_filter <- run_pca(reads_filter)
pca_reads_filter_plot <- plot_pca(pca_reads_filter$PCs, explained = pca_reads_filter$explained,
         metadata = anno_filter, color = "individual",
         shape = "replicate") +
  labs(title = "Filtered raw reads for single cells")
pca_reads_filter_plot

pca_molecules_filter <- run_pca(molecules_filter)
pca_molecules_filter_plot <- plot_pca(pca_molecules_filter$PCs, explained = pca_molecules_filter$explained,
         metadata = anno_filter, color = "individual",
         shape = "replicate") +
  labs(title = "Filtered raw molecules for single cells")
pca_molecules_filter_plot

Because we start with the union of observed genes in the single cell or bulk samples and then the expression cutoff is based on the expresssion in the single cells, it is possible that a gene is not observed at all in the bulk samples. This causes an error with the PCA because that gene is invariant, so they are filtered before performing PCA.

reads_bulk_zeros <- rowSums(reads_bulk_filter) == 0
rownames(reads_bulk_filter)[reads_bulk_zeros]

[1] "ERCC-00048"

pca_reads_bulk_filter <- run_pca(reads_bulk_filter[!reads_bulk_zeros, ])
pca_reads_bulk_filter_plot <- plot_pca(pca_reads_bulk_filter$PCs, explained = pca_reads_bulk_filter$explained,
         metadata = anno_bulk, color = "individual",
         shape = "replicate") +
  labs(title = "Filtered raw reads for bulk samples")
pca_reads_bulk_filter_plot