---
title: "scShardSplitRef: Algorithm Description"
author: "Irina Kuznetsova"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{scShardSplitRef: Algorithm Description}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
## 1. Purpose
This toy dataset provides a minimal example to demonstrate [the logic of the algorithm](#algorithm-walk-through) that prepares compatible reference files (GTF and FASTA) for `cellranger-arc mkref`.
It allows users to test the workflow using toy data and explains the algorithm logic in detail.
## 2. {scShardSplitRef} Guided Walkthrough Using a Toy Example
### Installation
```{r setup}
library("scShardSplitRef")
```
### Prerequisites
| Software | Version | Purpose |
|----------|---------|---------------------|
| bedtools | ≥ 2.30 | BED file operations |
### Input Files
| File | Format | Description | File Name |
|------------------|------------------|------------------|------------------|
| Gene annotations | GTF | Gene features and coordinates | extdata/AlgorithmToy.gtf |
| Centromere coordinates | BED-like | Centromere positions per chromosome | extdata/AlgorithmToy.bed |
### Workflow
#### 1. Determine Split Regions
The next function called `determine_split_regions()` prepares split coordinates in appropriate BED-like format.
This function ensures that chromosomes are not split in the middle of a gene region.
##### The Help Page for `determine_split_regions()`
```{r, eval=FALSE}
?determine_split_regions
```
**Input Data** This function takes a BED file and a GTF file as input.
- `BED`: 3 columns, tab-delimited, no header, no comments format (this centromere coordinates/`0` file is created manually)
- `GTF`: 9 columns, tab-delimited, lines beginning with `#` at the top of the file are treated as comments and allowed
- All files need to be decompressed
**Parameters:**
- `limit`: can be left as the default, which is the Cell Ranger ARC limit of `2^²⁹ (536,870,912 bp)`.
For the purpose of this tutorial, to better understand the algorithm, we set `limit=540`
- `shift_by`: the step size used to move the split coordinate to the right when searching for a safe location outside a genic region
- `clearance`: a buffer that checks whether there is a clear region free of genes on both sides of the proposed split coordinate
##### Usage Example: `determine_split_regions()`
In this example we use:
``` r
limit = 540L
shift_by = 30L
clearance = 40L
```
```{r}
# Read GTF file
gtf_path <- system.file(
"extdata/AlgorithmToy.gtf",
package = "scShardSplitRef"
)
# Read BED file
bed_path <- system.file(
"extdata/AlgorithmToy.bed",
package = "scShardSplitRef"
)
# Define output path
output_file <- file.path(tempdir(), "extdata")
dir.create(output_file, recursive = TRUE, showWarnings = FALSE)
output_file_name_dir <- file.path(
output_file,
"AlgorithmToy_DetermineSplitReg.bed"
)
# Determine split regions
get_split_reg <- determine_split_regions(
bed = bed_path,
gtf = gtf_path,
output_bed = output_file_name_dir,
limit = 540L,
shift_by = 30L,
clearance = 40L
)
get_split_reg
```
##### Algorithm Explanation: How `determine_split_regions()` works {#algorithm-walk-through}
#### No Overlap with Gene Coordinates
The chromosome is split at that coordinate, producing two regions:
- Left arm: from the beginning of the chromosome to the split coordinate
- Right arm: from the split coordinate to the end of the chromosome
> Note: If centromere locations are not provided (i.e., BED start position or 2nd column is `0`), the algorithm first creates provisional *even shards* based on `limit`.
> It then applies the same gene-clearance logic to each proposed split boundary, shifting boundaries as needed to avoid splitting through genic regions.
Algorithm Visual (Chr0): Centromere in BED File Defined
``` text
chr0: 450-800
Note, that each chromosome is defined by a coordinate range starting at 0. If a feature's start position is not 0, it indicates a split coordinate (centromere). In this example, 450 is the split coordinate (centromere), dividing the chromosome into two segments: 0–450 and 450–800.
limit = 540
shift_by = 30
clearance = 40
gene[G0]: 100-200
chr0: 0 100 200 800
|-----|====|---------|--------------|
[G0] ^
|
split (450)
|<---- Left arm ---->|<-- Right arm ->|
The centromere coordinate is defined in the BED file (chr0:450-800)
Split coordinate 450 is below the `limit` and does not overlap with any gene coordinates
Note: splitting first tries to use the centromere coordinate, not the `limit`!!!
Final: split at centromere coordinate 450 -> 1 split coordinate (2 shards)
```
Algorithm Visual (Chr1): Centromere in BED File Undefined
``` text
chr1: 0-820
Note: If centromere locations are NOT provided (i.e., BED start position is `0` for a chromosome), the algorithm first creates provisional even shards based on `limit`. It then applies the same gene-clearance logic to each proposed split boundary, shifting boundaries as needed to avoid splitting through genic regions.
limit = 540
shift_by = 30
clearance = 40
gene[G1]: 100-200
chr1: 0 100 200 410 820
|-----|====|---------|--------------|
[G1] ^
|
split (410)
|<---- Left arm ---->|<-- Right arm ->|
The centromere coordinate is not defined in BED file (chr1:0-820).
The `limit` set to 540.
The algorithm first creates provisional even shards based on `limit`:
ceiling(820/540) = 2 (shards) -> 820/2 = 410
Split coordinate 410 -> outside of G1[100-200], PASSING clearance = 40
Final: split at 410 -> 1 split coordinate (2 shards)
```
### Overlap with Gene Coordinates
In case of overlapping split coordinate with genic region, the split coordinate needs to be moved to a safe location that does not overlap any genes.
The algorithm proceeds as follows:
The split position is calculated from the limit parameter.
It looks to see if there is sufficient clearance between the split point and any neighbouring genes left or right.
If the clearance test does not pass the function uses `shift_by`, which essentially defines the size of the steps to iteratively take as it tries to walk away from any genes that are too close.
> Note: `clearancep` and `shift_by` operate independently.
> `clearance` determines whether a position is acceptable; `shift_by` determines how far to move when it is not.
Algorithm Visual (Chr2): the initial split coordinate is inside of a genic region
``` text
chr2: 0-1000
Note, that If centromere locations are not provided (i.e., BED start position is `0` for a chromosome), the algorithm first creates provisional *even shards* based on `limit`. It then applies the same gene-clearance logic to each proposed split boundary, shifting boundaries as needed to avoid splitting through genic regions
limit = 540
shift_by = 30
clearance = 40
gene[G2]: 450-550
chr2: 0 450 550 1000
|----------------|=======|--------------|
[G2]
The centromere coordinate is not defined in BED file (chr2:0-1000).
The `limit` is set to 540.
The algorithm first creates provisional even shards based on `limit`:
ceiling(1000/540) = 2 -> 1000/2 = 500
chr2: 0 450 550 1000
|----------------|=======|--------------|
[G2] ^
|
split (500)
Split coordinate 500 -> falls inside of G2[450-550]
Shift right by 30 bp -> 530 -> inside G2[450–550]
Shift right by 30 bp -> 560 -> outside G2[450–550], not passing clearance = 40
Shift right by 30 bp -> 590 -> outside G2[450–550], PASSING clearance = 40
Split coordinate = 590
chr2: 0 450 550 1000
|----------------|=======|--------------|
[G2] ^
|
split (590)
The distance between 0 and 590 exceeds the set `limit` of 540, find the next split coordinate:
ceiling(590/540) = 2 —> 590/2 = 295
Split coordinate = 295
chr2: 0 450 550 1000
|----------------|=======|--------------|
^ [G2] ^
| |
split (295) split (590)
Final: split at 295 and 590 -> 2 split coordinates (3 shards)
```
Algorithm Visual (Chr3): the initial split coordinate is inside of two overlapping genes
``` text
chr3: 0-1000
Note, that If centromere locations are not provided (i.e., BED start position is `0` for a chromosome), the algorithm first creates provisional *even shards* based on `limit`. It then applies the same gene-clearance logic to each proposed split boundary, shifting boundaries as needed to avoid splitting through genic regions
limit = 540
shift_by = 30
clearance = 40
gene[G3L]: 460-545
gene[G3R]: 490-560
chr3: 0 460 490 545 560 1000
|----------------|===|overlap|======|--------------|
[G3L] [G3R]
The centromere coordinate is not defined in BED file (chr3:0-1000)
The `limit` set to 540
The algorithm first creates provisional even shards based on `limit`:
ceiling(1000/540) = 2 -> 1000/2 = 500
chr3: 0 460 490 545 560 1000
|----------------|===|overlap|======|--------------|
[G3L] ^ [G3R]
|
split (500)
Split coordinate 500 -> falls inside of G3L[460-545] and G3R[490-560]
Shift to the right by 30bp -> 530 - inside of the G3R[490-560]
Shift to the right by 30bp -> 560 - still inside of the G3R[490-560]
Shift to the right by 30bp -> 590 - outside of the G3R[490-560], not passing clearance
Shift to the right by 30bp -> 620 - outside of the G3R[490-560] -> PASSING clearance
Split coordinate = 620
chr3: 0 460 490 545 560 1000
|----------------|===|overlap|======|--------------|
[G3L] [G3R] ^
|
split (620)
The distance between 0 and 620 exceeds the `limit` of 540, find the next split coordinate:
ceiling(620/540) = 2 —> 620/2 = 310 —> outside G3L[460-545] and G3R[490-560], PASSING clearance = 40
Split coordinate = 310
chr3: 0 460 490 545 560 1000
|----------------|===|overlap|======|--------------|
^ [G3L] [G3R] ^
| |
split (310) split (620)
Final: split at 620 and 310 -> 2 split coordinates (3 shards)
```
Algorithm Visual (Chr4): one gene in each arm (no collision)
``` text
chr4: 300-1500
Note, that each chromosome is defined by a coordinate range starting at 0. If a feature's start position is not 0, it indicates a split coordinate (centromere). In this example, 300 is the split coordinate (centromere), dividing the chromosome into two segments: 0–300 and 300–1500.
limit = 540
shift_by = 30
clearance = 40
gene[G4L]: 50-100
gene[G4R]: 800-900
chr4: 0 50 100 800 900 1500
|-----|===|--------------------------|======|------|
[G4L] ^ [G4R]
|
split (300)
The centromere coordinate is defined in the BED file (chr4:300-1500) ->
below the `limit` and does not overlap with any gene coordinates
Note: splitting first tries to use the centromere coordinate, not the `limit`!!!
The distance between 300 and 1500 exceeds the `limit` of 540, find how many shards:
ceiling( (1500-300)/540)=3 -> (1500-300)/3=400
Next split coordinate 300+400=700, outside of G4R[800-900], PASSING clearance = 40
Next split coordinate 700+400=1100, outside of G4R[800-900], PASSING clearance = 40
chr4: 0 50 100 800 900 1500
|-----|===|--------------------------|======|------|
[G4L] ^ ^ [G4R] ^
| | |
split (300) split (700) split (1100)
Final: split at 300, 700 and 1100 -> 3 split coordinates (4 shards)
```
Algorithm Visual (Chr5): a split coordinate cuts through a gene
``` text
chr5: 600-1500
# Note, that each chromosome is defined by a coordinate range starting at 0. If a feature's start position is not 0, it indicates a split coordinate. In this example, 600 is the split coordinate, dividing the chromosome into two segments: 0–600 and 600–1500.
limit = 540
shift_by = 30
clearance = 40
gene[G5BS]: 550-650
The centromere coordinate is defined in the BED file (chr5:600-1500)
Left arm: `[0, 600)`
Right arm: `[600, 1500)`
chr5: 0 550 650 1500
|-----------------|======|---------------------------------|
^ [G5BS]
|
split (600)
I. Left arm region
Note: splitting first tries to use the centromere coordinate, not the `limit`!!!
The distance between 0 and 600 exceeds the `limit` of 540, find how many shards:
ceiling((600-0)/540) = 2 -> 600/2=300
Split coordinate = 300
chr5: 0 550 650 1500
|-----------------|======|---------------------------------|
^ [G5BS]
|
split (300)
II. Right arm region
The right arm region (chr5:600-1500)
The distance between 600 and 1500 exceeds the `limit` of 540, find how many shards:
ceiling((1500-600)/540) = 2 -> (1500-600)/2=450
Split coordinates = 600 and 1050
chr5: 0 550 650 1500
|-----------------|======|---------------------------------|
^ ^ [G5BS] ^
| | |
split (300) split (600) split (1050)
Split coordinate 600 -> inside of the G5BS[550-650]
Shift to the right by 30bp -> 630 - inside of the G5BS[550-650]
Shift to the right by 30bp -> 660 - outside of the G5BS[550-650], not passing clearance
Shift to the right by 30bp -> 690 - outside of the G5BS[550-650] -> PASSING clearance
Split coordinate = 690
chr5: 0 550 650 1500
|-----------------|======|---------------------------------|
^ [G5BS] ^ ^
| | |
split (300) split (690) split (1050)
Final: split at 300, 690 and 1050 -> 3 split coordinates (4 shards)
```
Algorithm Visual (Chr6): gene structure with exons, CDS, and transcript: gene overlaps with cutting coordinate
``` text
chr6: 0-1000
limit = 540
shift_by = 30
clearance = 40
gene[G6R]: 400-500
transcript[G6R]: 400-500
exon[G6R]: 400-450
exon[G6R]: 460-480
exon[G6R]: 490–500
CDS[G6R]: 410-495
five_prime_utr[G6R]: 400-409
chr6: 0 400 500 1000
|-----------------|======|---------------------------------|
[G6R]^
|
split (500)
The centromere coordinate is not defined in BED file (chr6:0-1000)
The `limit` set to 540
The algorithm first creates provisional even shards based on `limit`:
ceiling(1000/540) = 2 -> 1000/2 = 500
Split coordinate 500 -> falls on edge of G6R[400-500]
Shift to the right by 30bp -> 530 - outside of the G6R[400-500], not passing clearance
Shift to the right by 30bp -> 560 - outside of the G6R[400-500] -> PASSING clearance
chr6: 0 400 500 1000
|-----------------|======|---------------------------------|
[G6R] ^
|
split (560)
The distance between 0 and 560 exceeds the `limit` of 540, find the next split coordinate:
ceiling(560/540) = 2 —> 560/2 = 280 —> outside G6R[400-500], PASSING clearance = 40
Split coordinate = 280
chr6: 0 400 500 1000
|-----------------|======|---------------------------------|
^ [G6R] ^
| |
split (280) split (560)
Final: split at 280, 560 -> 2 split coordinates (3 shards)
```
> Note: in real world scenarios users should consider picking `clearance = 900 (2 x 150 bp reads with up to a 600 bp insert size)`, and `shift_by in the region of 1000 to 4000` to ensure the algorithm can easily walk through a whole gene if needed.
#### 2. Process a GTF file by assigning split-region sequence names
Reads BED-like "split regions" and GTF files.
Assigns each GTF feature to exactly one split-region interval, and replaces the GTF Chr field with a split-region name of the form `Chr-RegionStart-RegionEnd`.
##### The Help Page for `process_gtf()`
```{r, eval=FALSE}
?process_gtf
```
**Input Data** This function takes a BED-like file containing split intervals as produced by `determine_split_regions()` and a GTF file as the input.
##### Usage Example: `process_gtf()`
```{r, eval=TRUE}
# BED
bed_modif <- system.file(
"extdata/AlgorithmToy_DetermineSplitReg.bed",
package = "scShardSplitRef"
)
# Function
gtf_processed <- process_gtf(
split_regions_bed = bed_modif,
gtf = gtf_path,
genome_name = "AlgorithmToyExample",
genome_version = "v1",
keep_attributes = NULL,
out_path = tempdir()
)
```
## 3. Acknowledgements
GRDC funded this project through the Analytics for the Australian Grains Industry (AAGI) investment, CUR2210-005OPX (AAGI-CU).
This work was supported by resources provided by the [Pawsey Supercomputing Research Centres](https://pawsey.org.au/) Nimbus Research Cloud (https://doi.org/10.48569/v0j3-qd51), with funding from the Australian Government and the Government of Western Australia.
This research used the Australian Research Data Commons (ARDC) ARDC Nectar Research Cloud Service [ARDC Nectar Research Cloud | Australian Research Data Commons](https://ardc.edu.au/services/ardc-nectar-research-cloud/). The ARDC is enabled by the Australian Government’s National Collaborative Research Infrastructure Strategy (NCRIS).