Patch-seq QC assessment for manuscript

Use the package patchSeqQC to assign QC scores (and in particular a contamination score) for each cell. To do this, we load the patchSeqQC library and then follow the relevant aspects of the tutorial here: https://github.com/PavlidisLab/patchSeqQC. The basic idea is as follows: (1) define markers for different broad classes, (2) determine expression of these genes in a reference data set, and (3) see how these genes are expressed in patch-seq data to assess contamination. The advantage of this method is it can be done independent of cluster calls, so long as you have a reliable set of marker genes for different broad classes. This requires both FACS and Patch-seq data to run.

It is important to note that this analysis is based almost entirely on the the Shreejoy Trapathy’s patchSeqQC library linked above, and which has now been published here: https://www.frontiersin.org/articles/10.3389/fnmol.2018.00363/full. Please be sure to cite this paper!

First load the required libraries and set the working directory. (NOTE: Probably most of these aren’t actually needed)

suppressPackageStartupMessages({
  library(feather)
  library(dplyr)
  library(bindrcpp)
  library(viridisLite)
  library(ggplot2)
  library(mfishtools)     # devtools::install_github("AllenInstitute/mfishtools")
  library(patchseqtools)  # devtools::install_github("AllenInstitute/patchseqtools")
  library(patchSeqQC)     # devtools::install_github('PavlidisLab/patchSeqQC')
  library(scrattch.io)    # devtools::install_github("AllenInstitute/scrattch.io")
})

package 㤼㸱dplyr㤼㸲 was built under R version 3.6.3package 㤼㸱ggplot2㤼㸲 was built under R version 3.6.3

options(stringsAsFactors = FALSE)

patchDir = "\\\\allen/programs/celltypes/workgroups/rnaseqanalysis/shiny/patch_seq/star/mouse_patchseq_VISp_20210403_collapsed40_cpm/"
facsDir  = "\\\\allen/programs/celltypes/workgroups/rnaseqanalysis/shiny/facs_seq/mouse_V1_ALM_20180520/"
outputFolder = "\\\\allen/programs/celltypes/workgroups/hct/HCT_RNAseq/Jeremy/patchseq_analysis/mouse_IVSCC_methods_paper/"
setwd(outputFolder)

Load and subsample the FACs data to 100 cells per cluster, excluding all ALM and outlier clusters.

Note that the anno.feather and data.feather files for the FACS data below can be generated by following the Code 1. Download and prepare the data script in the VENcelltypes R library available here: https://github.com/AllenInstitute/L5_VEN. Just follow the instructions for mouse VISp (the other parts can be ignored).

# Read in the data --> NOTE, this can be generated from the website
Samp.dat <- read_feather(paste(facsDir,"anno.feather",sep="")) 
Expr.dat <- feather(paste(facsDir,"data.feather",sep=""))   # FPKM
Samp.dat <- Samp.dat[match(Expr.dat$sample_id,Samp.dat$sample_id),] 

# Define a second annotation and data file with all VISp clusters
ld        <- sort(unique(Samp.dat$cluster_label))
useClust2 <- ld
for (val in c("ALM","Batch Grouping","Doublet","High Intron","Low Quality"))
  useClust2 <- as.character(useClust2[!grepl(val,useClust2)])
kpSamp2   <- subsampleCells(Samp.dat$subclass_label,100)
kpSamp2   <- kpSamp2&is.element(Samp.dat$cluster_label,useClust2)

gc()

          used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells 2153628 115.1    3609878 192.8  3609878 192.8
Vcells 7278972  55.6   13134686 100.3 10875259  83.0

annoFACs2 <- Samp.dat[kpSamp2,]
datFACs2  <- as.matrix(Expr.dat[kpSamp2,names(Expr.dat)!="sample_id"])
rownames(datFACs2) <- annoFACs2$sample_id
datFACs2  <- t(datFACs2)
annoFACs2$subclass_label = make.names(annoFACs2$subclass_label)

# Define class labels
classBr   <- annoFACs2$subclass_label
classBr[annoFACs2$class_label!="Non-Neuronal"] = annoFACs2$class_label[annoFACs2$class_label!="Non-Neuronal"]
classBr   <- factor(classBr)
clustersF <- factor(annoFACs2$subclass_label)

gc()

           used  (Mb) gc trigger   (Mb)   max used   (Mb)
Ncells  2156780 115.2    3609878  192.8    3609878  192.8
Vcells 92272499 704.0  829234220 6326.6 1294913086 9879.5

Load the Patch-seq data and filter to only include those cells that will be used in the IVSCC methods manuscript. This includes cells from Vip, Pvalb, Sst, and Rbp4 Cre lines.

# Read in the data
Samp.datp <- read_feather(paste(patchDir,"anno.feather",sep="")) 
Expr.datp <- feather(paste(patchDir,"data.feather",sep=""))   # FPKM
Samp.datp <- Samp.datp[match(Expr.datp$sample_id,Samp.datp$sample_id),] 

# Read in sample subset and check to make sure they are all there
cells_ivscc <- read.csv("cells_for_IVSCC_methods_all.csv")  # "cells_for_IVSCC_methods.csv"
print(c(length(as.character(cells_ivscc$specimen_id)),length(intersect(as.character(cells_ivscc$specimen_id),Samp.datp$spec_id_label))))

[1] 2717 2717

# Filter data
kpSampP    <- match(as.character(cells_ivscc$specimen_id),Samp.datp$spec_id_label)
#1:dim(Samp.datp)[1] # is.element(Samp.datp$collection_label,c("Patch-seq Production","Patch-seq Pre-production")) # Use all cells

annoPat_all<- Samp.datp[kpSampP,]
annoPat_all$dendcluster_color = annoPat_all$cluster_color  # I don't know why/if we need this line of code...
datPat_all <- as.matrix(Expr.datp[kpSampP,names(Expr.datp)!="sample_id"])
rownames(datPat_all) <- annoPat_all$sample_id
datPat_all <- t(datPat_all)

gc()

            used   (Mb) gc trigger   (Mb)   max used   (Mb)
Ncells   2226889  119.0    3609878  192.8    3609878  192.8
Vcells 221936526 1693.3  797735650 6086.3 1294913086 9879.5

Define marker genes for each broad class and contamination class (use 50 total for now). These are selected using some reasonable approach which involves a combination of median expression per class and proportion of cells per class expressing a given gene. The approach is specifically what was published in the patchSeqQC library, with the only difference being that the specific reference data set was changed to our mouse VISp/ALM data set. Here is the description of on and off markers: ‘On markers, are genes that are highly expressed in the cell type of interest with enriched expression relative to other cell types. The second class, Off markers, are expected to be expressed at low levels in a given patch-seq cell type.’

For this analysis, we build “on” markers by subclass, and “off” markers by subclass for Non-neuronal cells and by class for neuronal cells. This approximates what was done in Shreejoy’s paper.

markers     <- defineClassMarkers(datFACs2,clustersF,classBr,numMarkers = 50)
allMarkers  <- unique(unlist(markers))
markerTable <- NULL
for (i in 1:length(markers)) markerTable <- cbind(markerTable,markers[[i]])
colnames(markerTable) <- names(markers)
write.csv(markers, paste0(outputFolder,"class_markers.csv"),row.names=FALSE)

Format the data for FACS and patch-seq.

# Format the data for FACS and patch-seq
tmp                 <- datFACs2
rownames(tmp)       <- make.names(rownames(tmp))
facs_df             <- as.data.frame(t(tmp[allMarkers,])+1)
facs_df$sample_id   <- rownames(facs_df)
facs_df$major_type  <- as.character(classBr)
facs_df$contam_type <- as.character(clustersF)

tmp                 <- datPat_all
rownames(tmp)       <- make.names(rownames(tmp))
pat_df              <- as.data.frame(t(tmp[allMarkers,annoPat_all$sample_id])+1)
pat_df$sample_id    <- rownames(pat_df)

Define which subclass each patch-seq cell is assigned to, based on maximal marker expression.

nm          <- names(markers)
isOn        <- substr(nm,nchar(nm)-2,nchar(nm))=="_on"
useThese    <- nm[isOn&(!is.element(nm,paste0(nm,"_on")))]
useThese    <- setdiff(useThese,c("CR_on","Meis2_on"))        # These types are rare and unlikely to be actually patched.
subclassDat <- calcContamAllTypes(pat_df, markers[useThese])  # Identify subclass based on marker gene expression
subclass    <- colnames(subclassDat)[subclassDat %>% apply(1,which.max)]
subclass    <- gsub("_on","",subclass)

pat_df$contam_type <- subclass
tmp2 = match(pat_df$contam_type,annoFACs2$subclass_label)
pat_df$major_type  <- as.character(classBr)[tmp2]
pat_df$contam_type <- paste0(pat_df$contam_type,"_on")

Calculate the QC metrics

This is also a wrapper function for patchSeqQC which includes quality score, contamination score, and marker gene sum (both normalized and un-normalized). Here are some approximate definitions for these things:
* norm_marker_sum - This is a measure of how much expression of expected marker genes are found in a patch-seq cell relative to what is seen in the FACs data from which the markers were derived (more detail below.) * contam_sum - the contamination index for cell c (of type A), reflects off-target contamination across multiple broad cell types * quality_score - “we correlated each patch-seq sample’s expression of “on” and “off” marker genes with the average expression profile of dissociated cells of the same type’

(Details on norm marker sum: “We summarized the expression of multiple cell type specific markers specific to cell type B (MarkersB), in a cell c of type A as: Mc_A, B=∑m∈MarkersBlog2(cm). Where cm denotes the normalized expression of marker gene m in cell c. We further used the dissociated-cell reference data to quantify how much marker expression of cell type B’s markers one would typically expect in cells of type A as: dA_B=mediantypeA(Mc_A, B). Reflecting the median marker expression of cell type B’s markers in dissociated cells of type A.”)

qcMetrics <- calculatePatchSeqQCMetrics2(pat_df,facs_df,markers)

the standard deviation is zero

# We are using NMS score of 0.4 as a pass/fail call
qcMetrics$QC_pass <- c(TRUE,FALSE)[(qcMetrics$marker_sum_norm<0.40)+1]

Plot subclass vs. NMS pass

First group by NMS pass/fail.

# Merge sparse subclasses
class2 <- subclass
class2[is.element(class2,c("Serpinf1","Sncg","Lamp5"))] = "Inh_other"
class2[is.element(class2,c("L2.3.IT","L6.CT","L6.IT","NP"))] = "Exc_other"
class2 <- factor(class2,levels=c("Pvalb","Sst","Vip","Inh_other","L4","L5.IT","L5.PT","Exc_other"))
geno   <- factor(cells_ivscc$Cre, levels=c("Pvalb","Sst","Vip","Rbp4"))

datPlot <- data.frame(cell=annoPat_all$spec_id_label,
                      NMS_pass=factor(qcMetrics$QC_pass,levels=c(TRUE,FALSE)), 
                      genotype=geno, 
                      subclass=class2)
head(datPlot)

g1 <- ggplot(datPlot, aes(x= subclass,  group=NMS_pass)) + 
    geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") +
    geom_text(aes( label = scales::percent(..prop..,accuracy = 0.1), 
                   y= ..prop.. ), stat= "count", size=3, vjust = -.5) +
    labs(y = "Percent", fill="subclass") +
    facet_grid(genotype ~ NMS_pass) +
    scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0,0.1), add = 0)) +
    scale_x_discrete(expand = expansion(add = 1)) + 
    theme(axis.text.x = element_text(angle = 90))
g1

ggsave("subclass_plot.pdf",height=9,width=8)

Now subset for GABA cells only and bin by 0.1 bins.

kp    <- (geno!="Rbp4")
match <- as.character(geno)==as.character(class2)

datPlot <- data.frame(cell=annoPat_all$spec_id_label,
                      NMS_score=floor(10*qcMetrics$marker_sum_norm)/10, 
                      genotype=geno, 
                      match=match)[kp,]
head(datPlot)

tabPlot <- table(datPlot$NMS_score,datPlot$match,droplevels(datPlot$genotype))
for (i in 1:3) tabPlot[,,i] <- tabPlot[,,i]/rowSums(tabPlot[,,i])
write.csv(tabPlot[,2,],"heatmap_values.csv")

out <- tabPlot[,2,]
out <- out[dim(out)[1]:1,]
col <- WGCNA::blueWhiteRed(50)[30:50]

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

gplots::heatmap.2(out,scale="none", Rowv = NA, Colv=NA, col=col, cellnote=round(100*out),
                  notecol="black",trace="none",notecex=0.7,margins=c(5,15))

Discrepancy: Rowv is FALSE, while dendrogram is `both'. Omitting row dendogram.Discrepancy: Colv is FALSE, while dendrogram is `column'. Omitting column dendogram.

pdf("subclass_heatmap.pdf")

gplots::heatmap.2(out,scale="none", Rowv = NA, Colv=NA, col=col, cellnote=round(100*out),
                  notecol="black",trace="none",notecex=0.7,margins=c(8,20))

Discrepancy: Rowv is FALSE, while dendrogram is `both'. Omitting row dendogram.Discrepancy: Colv is FALSE, while dendrogram is `column'. Omitting column dendogram.

dev.off()

png 
  2 

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scrattch.io_0.1.0   rhdf5_2.28.0        patchSeqQC_0.1.0    patchseqtools_0.0.1 mfishtools_0.0.2    ggplot2_3.3.3       viridisLite_0.3.0  
 [8] bindrcpp_0.2.2      dplyr_1.0.5         feather_0.3.5      

loaded via a namespace (and not attached):
 [1] bitops_1.0-6          matrixStats_0.58.0    bit64_4.0.5           doParallel_1.0.16     RColorBrewer_1.1-2    dynamicTreeCut_1.63-1
 [7] backports_1.2.1       tools_3.6.1           bslib_0.2.4           utf8_1.1.4            R6_2.5.0              KernSmooth_2.23-18   
[13] rpart_4.1-15          Hmisc_4.5-0           DBI_1.1.1             lazyeval_0.2.2        BiocGenerics_0.30.0   colorspace_2.0-0     
[19] nnet_7.3-15           withr_2.4.1           tidyselect_1.1.0      gridExtra_2.3         preprocessCore_1.46.0 bit_4.0.4            
[25] compiler_3.6.1        WGCNA_1.70-3          Biobase_2.44.0        htmlTable_2.1.0       labeling_0.4.2        sass_0.3.1           
[31] caTools_1.18.1        scales_1.1.1          checkmate_2.0.0       stringr_1.4.0         digest_0.6.27         foreign_0.8-71       
[37] rmarkdown_2.7         base64enc_0.1-3       jpeg_0.1-8.1          pkgconfig_2.0.3       htmltools_0.5.1.1     fastmap_1.1.0        
[43] htmlwidgets_1.5.3     rlang_0.4.10          impute_1.58.0         rstudioapi_0.13       RSQLite_2.2.3         bindr_0.1.1          
[49] jquerylib_0.1.3       farver_2.1.0          generics_0.1.0        jsonlite_1.7.2        gtools_3.8.2          magrittr_2.0.1       
[55] GO.db_3.8.2           Formula_1.2-4         Matrix_1.3-2          Rcpp_1.0.6            munsell_0.5.0         S4Vectors_0.22.1     
[61] Rhdf5lib_1.6.1        fansi_0.4.2           lifecycle_1.0.0       stringi_1.5.3         yaml_2.2.1            gplots_3.1.1         
[67] grid_3.6.1            blob_1.2.1            parallel_3.6.1        crayon_1.4.1          lattice_0.20-41       splines_3.6.1        
[73] hms_1.0.0             knitr_1.31            pillar_1.5.1          fastcluster_1.1.25    codetools_0.2-18      stats4_3.6.1         
[79] glue_1.4.2            evaluate_0.14         latticeExtra_0.6-29   data.table_1.14.0     vctrs_0.3.6           png_0.1-7            
[85] foreach_1.5.1         gtable_0.3.0          purrr_0.3.4           assertthat_0.2.1      cachem_1.0.4          xfun_0.23            
[91] survival_3.2-7        tibble_3.1.0          iterators_1.0.13      AnnotationDbi_1.46.1  memoise_2.0.0         IRanges_2.18.3       
[97] cluster_2.1.1         ellipsis_0.3.1