Introduction

This page will walk through the steps to go from the raw 10x sequencing fastq files to a count matrix and demuxlet assignment of droplets to individuals. This involves running a Snakemake pipeline located in the code directory and some code in R. The end product are seurat objects containing the raw count matrices and assignment to one of the three individuals and two treatment groups (strain or unstrain) for each of the two 10x lanes used in the sequencing experiment.

Run the Snakemake Pipeline after downloading necessary files

1. Move the fastq files from all samples (these are output files from 10x runs, containing both forward and reverse sequences) into the folder code/single_cell_preprocessing/fastq/. Undetermined data files are not required.

2. Modify the code/single_cell_preprocessing/Pipeline/cluster_solo.json file to correspond to the computing cluster you are working with.

3. Unzip the whitelist file in the code/single_cell_preprocessing/ directory.

4. Place the human.YRI.cellranger.exons.vcf file into the code/single_cell_preprocessing/ directory.

5. Install the conda working by running “conda env create –file environment.yaml”

6. Run the file “submit.sh”.

Load files into R

After running the pipeline, two directories will be created corresponding to the two 10x lanes involved in the sequencing experiment (“YG-AH-2S-ANT-1_S1_L008” and “YG-AH-2S-ANT-2_S2_L008”). These directories contain the demuxlet and STAR SOLO outputs for each 10x lane.

Details about the individualXtreatment status of the samples that were pooled for each of the two lanes: ANT1: (NA19160 unstrained; NA18856 unstrained; NA18855 strained) ANT2: (NA19160 strained; NA18855 unstrained)

library(data.table)
library(Matrix)
library(Seurat)
library(readr)
library(stringr)
library(plyr)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:plyr':

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize

The following objects are masked from 'package:data.table':

    between, first, last

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

## link to directories containing data files (count matrices)
proj_dir <- "code/single_cell_preprocessing/"
ANT1_dir <- paste0(proj_dir, "YG-AH-2S-ANT-1_S1_L008/")
ANT2_dir <- paste0(proj_dir, "YG-AH-2S-ANT-2_S2_L008/")

#read in data

##Gene Output from STARSOLO
#ANT1
demuxlet1 <- fread(paste0(ANT1_dir, "demuxlet.best", sep = ""))
count_data1 <- readMM(paste0(ANT1_dir, "Gene/filtered/matrix.mtx"))
genes1 <- read_tsv(paste0(ANT1_dir, "Gene/filtered/genes.tsv"), col_names = F)

Parsed with column specification:
cols(
  X1 = col_character(),
  X2 = col_character()
)

barcodes1 <- as.data.frame(read_tsv(paste0(ANT1_dir, "Gene/filtered/barcodes.tsv"), col_names = F))

Parsed with column specification:
cols(
  X1 = col_character()
)

#ANT2
demuxlet2 <- fread(paste0(ANT2_dir, "demuxlet.best", sep = ""))
count_data2 <- readMM(paste0(ANT2_dir, "Gene/filtered/matrix.mtx"))
genes2 <- read_tsv(paste0(ANT2_dir, "Gene/filtered/genes.tsv"), col_names = F)

Parsed with column specification:
cols(
  X1 = col_character(),
  X2 = col_character()
)

barcodes2 <- as.data.frame(read_tsv(paste0(ANT2_dir, "Gene/filtered/barcodes.tsv"), col_names = F))

Parsed with column specification:
cols(
  X1 = col_character()
)

Based on the demuxlet output, assign label for barcodes based on “BEST” output and filter for “SNG-” barcodes

#this function returns a dataframe with two columns, one corresponding to the barcodes and one corresponding to the label given by demuxlet
return_singlet_label <- function(barcodes, demuxlet.out){
  labels <- demuxlet.out$BEST[match(unlist(barcodes), demuxlet.out$BARCODE)]
  return(cbind(barcodes, labels))
}

barcodes1_labeled <- return_singlet_label(barcodes1, demuxlet1)
barcodes2_labeled <- return_singlet_label(barcodes2, demuxlet2)

#table of singlets/multiplets in the filtered data based on demuxlet
table(barcodes1_labeled$labels)

DBL-NA18855-NA18856-0.500 DBL-NA19160-NA18855-0.500 
                       13                         6 
DBL-NA19160-NA18856-0.500               SNG-NA18855 
                        7                       411 
              SNG-NA18856               SNG-NA19160 
                      452                       260 

table(barcodes2_labeled$labels)

DBL-NA18855-NA18856-0.500 DBL-NA18855-NA19160-0.500 
                       20                      2719 
DBL-NA18856-NA18855-0.500 DBL-NA19160-NA18855-0.500 
                        4                      5650 
DBL-NA19160-NA18856-0.500               SNG-NA18855 
                       22                      4990 
              SNG-NA19160 
                     1506 

## filter for droplets in the count data that are singlets (remove multiplets)
#ANT1
demuxlet_single1 <- demuxlet1 %>%
     dplyr::filter(grepl("SNG-", BEST))
singlets_index1 <- unlist(lapply(barcodes1_labeled$X1,"%in%", table = demuxlet_single1$BARCODE), use.names = F) #get index of barcodes that are singlets
barcodes_singlets1 <- barcodes1_labeled[singlets_index1,] #use index to subset matrix + barcode names
count_data_singlets1 <- count_data1[,singlets_index1]

#ANT2
demuxlet_single2 <- demuxlet2 %>%
     dplyr::filter(grepl("SNG-", BEST))
singlets_index2 <- unlist(lapply(barcodes2_labeled$X1,"%in%", table = demuxlet_single2$BARCODE), use.names = F) #get index of barcodes that are singlets
barcodes_singlets2 <- barcodes2_labeled[singlets_index2,] #use index to subset matrix + barcode names
count_data_singlets2 <- count_data2[,singlets_index2]

Create Seurat object for each dataset (for singlet barcodes) and add metadata in the form of singlet identity for each barcode.

#Change labels to reflect strain/unstrain based on prior knowledge of which strainXindividual combinations went into each pool

strainIndlabels1 <- revalue(barcodes_singlets1$labels,
                            c("SNG-NA18856"= "NA18856_Unstrain",
                              "SNG-NA18855" = "NA18855_Strain",
                              "SNG-NA19160" = "NA19160_Unstrain"))

strainIndlabels2 <- revalue(barcodes_singlets2$labels,
                            c("SNG-NA18855" = "NA18855_Unstrain",
                              "SNG-NA19160" = "NA19160_Strain"))


rownames(count_data_singlets1) <- genes1$X2
colnames(count_data_singlets1) <- barcodes_singlets1$X1

ANT1_seurat <- CreateSeuratObject(counts = count_data_singlets1, project = "ANT1") %>%
  AddMetaData(strainIndlabels1, col.name = "labels")

Warning: Non-unique features (rownames) present in the input matrix, making
unique

rownames(count_data_singlets2) <- genes2$X2
colnames(count_data_singlets2) <- barcodes_singlets2$X1

ANT2_seurat <- CreateSeuratObject(counts = count_data_singlets2, project = "ANT2") %>%
  AddMetaData(strainIndlabels2, col.name = "labels")

Warning: Non-unique features (rownames) present in the input matrix, making
unique

Save Seurat objects

saveRDS(ANT1_seurat, "data/ANT1.rds")
saveRDS(ANT2_seurat, "data/ANT2.rds")

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] dplyr_1.0.2       plyr_1.8.6        stringr_1.4.0     readr_1.3.1      
[5] Seurat_3.2.3      Matrix_1.2-18     data.table_1.13.0

loaded via a namespace (and not attached):
  [1] Rtsne_0.15            colorspace_2.0-0      deldir_0.1-23        
  [4] ellipsis_0.3.1        ggridges_0.5.1        rprojroot_2.0.2      
  [7] fs_1.3.1              spatstat.data_1.7-0   leiden_0.3.1         
 [10] listenv_0.8.0         npsurv_0.4-0          ggrepel_0.9.0        
 [13] codetools_0.2-16      splines_3.6.1         lsei_1.2-0           
 [16] knitr_1.23            polyclip_1.10-0       jsonlite_1.7.2       
 [19] workflowr_1.6.2       ica_1.0-2             cluster_2.1.0        
 [22] png_0.1-7             uwot_0.1.10           shiny_1.3.2          
 [25] sctransform_0.3.2     compiler_3.6.1        httr_1.4.2           
 [28] lazyeval_0.2.2        later_1.1.0.1         htmltools_0.5.0      
 [31] tools_3.6.1           rsvd_1.0.1            igraph_1.2.4.1       
 [34] gtable_0.3.0          glue_1.4.2            RANN_2.6.1           
 [37] reshape2_1.4.3        rappdirs_0.3.1        Rcpp_1.0.5           
 [40] spatstat_1.64-1       scattermore_0.7       vctrs_0.3.6          
 [43] gdata_2.18.0          nlme_3.1-140          lmtest_0.9-37        
 [46] xfun_0.8              globals_0.12.5        mime_0.9             
 [49] miniUI_0.1.1.1        lifecycle_0.2.0       irlba_2.3.3          
 [52] gtools_3.8.1          goftest_1.2-2         future_1.18.0        
 [55] MASS_7.3-52           zoo_1.8-8             scales_1.1.1         
 [58] hms_0.5.3             promises_1.1.1        spatstat.utils_1.17-0
 [61] parallel_3.6.1        RColorBrewer_1.1-2    yaml_2.2.1           
 [64] reticulate_1.16       pbapply_1.4-0         gridExtra_2.3        
 [67] ggplot2_3.3.3         rpart_4.1-15          stringi_1.4.6        
 [70] caTools_1.17.1.2      rlang_0.4.10          pkgconfig_2.0.3      
 [73] matrixStats_0.57.0    bitops_1.0-6          evaluate_0.14        
 [76] lattice_0.20-41       ROCR_1.0-7            purrr_0.3.4          
 [79] tensor_1.5            patchwork_1.1.0       htmlwidgets_1.5.2    
 [82] cowplot_1.1.0         tidyselect_1.1.0      RcppAnnoy_0.0.18     
 [85] magrittr_2.0.1        R6_2.5.0              gplots_3.0.1.1       
 [88] generics_0.0.2        pillar_1.4.7          whisker_0.3-2        
 [91] mgcv_1.8-28           fitdistrplus_1.0-14   survival_2.44-1.1    
 [94] abind_1.4-5           tibble_3.0.4          future.apply_1.3.0   
 [97] crayon_1.3.4          KernSmooth_2.23-15    plotly_4.9.2.1       
[100] rmarkdown_1.13        grid_3.6.1            git2r_0.26.1         
[103] digest_0.6.27         xtable_1.8-4          tidyr_1.1.2          
[106] httpuv_1.5.1          munsell_0.5.0         viridisLite_0.3.0