Last updated: 2019-08-08

Checks: 7 0

Knit directory: scRNA-seq-workshop-Fall-2019/

This reproducible R Markdown analysis was created with workflowr (version 1.4.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20190718) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  data/pbmc10k/
    Untracked:  data/pbmc5k/

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File Version Author Date Message
Rmd 0cf484c Ming Tang 2019-08-08 wflow_publish(c(“analysis/scRNAseq_workshop_0.Rmd”, “analysis/scRNAseq_workshop_1.Rmd”,
html b27f0e8 Ming Tang 2019-07-30 Build site.
html ac73193 Ming Tang 2019-07-30 Build site.
html 614d882 Ming Tang 2019-07-29 Build site.

From sequencer to cellranger

In this section, I will show you how to prepare the fastq files and count the scRNAseq matrix by cellranger. After sequencing, one usually gets a folder from the sequencing core with a folder structure like:

The bcl (Binary Base Call) files in the Data folder contains the raw data generated from the illumina sequencers. cellranger wraps the illumina bcl2fastq command into cellranger mkfastq to convert it to fastq files for single-cell RNAseq data.

cellranger mkfastq

For details, check the tutorial from 10x Genoimcs.

On Odyssey computing cluster:

module load bcl2fastq2
cellranger mkfastq --id=test \
                   --run=/path/to/the/run/folder \
                   --csv=test.csv \
                   --jobmode=local \
                   --localmem=40 \
                   --localcores=12

test.csv is a comma seprated file with three columns:

Lane,Sample,Index
1,test_sample,SI-GA-A3

cellranger count

After cellranger mkfastq, we are ready to align the fastqs to the reference genome and count how many reads per gene per cell. These steps are wraped in cellranger count command.

cellranger count --id=sample345 \
                   --transcriptome=/opt/refdata-cellranger-GRCh38-3.0.0 \
                   --fastqs=/home/test/outs/fastq_path/HAWT7ADXX/test_sample/ \
                   --sample=mysample \
                   --expect-cells=6000

What does the output of cellranger count look like?

In the sample345 folder there is an outs folder, and you will find the files Seurat works with in the filtered_feature_bc_matrix folder. There are 3 files in the folder:

ls -sh filtered_feature_bc_matrix/
total 90M
 60K barcodes.tsv.gz  300K features.tsv.gz   90M matrix.mtx.gz
 
# The `barcodes.tsv.gz` contains the cell barcode that passed the `cellranger` filter.
zcat barcodes.tsv.gz | head -5
AAACCCAAGCGCCCAT-1
AAACCCAAGGTTCCGC-1
AAACCCACAGAGTTGG-1
AAACCCACAGGTATGG-1
AAACCCACATAGTCAC-1

# how many cells (barcodes)?
zcat barcodes.tsv.gz | wc -l
11769

# The `features.tsv.gz` contains the ENSEMBLE id and gene symbol
zcat features.tsv.gz | head -5
ENSG00000243485 MIR1302-2HG     Gene Expression
ENSG00000237613 FAM138A Gene Expression
ENSG00000186092 OR4F5   Gene Expression
ENSG00000238009 AL627309.1      Gene Expression
ENSG00000239945 AL627309.3      Gene Expression

## how many genes?
zcat features.tsv.gz | wc -l
33538


# matrix.mtx.gz is a sparse matrix which contains the non-zero counts
zcat matrix.mtx.gz | head -10
%%MatrixMarket matrix coordinate integer general
%metadata_json: {"format_version": 2, "software_version": "3.0.0"}
33538 11769 24825783
33509 1 1
33506 1 4
33504 1 2
33503 1 10
33502 1 5
33500 1 20
33499 1 9

Most of the entries in the final gene x cell count matrix are zeros. Sparse matrix efficiently save the disk space by only recording the non-zero entries.

You see the dimension of the matrix is 33538 x 11769 and the number of non-zero entries is 24825783

e.g. for the subsequent two rows in the sparse matrix:

33509 1 is the index of the row (gene) and column(cell) of that non-zero entry in the matrix, and 1 is the count number.

33506 1 is the index of the row and column of that non-zero entry in the matrix, and 4 is the count number.

Alternatives to cellranger

cellranger is very slow. It can take several days to run a mouse single-cell RNAseq data set with even 20 CPUs. There are other tools which can process single-cell RNAseq data set much faster and accurate as well.

Alevin

Paper : Alevin efficiently estimates accurate gene abundances from dscRNA-seq data

It supports Drop-seq and 10x-Chromium v1/2/3.

Tutorial

Kallisto/bustools

Paper: Modular and efficient pre-processing of single-cell RNA-seq

the bustoools commands we implemented are generic and will work with any BUS file, generated with data from any scRNA-seq technology. Distinct technology encodes barcode and UMI information differently, but the kallisto bus command can accept custom formatting rules. While the pre-processing steps for error correction and counting may need to be optimized for the distinguishing characteristics of different technologies, the modularity of the bustools based workflow makes such customization possible and easy.

Tutorial

Scumi

scumi is a flexible Python package to process fastq files generated from different single-cell RNA-sequencing (scRNA-seq) protocols to produce a gene-cell sparse expression matrix for downstream analyses, e.g., discovering cell types and inferring cell lineages.

It supports CEL-Seq2, 10x Chromium, Drop-seq, Seq-Well, CEL-Seq2, inDrops, and SPLiT-seq.

Public data sets

If you want to analyze public single-cell data, there are databases you can go to:

  • scRNAseq bioc package Gene-level counts for a collection of public scRNA-seq datasets, provided as SingleCellExperiment objects with cell- and gene-level metadata.
  • human cell atlas database
  • SRA
  • EMBL-EBI atlas
  • PanglaoDB is a database for the scientific community interested in exploration of single cell RNA sequencing experiments from mouse and human. We collect and integrate data from multiple studies and present them through a unified framework.
  • scRNASeqDB database, which contains 36 human single cell gene expression data sets collected from Gene Expression Omnibus (GEO)
  • JingleBell A repository of standardized single cell RNA-Seq datasets for analysis and visualization at the single cell level.
  • Broad single cell portal
  • The conquer (consistent quantification of external rna-seq data) repository is developed by Charlotte Soneson and Mark D Robinson at the University of Zurich, Switzerland. It is implemented in shiny and provides access to consistently processed public single-cell RNA-seq data sets.

Steps to follow this tutorial

Click the Terminal tab in your Rstudio

# go to your home directory
cd ~

# clone the github repo
git clone https://github.com/crazyhottommy/scRNA-seq-workshop-Fall-2019

go to Files —> scRNA-seq-workshop-Fall-2019 folder and double click the scRNA-seq-workshop-Fall-2019.Rproj file. Inside the analysis folder, you can find all the Rmd files which you can execute the commands.


sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

loaded via a namespace (and not attached):
 [1] workflowr_1.4.0 Rcpp_1.0.0      digest_0.6.18   rprojroot_1.3-2
 [5] backports_1.1.3 git2r_0.23.0    magrittr_1.5    evaluate_0.12  
 [9] stringi_1.2.4   fs_1.2.6        whisker_0.3-2   rmarkdown_1.11 
[13] tools_3.5.1     stringr_1.3.1   glue_1.3.0      xfun_0.4       
[17] yaml_2.2.0      compiler_3.5.1  htmltools_0.3.6 knitr_1.21