The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

scShardSplitRef: Algorithm Description

Irina Kuznetsova

1. Purpose

This toy dataset provides a minimal example to demonstrate the logic of the algorithm 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

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()
?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:

limit = 540L 
shift_by = 30L 
clearance = 40L
# 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
)
#> Preparing split regions from 'bed' and 'gtf' into 'output_bed' (`limit` = 540,
#> `shift_by` = 30, `clearance` = 40).
#> Built 21 split intervals across 7 chromosomes and wrote them to 'output_bed'.

get_split_reg
#> [1] "/tmp/RtmpkhUpln/extdata/AlgorithmToy_DetermineSplitReg.bed"


Algorithm Explanation: How determine_split_regions() works

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 
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 
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 
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 
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) 
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 
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 
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()
?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()
# 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()
)
#> ℹ Validating inputs and reading files...
#> All files loaded and validated successfully.
#> ℹ Finished writing processed GTF: /tmp/RtmpkhUpln/B1_FINAL_MODIFIED_GTF_AlgorithmToyExample_v1.gtf

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 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. The ARDC is enabled by the Australian Government’s National Collaborative Research Infrastructure Strategy (NCRIS).

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.