/
scRNAseq_workshop_0.Rmd
168 lines (117 loc) · 7.12 KB
/
scRNAseq_workshop_0.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
---
title: "Part 0"
output: html_document
---
# 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`](https://support.illumina.com/sequencing/sequencing_software/bcl2fastq-conversion-software.html) command into `cellranger mkfastq` to convert it to fastq files for single-cell RNAseq data.
### cellranger mkfastq
For details, check the [tutorial](https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger) from 10x Genoimcs.
On `Odyssey` computing cluster:
```{bash eval =FALSE}
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.
```{bash eval =FALSE}
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:
```{bash eval =FALSE}
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](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1670-y)
It supports Drop-seq and 10x-Chromium v1/2/3.
[Tutorial](https://salmon.readthedocs.io/en/latest/alevin.html)
### Kallisto/bustools
Paper: [Modular and efficient pre-processing of single-cell RNA-seq](https://www.biorxiv.org/content/10.1101/673285v2)
>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](https://www.kallistobus.tools/)
### Scumi
[scumi](https://bitbucket.org/jerry00/scumi-dev/src/master/) 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](https://bioconductor.org/packages/devel/data/experiment/html/scRNAseq.html) 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](https://staging.data.humancellatlas.org/)
* [SRA](https://www.ncbi.nlm.nih.gov/sra)
* [EMBL-EBI atlas](https://www.ebi.ac.uk/gxa/sc/home)
* [PanglaoDB](https://panglaodb.se/) 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](https://bioinfo.uth.edu/scrnaseqdb/) database, which contains 36 human single cell gene expression data sets collected from Gene Expression Omnibus (GEO)
* [JingleBell](http://jinglebells.bgu.ac.il/) A repository of standardized single cell RNA-Seq datasets for analysis and visualization at the single cell level.
* [Broad single cell portal](https://portals.broadinstitute.org/single_cell)
* The [conquer](http://imlspenticton.uzh.ch:3838/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
```bash
# 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.