DropletQC

https://waltermuskovic.github.io/DropletQC/
Admin
DropletQC is an R package that calculates nuclear fraction scores from scRNA-seq BAM files to...
Tokens:6,598
Snippets:43
Trust Score:9.7
Update:5 months ago
Show doc for...
Tokens:
Raw
### Install DropletQC Package

Source: https://waltermuskovic.github.io/DropletQC/index.html

Install the DropletQC package from GitHub using devtools. Ensure devtools is installed first.

```R
# install.packages("devtools")
devtools::install_github("WalterMuskovic/DropletQC", build_vignettes = TRUE)
```

--------------------------------

### Example Usage of identify_empty_drops

Source: https://waltermuskovic.github.io/DropletQC/reference/identify_empty_drops.html

This example demonstrates how to use the `identify_empty_drops` function with sample data. It preprocesses the data, calls the function, and then displays the head of the resulting data frame and a table of cell statuses.

```R
data("qc_examples")
gbm <- qc_examples[qc_examples$sample=="GBM",]
gbm.ed <- gbm[,c("nuclear_fraction_droplet_qc","umi_count")]
gbm.ed <- identify_empty_drops(nf_umi = gbm.ed)
head(gbm.ed)

#>                    nuclear_fraction_droplet_qc umi_count cell_status
#> AAACCCAAGGCGATAC-1                   0.1947243      2226        cell
#> AAACCCAAGGCTGTAG-1                   0.2766798      1063        cell
#> AAACCCACAAGTCCCG-1                   0.1843824     17883        cell
#> AAACCCACAGATGCGA-1                   0.2919902      8172        cell
#> AAACCCACAGGTGAGT-1                   0.3295617      9057        cell
#> AAACCCAGTCTTGCGG-1                   0.3795893      5612        cell
table(gbm.ed$cell_status)

#> 
#>          cell empty_droplet 
#>          5296           220 


```

--------------------------------

### Example: Identify Damaged Cells

Source: https://waltermuskovic.github.io/DropletQC/reference/identify_damaged_cells.html

This example demonstrates how to use the identify_damaged_cells function. It first identifies empty droplets, then uses the output to identify damaged cells, and finally displays the head of the resulting data frame and a table of cell statuses.

```R
data("qc_examples")
gbm <- qc_examples[qc_examples$sample=="MB",]
gbm.ed <- gbm[,c("nuclear_fraction_droplet_qc","umi_count")]
gbm.ed <- identify_empty_drops(nf_umi = gbm.ed)
gbm.ed$cell_type <- gbm$cell_type
gbm.ed.dc <- identify_damaged_cells(gbm.ed, verbose=FALSE)
gbm.ed.dc <- gbm.ed.dc[[1]]
head(gbm.ed.dc)

#>                    nuclear_fraction_droplet_qc umi_count cell_status
#> AAACCCACAAGAATAC-1                   0.4176170     12557        cell
#> AAACCCACAATAAGGT-1                   0.3744720     10958        cell
#> AAACCCACAGCCCACA-1                   0.6492798      6095        cell
#> AAACCCACAGTAACGG-1                   0.3900973     12802        cell
#> AAACCCACATAAGATG-1                   0.1312698     13672        cell
#> AAACCCACATAATCGC-1                   0.4443030     26653        cell
#>                                cell_type
#> AAACCCACAAGAATAC-1   neuron_unresolved_2
#> AAACCCACAATAAGGT-1 migrating_interneuron
#> AAACCCACAGCCCACA-1    neuron_hippocampus
#> AAACCCACAGTAACGG-1 migrating_interneuron
#> AAACCCACATAAGATG-1   neuron_unresolved_1
#> AAACCCACATAATCGC-1 migrating_interneuron
table(gbm.ed.dc$cell_status)

#> 
#>          cell  damaged_cell empty_droplet 
#>          7389          1349           920 



```

--------------------------------

### Calculate Nuclear Fraction from Outs Directory

Source: https://waltermuskovic.github.io/DropletQC/reference/nuclear_fraction_tags.html

Calculates nuclear fraction using an 'outs' directory. Set tiles and cores to 1 for simplicity and verbose to FALSE for quiet execution.

```R
nf1 <- nuclear_fraction_tags(
    outs = system.file("extdata", "outs", package =
    "DropletQC"),
     tiles = 1, cores = 1, verbose = FALSE)
head(nf1)
```

--------------------------------

### Calculate Nuclear Fraction using Annotation

Source: https://waltermuskovic.github.io/DropletQC/index.html

Use nuclear_fraction_annotation for a more flexible calculation of the nuclear fraction score when BAM files lack region tags. Requires annotation, BAM, and barcode files.

```R
nf2 <- nuclear_fraction_annotation(
 annotation_path = system.file("extdata/outs/chr1.gff3",package = "DropletQC"),
 bam = system.file("extdata/outs/possorted_genome_bam.bam",package = "DropletQC"),
 barcodes = system.file("extdata/outs/filtered_feature_bc_matrix/barcodes.tsv.gz",package = "DropletQC"),
 tiles = 1, cores = 1, verbose = FALSE)
head(nf2)
```

--------------------------------

### Identify Empty Droplets Function Signature

Source: https://waltermuskovic.github.io/DropletQC/reference/identify_empty_drops.html

This is the function signature for `identify_empty_drops`. It takes nuclear fraction and UMI data, along with optional parameters for rescue thresholds and plotting.

```R
identify_empty_drops(
  nf_umi,
  nf_rescue = 0.05,
  umi_rescue = 1000,
  include_plot = FALSE,
  plot_name = NULL,
  plot_path = NULL,
  plot_width = 18,
  plot_height = 13,
  pdf_png = "png"
)
```

--------------------------------

### identify_empty_drops Function

Source: https://waltermuskovic.github.io/DropletQC/reference/identify_empty_drops.html

Identifies empty droplets by calculating the kernel density estimate of nuclear fraction scores and finding the trough after the first peak.

```APIDOC
## identify_empty_drops

### Description
This function is used to identify a suitable nuclear fraction cut-off point to guide the identification of empty droplets. To do this it calculates the kernel density estimate of the input nuclear fraction scores and identifies the trough after the first peak, which is assumed to represent the population of empty droplets.

### Arguments
- **nf_umi** (data frame): A data frame containing two columns; the nuclear fraction estimates in the first column and the total UMI count for each barcode in the second column.
- **nf_rescue** (numeric, optional): A rescue parameter defining a minimum nuclear fraction score between zero and one. This is used in combination with `umi_rescue` to identify cells that were misidentified as empty droplets. Defaults to 0.05.
- **umi_rescue** (integer, optional): A rescue parameter defining a minimum UMI count. This is used in combination with `nf_rescue` to identify cells that were misidentified as empty droplets. Defaults to 1000.
- **include_plot** (logical, optional): Whether or not to produce a plot illustrating how the nuclear fraction threshold was identified and which barcodes have been called as empty droplets. In the plot of nuclear fraction vs log10(UMI counts), empty droplets are expected to occupy the lower left corner of the plot. Defaults to FALSE.
- **plot_name** (character, optional): If provided, a plot will be saved with the provided name.
- **plot_path** (character, optional): If provided, a plot will be saved to the provided path.
- **plot_width** (numeric, optional): Plot width in cm. Defaults to 18.
- **plot_height** (numeric, optional): Plot height in cm. Defaults to 13.
- **pdf_png** (character, optional): Either "png" or "pdf". Defaults to "png".

### Value
data frame: The original data frame is returned plus an additional column identifying each barcode as a "cell" or "empty_droplet".

### Examples
```R
data("qc_examples")
gbm <- qc_examples[qc_examples$sample=="GBM",]
gbm.ed <- gbm[,c("nuclear_fraction_droplet_qc","umi_count")]
gbm.ed <- identify_empty_drops(nf_umi = gbm.ed)
head(gbm.ed)

#\>                    nuclear_fraction_droplet_qc umi_count cell_status
#\> AAACCCAAGGCGATAC-1                   0.1947243      2226        cell
#\> AAACCCAAGGCTGTAG-1                   0.2766798      1063        cell
#\> AAACCCACAAGTCCCG-1                   0.1843824     17883        cell
#\> AAACCCACAGATGCGA-1                   0.2919902      8172        cell
#\> AAACCCACAGGTGAGT-1                   0.3295617      9057        cell
#\> AAACCCAGTCTTGCGG-1                   0.3795893      5612        cell
table(gbm.ed$cell_status)

#\> 
#\>          cell empty_droplet 
#\>          5296           220 

```
```

--------------------------------

### Calculate Nuclear Fraction from BAM and Barcodes

Source: https://waltermuskovic.github.io/DropletQC/reference/nuclear_fraction_tags.html

Calculates nuclear fraction by directly providing BAM file path and a vector of barcodes. Set tiles and cores to 1 for simplicity and verbose to FALSE for quiet execution.

```R
nf2 <- nuclear_fraction_tags(
   bam = system.file("extdata", "outs","possorted_genome_bam.bam", package =
   "DropletQC"),
   barcodes = c("AAAAGTCACTTACTTG-1",
                "AAAAGTGGATCTCTAA-1",
                "AAACACGTTCTCATCG-1"),
   tiles = 1, cores = 1,
   verbose = FALSE)
nf2
```

--------------------------------

### Calculate Nuclear Fraction using Tags

Source: https://waltermuskovic.github.io/DropletQC/index.html

Use nuclear_fraction_tags to efficiently calculate the nuclear fraction score when BAM files contain region tags. Point to the output directory of Cell Ranger.

```R
library(DropletQC)
nf1 <- nuclear_fraction_tags(
    outs = system.file("extdata", "outs", package = "DropletQC"),
     tiles = 1, cores = 1, verbose = FALSE)
head(nf1)
```

--------------------------------

### identify_damaged_cells

Source: https://waltermuskovic.github.io/DropletQC/reference/identify_damaged_cells.html

Identifies damaged cells by analyzing UMI counts and nuclear fraction. It uses a Gaussian mixture model to differentiate between 'cell' and 'damaged_cell' populations.

```APIDOC
## identify_damaged_cells

### Description
This function uses a combination of the cell UMI counts and the nuclear fraction score to assign each cell one of two values; "cell" or "damaged_cell". This is based on the idea that damaged cells have a lower UMI count and higher nuclear fraction than whole cells. The expected input is a data frame with four columns. The first three columns should contain; the nuclear fraction score, total UMIs and a character vector describing each cell as "cell" or "empty_droplet". This is the format output by the `identify_empty_drops` function. The fourth column should be a character vector with user-assigned cell types. Internally, the provided data frame is split by cell type and a Gaussian mixture model with a maximum of two components is fit to the umi counts and nuclear fraction scores. The parameters of the model are estimated using expectation maximisation (EM) with the `mclust` package. The best model is selected using the Bayesian Information Criterion (BIC). The two populations (cells and damaged cells) are assumed to have equal variance (mclust model name "EEI").

### Arguments
- **nf_umi_ed_ct** (data frame) - Required - A data frame with four columns. The first three columns should match the output from the `identify_empty_drops` function. The fourth column should contain cell type names.
- **nf_sep** (numeric) - Optional - The minimum separation of the nuclear fraction score required between the cell and damage cell populations.
- **umi_sep_perc** (numeric) - Optional - This is the minimum percentage of UMIs which the damaged cell population is required to have compared to the cell population. For example, if the mean UMI of the distribution fit to the whole cell population is 10,000 UMIs, the mean of the distribution fit to the damaged cell population must be at less than 7,000 UMIs if the umi_sep parameter is 30 (%)
- **output_plots** (logical) - Optional - Whether or not to return plots.
- **verbose** (logical) - Optional - Whether to print updates and progress while fitting with EM.

### Value
list - A list of length two. The first element in the list contains a data frame with the same dimensions input to the `nf_umi_ed_ct` argument, with "damaged_cell" now recorded in the third column. The second element is NULL unless `output_plots`=TRUE. If requested, three plots are returned for each cell type in a named list, combined using `ggpubr::ggarrange`. For each cell type, the first plot illustrates the cell and damaged cell populations (if any) in a plot of nuclear fraction vs log10(UMI counts). Damaged cells are expected to be in the lower right portion of the plot(lower UMI counts and higher nuclear fraction). The second and third plots show the model fits to the nuclear fraction and UMI count distributions respectively. Solid lines indicate the distribution mean, while dashed lines indicate the positions of the thrsholds controlled by the `nf_sep` and `umi_sep` parameters.

### Example
```R
#1
data("qc_examples")
gbm <- qc_examples[qc_examples$sample=="MB",]
gbm.ed <- gbm[,c("nuclear_fraction_droplet_qc","umi_count")]
gbm.ed <- identify_empty_drops(nf_umi = gbm.ed)
gbm.ed$cell_type <- gbm$cell_type
gbm.ed.dc <- identify_damaged_cells(gbm.ed, verbose=FALSE)
gbm.ed.dc <- gbm.ed.dc[[1]]
head(gbm.ed.dc)

#>                    nuclear_fraction_droplet_qc umi_count cell_status
#> AAACCCACAAGAATAC-1                   0.4176170     12557        cell
#> AAACCCACAATAAGGT-1                   0.3744720     10958        cell
#> AAACCCACAGCCCACA-1                   0.6492798      6095        cell
#> AAACCCACAGTAACGG-1                   0.3900973     12802        cell
#> AAACCCACATAAGATG-1                   0.1312698     13672        cell
#> AAACCCACATAATCGC-1                   0.4443030     26653        cell
#>                                cell_type
#> AAACCCACAAGAATAC-1   neuron_unresolved_2
#> AAACCCACAATAAGGT-1 migrating_interneuron
#> AAACCCACAGCCCACA-1    neuron_hippocampus
#> AAACCCACAGTAACGG-1 migrating_interneuron
#> AAACCCACATAAGATG-1   neuron_unresolved_1
#> AAACCCACATAATCGC-1 migrating_interneuron
table(gbm.ed.dc$cell_status)

#> 
#>          cell  damaged_cell empty_droplet 
#>          7389          1349           920 


```
```

--------------------------------

### nuclear_fraction_annotation

Source: https://waltermuskovic.github.io/DropletQC/reference/nuclear_fraction_annotation.html

Calculates the nuclear fraction score by parsing BAM file reads and determining overlap with intron/exon intervals from an annotation file.

```APIDOC
## nuclear_fraction_annotation

### Description
Calculates the nuclear fraction score by parsing all of the reads in the provided BAM file and determining their overlap with intron and exon intervals obtained from the provided annotation file.

### Arguments
- **annotation_path** (character) - Path to the annotation file.
- **annotation_format** (character, optional) - Format of the annotation file. Can be "auto", "gff3", or "gtf". Defaults to "auto".
- **bam** (character) - Path to the BAM file.
- **bam_index** (character, optional) - Path to the BAM file index. Defaults to `paste0(bam, ".bai")`.
- **barcodes** (character) - Either a vector of barcode names or the path to a `barcodes.tsv.gz` file.
- **cell_barcode_tag** (character, optional) - BAM tag containing the cell barcode. Defaults to "CB".
- **cores** (numeric, optional) - Number of cores for parallel processing. Defaults to `future::availableCores() - 1`.
- **tiles** (integer, optional) - Number of tiles to split transcripts into for processing chunks. Defaults to 1000.
- **verbose** (logical, optional) - Whether to print progress. Defaults to TRUE.

### Value
data.frame - A data frame containing the nuclear fraction score. The row names match the supplied barcodes.

### Example
```R
nf3 <- nuclear_fraction_annotation(
  annotation_path = system.file("extdata/outs/chr1.gff3",
    package = "DropletQC"),
  bam = system.file("extdata/outs/possorted_genome_bam.bam",
    package = "DropletQC"),
  barcodes = system.file(
    "extdata/outs/filtered_feature_bc_matrix/barcodes.tsv.gz",
    package = "DropletQC"),
  tiles = 1, cores = 1, verbose = FALSE
)
head(nf3)
```
```

--------------------------------

### Calculate Nuclear Fraction Annotation

Source: https://waltermuskovic.github.io/DropletQC/reference/nuclear_fraction_annotation.html

Use this function to calculate the nuclear fraction score. It requires paths to annotation and BAM files, and optionally cell barcodes. Adjust cores, tiles, and verbosity for performance and output.

```R
nuclear_fraction_annotation(
  annotation_path,
  annotation_format = "auto",
  bam,
  bam_index = paste0(bam, ".bai"),
  barcodes,
  cell_barcode_tag = "CB",
  cores = future::availableCores() - 1,
  tiles = 1000,
  verbose = TRUE
)
```

```R
nf3 <- nuclear_fraction_annotation(
 annotation_path = system.file("extdata/outs/chr1.gff3",
  package = "DropletQC"),
 bam = system.file("extdata/outs/possorted_genome_bam.bam",
  package = "DropletQC"),
 barcodes = system.file(
 "extdata/outs/filtered_feature_bc_matrix/barcodes.tsv.gz",
  package = "DropletQC"),
 tiles = 1, cores = 1, verbose = FALSE)
head(nf3)
```

--------------------------------

### nuclear_fraction_tags Function Definition

Source: https://waltermuskovic.github.io/DropletQC/reference/nuclear_fraction_tags.html

Defines the nuclear_fraction_tags function with parameters for input data, BAM tags, and processing options.

```R
nuclear_fraction_tags(
  outs = NULL,
  bam = NULL,
  bam_index = paste0(bam, ".bai"),
  barcodes = NULL,
  cores = future::availableCores() - 1,
  tiles = 100,
  cell_barcode_tag = "CB",
  region_type_tag = "RE",
  exon_tag = "E",
  intron_tag = "N",
  verbose = TRUE
)
```

--------------------------------

### Identify Damaged Cells Function Call

Source: https://waltermuskovic.github.io/DropletQC/reference/identify_damaged_cells.html

Call the identify_damaged_cells function with specified parameters. Ensure the input data frame `nf_umi_ed_ct` has the correct format.

```R
identify_damaged_cells(
  nf_umi_ed_ct,
  nf_sep = 0.15,
  umi_sep_perc = 50,
  output_plots = FALSE,
  verbose = TRUE
)
```

--------------------------------

### nuclear_fraction_tags Function

Source: https://waltermuskovic.github.io/DropletQC/reference/nuclear_fraction_tags.html

Calculates the nuclear fraction statistic for each input cell barcode. It uses region type tags in a BAM file to determine the fraction of reads that are intronic.

```APIDOC
## nuclear_fraction_tags

### Description
Uses region type tags in a provided BAM file to calculate for each input cell barcode the nuclear fraction statistic. This is the fraction of reads that are intronic: nuclear fraction = # intronic reads / (# intronic reads + # of exonic reads).

### Parameters
#### Path Parameters
- **outs** (character) - Required - The path to the 'outs' directory created by Cell Ranger. If not provided, `bam`, `bam_index`, and `barcodes` must be supplied.
- **bam** (character) - Optional - The path to the input bam file. Not required if an 'outs' directory is provided.
- **bam_index** (character) - Optional - The path to the input bam file index. Defaults to `paste0(bam, ".bai")`. Not required if an 'outs' directory is provided.
- **barcodes** (character) - Optional - Either a vector of barcode names or the path to the barcodes.tsv.gz file output by Cell Ranger. Not required if an 'outs' directory is provided.

#### Query Parameters
- **cores** (numeric) - Optional - Runs the function in parallel using furrr:future_map() with the requested number of cores. Defaults to `future::availableCores() - 1`.
- **tiles** (numeric) - Optional - Splits the genome into tiles to speed up BAM file processing. Defaults to 100.
- **cell_barcode_tag** (character) - Optional - The BAM tag containing the cell barcode sequence. Defaults to "CB".
- **region_type_tag** (character) - Optional - The BAM tag containing the region type. Defaults to "RE".
- **exon_tag** (character) - Optional - The character string that defines a read as exonic. Defaults to "E".
- **intron_tag** (character) - Optional - The character string that defines a read as intronic. Defaults to "N".
- **verbose** (logical) - Optional - Whether or not to print progress. Defaults to TRUE.

### Returns
data.frame - A 1-column data frame containing the calculated nuclear fraction statistic for each input barcode. The order and names of the rows will match those of the input cell barcodes.

### Examples
```R
# Example using outs directory
nf1 <- nuclear_fraction_tags(
    outs = system.file("extdata", "outs", package = "DropletQC"),
     tiles = 1, cores = 1, verbose = FALSE)
head(nf1)

# Example using bam and barcodes directly
nf2 <- nuclear_fraction_tags(
   bam = system.file("extdata", "outs", "possorted_genome_bam.bam", package = "DropletQC"),
   barcodes = c("AAAAGTCACTTACTTG-1",
                "AAAAGTGGATCTCTAA-1",
                "AAACACGTTCTCATCG-1"),
   tiles = 1, cores = 1,
   verbose = FALSE)
nf2
```
```
Tokens:
Raw