TREX is an experimental workflow that enables simultaneous lineage TRacking and EXpression profiling of single cells using RNA-sequencing. The method is described in the paper Clonal relations in the mouse brain revealed by single-cell and spatial transcriptomics.
An essential part of this workflow is presented here: the extraction of genetic barcodes or "cloneIDs" from single-cell or spatial transcriptomes and the reconstruction of related cells/spots.
The tool uses BAM files of one or multiple sequencing libraries as an input for the generation of cloneID count matrices and identifies clonally related cells based on Jaccard similarity between each pair of cloneID+cells.
Currently, TREX is compatible with common RNA-sequencing library preparation methods and data formats provided by 10X Chromium, 10X Visium, Smart-seq2 and 3.
TREX requires Python 3.7 or newer. We recommend that you install TREX into a "virtual environment", which can be done by running the following commands:
python3 -m venv trex-venv
trex-venv/bin/pip install git+https://github.com/frisen-lab/TREX.git
trex-venv
is the name of the directory that will be created and which will contain the virtual environment.
You can choose a different name.
Activate the virtual environment by running
source trex-venv/bin/activate
.
You need to repeat this step in every new shell in order to be able to run TREX.
Finally, test the installation by running trex --version
.
See Changelog.
Clone the Git repository or download it as a ZIP file and unpack it.
The directory tests/data/
contains a test dataset in the proper format.
Run TREX on it:
trex run10x -s 695 -e 724 tests/data/
This will create a folder trex_run/
in the current directory
(use the -o
option to choose a different folder) with the results.
See the "Runnning TREX" section below for further details.
Here we show how to run TREX on the data we analyzed in our Nature Neuroscience paper (https://doi.org/10.1038/s41593-022-01011-x).
These instructions have been tested with Cell Ranger 6.1.2, but the original analysis was done with Cell Ranger 2.2.0.
The data is available under GEO accession GSE153424. The instructions below analyze sample GSM4644060 because it is relatively small.
GSM4644060 is also called "brain1_str" (str stands for striatum) in the GEO description, and it has the ID "10X_41" in Supplementary Table 4 in the paper.
As can be seen on the overview page for GSM4644060, reads for the dataset are available from SRA experiment accession SRX8627776, which in turn links to two run accessions:
We use the run accessions to retrieve the data.
Go to the Cell Ranger download page. Download and install Cell Ranger.
Download and extract the mouse reference dataset refdata-gex-mm10-2020-A.tar.gz
:
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz tar xvf refdata-gex-mm10-2020-A.tar.gz
Create a custom reference by adding an extra chrH2B-EGFP-N contig (EGFP-cloneID virus) to the above mouse reference dataset.
references/
refers to the directory in this repository. cellranger mkref
takes about one hour.
You may choose to continue to the next step (downloading reads) while this is running.
cat refdata-gex-mm10-2020-A/fasta/genome.fa references/chrH2B-EGFP-N.fa > mm10_H2B-EGFP-30N_genome.fa cat refdata-gex-mm10-2020-A/genes/genes.gtf references/chrH2B-EGFP-N.gtf > mm10_H2B-EGFP-30N_genes.gtf cellranger mkref --genome=mm10_H2B-EGFP-30N --fasta=mm10_H2B-EGFP-30N_genome.fa --genes=mm10_H2B-EGFP-30N_genes.gtf > mkref_mm10_H2B-EGFP-30N.out
Create a fastq
directory and change into it:
mkdir fastq
cd fastq
If you do not have it, install fastq-dump
. For example, if you have Conda (with the Bioconda channels activated), run
conda create -n trex sra-tools
conda activate trex
Download the reads:
fastq-dump --gzip --split-3 --defline-qual '+' --defline-seq '@$ac.$sn' SRR12103475 SRR12103476
Rename the files so that Cell Ranger can find them:
mv SRR12103475_1.fastq.gz SRR12103475_S1_L001_R1_001.fastq.gz
mv SRR12103475_2.fastq.gz SRR12103475_S1_L001_R2_001.fastq.gz
mv SRR12103476_1.fastq.gz SRR12103476_S1_L001_R1_001.fastq.gz
mv SRR12103476_2.fastq.gz SRR12103476_S1_L001_R2_001.fastq.gz
cd ..
Continue with this step only when the above steps have finished.
Run cellranger count
:
cellranger count --transcriptome=mm10_H2B-EGFP-30N --id=brain1_str --fastqs=fastq/ --sample=SRR12103475 --sample=SRR12103476 --expect-cells=2299
The --expect-cells
parameter is set to the number of cells loaded per well in the 10X chip used for droplet generation (4000 in this case, see the "cells recovered for 10X" column in Suppl. Table 4) divided by 1.74.
With trex
installed (as described above), run
trex run10x -o trex_brain1_str brain1_str
Results will be written to a new directory named trex_brain1_str
.
The input directory for TREX must be a Cell Ranger output directory. In case of Smart-Seq2 / 3 data, one BAM file with all cells or a folder with one BAM file per cell is expected (see zUMIs output)
See the contents of the tests/data/
directory to learn which are the minimum files necessary.
Cell Ranger/zUMIs must have been configured to map against a reference augmented by an extra chromosome that contains the cloneID. By default, that extra chromosome is assumed to be the last in the BAM file (use --chromosome
to choose a different one).
The options -s
and -e
set where on the extra chromosome the cloneID is located (-s
gives start and -e
gives end in 1-based coordinates).
Please also run trex run10x --help
(or trex smartseq2 --help
and trex smartseq3 --help
respectively ) to see the other available command-line options.
This is an overview of the steps that the trex run10x
/ trex smartseq3
command performs.
--chr
, -s
and
-e
flags or, if the flags are not given, to the region that
has been automatically identified to be the region containing
the variable cloneID sequence,CB
and UB
in case of 10x data , SAM tags BC
and UB
in case of Smart-seq3 data),--filter-cellid
or -f
).trex smartseq2
follows a similar pipeline with the following differences:
The Jaccard index measures the similarity of two sample sets, in this case
the similarity of two sets of cloneIDs. It is calculated by dividing
the number of overlapping, unique cloneIDs between cell A and B by the total
number of unique cloneIDs in cell A and B. An index of 0.0 indicates no
overlapping cloneIDs and an index of 1.0 a perfect match. The Jaccard
threshold is the Jaccard index above which two cells are merged into one
clone. It can be set with the --jaccard-threshold
flag and is 0.7
by default, meaning cell A and B are merged into one clone if they have more
than 70% of cloneIDs in common.
Tab-separated file of cell IDs to keep in the TREX run. Adding this file via
the --filter-cellid
or -f
option allows to focus the analysis on specific cells
and to filter out low quality cells or doublets.
Example:
0 CACTCGTGGTACACACTCCG
1 CACTCGTGGTACCACAAGCA
Text file with cloneIDs to ignore. The format is one cloneID per line.
Adding this file via the --filter-cloneids
option allows to ignore cloneIDs that have been identified as overrepresented or artefactual during library characterization.
Example file:
GGTCTCCCTATACCAACAGTATCGTCTCAA
GGGTTCTGGGATATTACGTTGACTTGAGAG
TREX by default writes its results to a newly created directory
named trex_run
.
The name of the output directory can be changed with --output
(or -o
).
The files created in the output directory are described below.
They are listed in the order in which the program creates them,
see the "Pipeline steps overview" section for more details about what is
done in each step.
Many of the files are tables in either tab-separated value (TSV) format or in comma-separated values (CSV) format.
log.txt
This file contains a copy of the output that a TREX run prints to the terminal.
entries.bam
All usable reads from the input BAM file. (See the pipeline steps overview for a description of "usable reads".)
reads.txt
A table listing cell_id, umi and clone_id for each usable input read. Example:
#cell_id umi clone_id
TGACGGCGTTACCAGT AAAAAACTGT TGTCAATCGTTCGGTTGAGCAAGATCTTAG
0
in a cloneIDs signals a deleted base (CIGAR operation
D
in the input BAM file).-
for each missing base.molecules.txt
A table listing cell_id, umi and clone_id for each detected molecule. Example:
#cell_id umi clone_id
AAACCTGAGAGGTACC AGTTAAAGTA TGTCAATCGTTCGGTTGAGCAAGATCTTAG
molecules_corrected.txt
Same as molecules.txt
, but after filtering and error correction:
For the run10x
subcommand, this table contains an additional original_clone_id
column that shows how the original cloneID before correction.
cells.txt
A table listing the detected cells.
The fields are cell_id, a colon (:
) and then
pairs of columns clone_id1 and count1 for each cloneID found in that cell.
Example:
#cell_id : clone_id1 count1 clone_id2 count2 ...
AAACCTGAGAGGTACC : TGTCAATCGTTCGGTTGAGCAAGATCTTAG 1
AAACCTGAGTAGCGGT : TGTCAATCGTTCGGTTGAGCAAGATCTTAG 3
cells_filtered.txt
The same as cells.txt
, but with error-corrected cloneID counts.
Cells that end up without any cloneIDs after error correction are removed.
umi_count_matrix.csv
A matrix of UMI counts with cells as columns and cloneIDs as rows.
This is a different representation of the data written to cells_filtered.txt
.
Only created if option --umi-matrix
is used.
read_count_matrix.csv
Instead of UMI count matrix, running trex smartseq2
produces matrix of read counts with cells as columns and cloneIDs as rows.
This is a different representation of the data written to cells_filtered.txt
.
Only created if option --read-matrix
is used.
components.txt
Connected components of the clone graph.
components_corrected.txt
Connected components of the clone graph after error correction.
graph.pdf
, graph_corrected.pdf
Plots of the clone graph and the error-corrected clone graph.
These are only created if option --plot
is used.
graph_corrected.gv
and graph.gv
are textual descriptions of the graphs in GraphViz format.
doublets.txt
A list of the cell IDs of the cells that were detected to be doublets.
clones.txt
A table listing the cell IDs belonging to each clone. The columns are clone_id and cell_id where clone_id is a number that identifies the clone.
clone# cell_id
1 TGGCGCAAGAATAGGG
clone_details.txt
A table listing all clones (a clone is a set of cells), one line per clone. The columns are:
The columns are clone# and clone_seq where clone_id is a number that identifies the clone and clone_seq its nucleotide sequence.
clone_nr clone_id n_cells clone_ids_per_cell
1 ACTAGGAGATTGACGGATCACCTTTGGTCG 3 1.00
data.loom
A loom file.
This file is created only if option --loom
(or -l
) is used.
trex qc --plot-jaccard-matrix --plot-hamming-distance DIRECTORY
qc takes as an input the directory (or directories) of trex output.
Plotting the jaccard similarity matrix between cells requires some time as jaccard similarity is calculated pairwise amongst all cells.
This can be activated adding the optional flag --plot-jaccard-matrix
.
Hamming distance between all viral cloneIDs found in the dataset after each step can be plotted by means of the optional flag --plot-hamming-distance
.
This will add a PDF file named quality_report.pdf describing the quality of the TREX run inside the same folder with the TREX output.
This report contains:
Each of these plots has four subplots corresponding to different steps of the TREX pipeline.
It is highly recommended that you develop TREX within a separate virtual environment:
python3 -m venv --prompt trex .venv
source .venv/bin/activate
Install TREX in "editable" mode:
pip install -e .
Install pre-commit
and install the pre-commit hooks
(these run at git commit
time and do some checks on
the to-be-committed files):
pip install pre-commit
pre-commit install