# Details of the pipeline steps

Note that all rules use uncompressed intermediate files for efficiency. The size of these files is less of a concern as these are temporary files which are automatically deleted when no longer needed. The input FASTQ files may either be uncompressed or gzipped.

![Workflow rulegraph](../images/rulegraph.svg)

## `determine_bc_cutoff`
When running `determine_bc_cutoff`, three rules are run, `quality_trim`, `demultiplex`, and `fastqc`.

## `fastqc`
Completly separately from the other jobs `fastqc` is run on all input files, the results of which are found in `results/<sample>/qc`.

## `quality_trim`
`quality_trim` is a `cutadapt` command filtering out low-quality reads pairs and reads being too short. The rule is invoked once per input file pair and hence gives a pair of output files per lane. The output filenames are `results/<sample>/<sample>_QF_<lane>_R1.fastq` and `results/<sample>/<sample>_QF_<lane>_R2.fastq`. The `R1` files are fed into the demultiplexer and are typically never seen by the user, whereas the `R2` files are to be alignment and are hence kept until the `all` rule completes the alignment.

## `demultiplex`
In this step, the `posDemux` package demultiplexes the forward reads for each sample. If the sample has multiple lanes, the demultiplexer sequentially streams all the lane files in the same rule invokation.

## `all`

The rule `all` is only supposed to be run after the barcode threshold is set in the config file. This rules completes the pipeline. If the results from `determine_bc_cutoff` exists they will be used, otherwise all rules are carried out.

## `select_reads`

R scripts which filters the barcode table such that only the barcode combinations above the barcode cutoff are kept. The frequency table is small enough to keep in memory, but the barcode table might take large amounts of RAM. Therefore, the file with the barcode table is streamed and the filtered version ends up in a SQLite database which is later read by `count_genes`.


## `create_bc_plots`

R script whic creates the knee and density plots of the barcode distributions by the help of the `posDemux` package. Keep in mind that the density plot is scaled by the number of reads on the y-axis. 

## `bwa_index`
Prepares the reference genome for alignment by indexing the genome. This rule in invoked once per reference genome and the index files are created in the same directory as the reference genome.

## `bwa_align`
Performs the alignment of the reverse reads to the reference genome by `bwa mem`. The result of this step is a SAM file with the alignments.

## `remove_xs_tags`
The `XS` tag generated by `bwa mem` tells how many alternative alignments there are of a read. We need this information for `count_genes`, but `feature_counts` overwrites this tag causing confusion. Therefore, this rule renames the tag from `XS` to `XG`.

## `sam_to_bam`

Converts the SAM file to a BAM file as it is required by `feature_counts`.

## `bam_sort`

Sort the BAM file by index coordinate as required by `feature_counts`.

## `index_bam`

Effectively interating over the BAM files in the rules `feature_counts` and `count_genes` requires an BAM file index with the `.bai` extension. This rule ensures this index is being generated.

## `feature_counts`

Annotates the alignments with their corresponding features (i.e. genes, operons) using the reference annotation.

## `count_genes`

The final stage of the pipeline. It is served by a Python script which processes the data in three stages:

* Scans through the feature annotated alignments and compares with the barcode table. Any alignment which is not a primarily alignment of a read, any read not aligning to the reference genome, or any read for which no feature could be found, are discarded. Additionally, reads which align ambiguously, this is, with the highest alignment score to two or more genomic regions, are also discared with the expection of rRNA features for which ambiguously alignments are allowed. The read's barcode is looked up to see if it matches one of the selection reads from `select_reads`. For the reads which pass all these filters, the cell barcode, the UMI and gene annotation is recorded.
* Deduplicates barcode-gene groups, this is: For each combination of cell barcode and gene, the unique UMIs are extracted and counted. For this process [UMI-tools](https://github.com/CGATOxford/UMI-tools) `directional` algorithm is used for fault tolerance.
* Creates the gene count table and UMI count table which are the principal end products of the pipeline.


# Customizing the pipeline

When creating such a pipeline, there is always a compromise between making the pipeline flexible enough to accomodate new data and avoiding bloating the pipeline with options and features which are not really needed. This pipeline is therefore written for the existing PETRI-seq data in mind, but we still think the pipeline can be customized for similar assays. For this reason, we give a brief guide of the elements of the pipeline which can be adopted if needed:

## Parallelization

The rules `bwa_align`, `count_genes`, `preprocess` are the only ones where parallelism is used inside a single rule invocation (independent rule invocations can still run in parallel). Sometimes the `cores:` directive should be changed to allow better scheduling of the jobs. This is especially prominent for `count_genes` where only the second stage of the pipeline in parallelized, while the first step is usually the slowest one.

## Quality trimming

In the rule `quality_trim` in `rules/preprocess.smk`, the low quality bases with a quality below 10 are trimmed from both ends of the reads and reads with more than 3 N nucleteotides are removed. In addition, forward reads with less than 58 nucleotides and reverse reads with less than 14 nucleotides are discarded. If you want to change this behavior, please consult Cutadapt's reference guide (https://cutadapt.readthedocs.io/en/stable/reference.html).


## Barcodes

The pipeline is currently hard-coded to use the same barcodes as Blattman et al. 2020. The barcode files are located in `../resources/sc_barcodes_v2`. In order to change the barcode files being used, update the fields `bc1`, `bc2`, `bc3` for the rule `demultiplex` in `rules/sc_pipeline`. Note that if the barcodes are of different size than in the original setup (7 nt), the segment lengths must be updated. Also, if the location of the barcodes in the reads changes, this must also be accounted for by the specification of the segmentation.

## Segmentation

The PETRI-seq protocol mandates the following segmentation of the forward reads (from 5' to 3'):

*   UMI: 7 nucletides
*   Barcode 3(`bc3`) : 7 nucleotides
*   Linker: 15 nucletides
*   Barcode 2(`bc2`): 7 nucleotides
*   Linker: 14 nucletides
*   Barcode 1(`bc1`): 7 nucleotides
*   The rest of the read is ignored

For altering this sequence, you need to modify `scripts/demultiplexer.R`, by changing the variables `sequence_annotation`, `segment_lengths`, and `BARCODE_WIDTH` with help from the documentation for `posDemux::combinatorial_demultiplex()`.
