isopod

KO | EN

gitlite — search

#python #java #python3 #arduino #golang #machine-learning #rust #html #flask #javascript #seismology #nodejs

★ 11 Open GitHub ↗

An R package to simulate and detect differential transcript usage in transcript-level single-cell datasets.

// repository documentation

Was this content helpful?

★ 0 (0 ratings)

isopod
================

# Overview

isopod is an R package with two distinct purposes centered around
long-read single-cell data:

1.  Simulation of long-read single-cell datasets with
    gene/transcript-level diifferential transcript usage (DTU) at the
    counts table level.

2.  Identification of genes and transcripts exhibiting statistically
    significant DTU in real long-read single-cell datasets.

Isopod’s DTU detection requires two main inputs, while isopod’s
simulation component creates these as outputs:

- An isoform-level counts table as a dataframe, with isoform names, gene
  names, and counts data.
- A dataframe that specifies which group each cell belongs to.

When performing DTU analysis, provided these inputs, isopod first
filters the counts table to exclude low count or sparsely represented
genes. In addition, low-count isoform counts are collapsed (but not
removed) into an artificial isoform on a per-gene basis, allowing us to
remove from analysis isoforms resulting from noisy read assignment,
while retaining their counts data for gene-level calculations required
by other isoforms in the gene.

The filtered dataset is then passed to a permutation analysis, which
finds transcripts exhibiting DTU between the specified cell group
against the rest of the dataset, in a one-vs-all fashion. If you want to
compare two specific cell groups, or to run each cell group against the
rest of the dataset automatically, check both [2.
get_permutation_pvals](#get_permutation_pvals) and the function
documentation at `get_permutation_pvals()`.

The final output from the pipeline is a list containing different
dataframes, one of which is information on whether isoforms exhibit DTU
in a cell group.

Isopod’s simulation component instead works to simulate the counts table
and cell group assignments given certain parameters by the user, such as
the gamma distribution parameter to draw mean gene counts from or the
number of cell groups. Users are able to control multiple generation
parameters, including:

- The number of genes, transcripts, cell groups and their sizes

- Gamma distribution parameters to draw average gene counts from

- Distribution of transcripts among genes

- Proportion of genes to exhibit DTU, along with which transcripts
  should show DTU

- Proportion of genes to exhibit DGE, along with the magnitude of DGE

- Proportion of cells and genes to exhibit outlier counts, along with
  the outlier magnitude

- Counts generation from either a Poisson or negative binomial
  distribution.

Isopod then uses these parameters to create a simulated dataset,
outputting a list that contains the main output tables. The list also
keeps information on all input parameters, DTU/DGE/outlier genes, and
other generation details for posterity, allowing for easier
reconstruction of similar datasets.

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

# Installing Isopod

The development version of Isopod can be installed and loaded from this
repository through the following code:

``` r
remotes::install_github('michael-nakai/isopod')
library(isopod)
```

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

# Simulating a long-read single-cell dataset

Isopod’s transcript-level counts table simulation allows users to create
a long-read single-cell dataset to their specifications. Here we’ll
generate a sample dataset with the following parameters:

- 3 cells groups, with the groups containing 100, 200, and 300 cells
  respectively.

- 1,000 genes, half containing two transcripts per gene and half
  containing three transcripts per gene.

- The default distribution parameters for gene counts and counts
  generation.

- 10% of genes exhibiting DTU.

- 10% of genes exhibiting DGE.

- 1% of genes and cells should exhibit outlier counts.

We first create a dataframe that defines the transcripts per gene
distribution.

``` r
transcripts_per_gene <- data.frame('isoforms_in_gene' = c(2, 3), 'proportion' = c(0.5, 0.5))
```

We can then immediately simulate the dataset in one line of code, using
the wrapper function `create_simulated_dataset()`.

``` r
dataset_name <- 'test_simulation'

simulated_dataset <- create_simulated_dataset(dataset_name,
 number_of_genes = 1000,
 cells_in_groups = c(100, 200, 300),
 isoform_number_distribution = transcripts_per_gene,
 proportion_of_genes_with_DTU = 0.1,
 proportion_of_genes_with_DGE = 0.1,
 outlier_gene_proportion = 0.01,
 proportion_of_cells_with_inflated_counts = 0.01)
```

## Outputs

The function produces a list containing all relevant data from the
simulation. The important outputs are:

- `cell_designations`: A dataframe tying cells to cell groups. DTU and
  DGE are generated between these groups.

- `counts_table`: The long-read single-cell transcript-level counts
  table.

- `isoform_means_table`: A dataframe of average isoform counts for every
  gene. Average gene counts are randomly drawn from a gamma
  distribution. Isoforms receive a proportion of the average gene
  counts, as seen in `isoform_props_table`.

- `cell_outliers` and `gene_expression_outliers`: Dataframes describing
  outlier cells and genes respectively, and the inflation factor applied
  to the cell or gene counts.

- `swapped_genes`: A list of vectors describing the genes exhibiting DTU
  in every subsequent cell group. DTU is mainly simulated against the
  first provided group of cells.

- `DGE_details`: A dataframe of genes exhibiting the genes exhibiting
  DGE, the group they exhibit them in, and the inflation factor assigned
  to every gene.

Other details regarding simulation parameters and inputs are also stored
in the output list for easy query of the simulation.

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

# Running isopod’s DTU detection on a dataset

To run isopod’s DTU detection, we need two dataframes:

- An isoform-level counts table with isoform and gene names
- A dataframe that ties each cell to a cell group.

Here’s an example of what the tables might look like:

``` r
counts_table <- data.frame('transcript_id' = c('transcript1', 'transcript2', 'transcript3', 'transcript4', 'transcript5'),
 'gene_id' = c('gene1', 'gene1', 'gene1', 'gene2', 'gene2'),
 'cell_1' = c(5, 2, 0, 10, 5),
 'cell_2' = c(5, 2, 0, 10, 5),
 'cell_3' = c(0, 2, 5, 5, 10),
 'cell_4' = c(0, 2, 5, 5, 10))
cell_groups <- data.frame('cell_id' = c('cell_1', 'cell_2', 'cell_3', 'cell_4'),
 'cell_group' = c('group_1', 'group_1', 'group_2', 'group_2'))

counts_table
#>   transcript_id gene_id cell_1 cell_2 cell_3 cell_4
#> 1   transcript1   gene1      5      5      0      0
#> 2   transcript2   gene1      2      2      2      2
#> 3   transcript3   gene1      0      0      5      5
#> 4   transcript4   gene2     10     10      5      5
#> 5   transcript5   gene2      5      5     10     10

cell_groups
#>   cell_id cell_group
#> 1  cell_1    group_1
#> 2  cell_2    group_1
#> 3  cell_3    group_2
#> 4  cell_4    group_2
```

To quickly run the pipeline and get the results of the permutation
analysis, we can run the following:

``` r
permutation_results <- run_everything(counts_table, 
 cell_groups, 
 transcript_id_colname = 'transcript_id', 
 gene_id_colname = 'gene_id',
 cell_labels_colname = 'cell_group', 
 cell_group_to_analyse = 'group_1', 
 output_folder = 'permutation_results')
```

The permutation analysis is computationally expensive, and may take some
time depending on the size of your dataset and the computational power
availble to you. Luckily, it’s also multithreaded, and will try to use
all available cores (except for one) on your computer or cluster to
speed up the permutations. This means that it can take anywhere from 5
minutes (for a small dataset) to a few hours (for a large dataset with
moderate compute power). It’s probably best to check the ETA listed on
the progress bar and have a coffee or dinner while isopod runs.

## Outputs

`run_everything()` creates a few files containing the permutation
results in the output folder. These files can be read back into R with
`readRDS()`, and are:

- `filtered_counts_table.rds` - A list containing the counts table after
  the filtering step, and a list of collapsed isoforms for each gene.
  See [1. filter_counts_table](#filter_counts_table).
- `permutation_results.rds` - A list of all results from the DTU
  permutation analysis.

The permutation results object contains a few items. The main outputs
are:

- `permutation_pvalues` - The isoform-level permutation *p*-values for
  the group specified. Isoforms exhibiting DTU between the group and all
  other cells will be significant.

- `permutation_pvalues_gene` - The gene-level permutation *p*-values for
  the group. This tests whether any change between the group and all
  other groups can be seen for each gene. Put simply, this tests whether
  any isoform proportions are significantly changing within the gene,
  between the specified group and all other cells.

- `unadjusted_permutation_pvalues` and
  `unadjusted_permutation_pvalues_gene`, which contain the raw
  *p*-values at the isoform and gene-level.

- `genes_filtered_from_cutoff` - A list of genes filtered during the
  permutation analysis. See [Permutation
  analysis](#permutation-analysis) for more information on the cutoff.

You’re probably looking for `permutation_pvalues` and
`permutation_pvalues_gene`, which gives you the permutation *p*-values
at the isoform and gene-level. These *p*-values have been corrected
using the Benjamini-Hochberg adjustment method. If you’d like raw
*p*-values instead to correct at your leisure, look at
`unadjusted_permutation_pvalues` and
`unadjusted_permutation_pvalues_gene` instead.

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

# Running individual functions for DTU detection

Running each function lets you access additional options for each, such
as running **all** groups at once instead of one at a time. Each
function invoked in `run_everything()` is listed below, with
explanations of all additonal options.

## 1. filter_counts_table

- `gene_count_threshold` : A gene-level filtering step. If the gene
  counts (found by taking the sum of all transcripts in a gene) are
  below this number, the transcripts associated with the gene will be
  filtered out of the table.

- `autofiltering` : A gene-level filtering step that can be set to
  `TRUE` or `FALSE`. Removes genes where at least N cells per group do
  not contain counts for the gene, where N is a proportion of the cell
  count of the smallest group (see `autofilter_threshold`). Put simply,
  filters genes that are not represented in many cells.

- `autofilter_threshold` : A decimal between 0 and 1 that determines the
  proportion of cells in the smallest group that must have a count
  during autofiltering.

- `cell_labels_table` and `cell_labels_colname` : Must be supplied if
  `autofiltering` is set. A table containing the cell IDs, and which
  cell group each ID belongs to. The `cell_labels_colname` is the name
  of the column that contains the cell groups.

- `collapse_isoforms_with_counts_below` : A transcript-level filtering
  step. Transcripts are removed from further analysis if they have total
  counts below this number. If two or more transcripts in a gene are
  removed from analysis, their counts are collapsed into an artificial
  isoform and the transcripts are removed from the counts table. This
  preserves the transcript counts (which is important during permutation
  analysis) while removing counts resulting from noisy or misaligned
  reads. Transcript IDs that are removed from further analysis are
  reported in the resulting object.

## 2. get_permutation_pvals

To accomodate and make easier different types of comparisons in the
dataset, isopod offers a few different permutation analysis run modes.

### Run modes

The permutation function has three distinct run modes, outlined briefly
below.

1.  One group against the rest of the dataset
2.  One group against another specific group
3.  All groups against the rest of the dataset

We generally recommend running the cell group of interest against either
the rest of the dataset (**run mode 1**) or another cell group (**run
mode 2**), as calculating permutations for all clusters at once is
computationally expensive and can take some time without access to a
cluster computing environment. However, we provide the option to
calculate permutation *p*-values for all clusters (**run mode 3**) if
needed, which can help with certain other analyses.

If you’d like to run isopod using **run modes 1 or 2**, then set the
following function arguments:

- `analysis_group_1` : This should be the name of the cell group of
  interest. If you’d like to compare one group against the rest of the
  dataset (**run mode 1)**, just set this argument and ignore
  `analysis_group_2`.

- `analysis_group_2` : If you’d like to compare one group directly
  against another (**run mode 2**), set this to the name of the second
  cell group of interest.

For **run mode 3,** set the following function argument:

- `run_on_all_groups` : Set this to `TRUE` to calculate permuation
  *p*-values for all clusters (**run mode 3**), and leave the
  `analysis_group_1` and `analysis_group_2` arguments as `NA`.

### Other permutation options

The permutation function also has some other options, which are
explained below:

- `cores` : The number of CPU cores made available to isopod. Generally,
  it’s OK to leave this as the default `0`, since this will let isopod
  try to use all available cores (except for one, which it leaves for
  you). The more cores available, the faster each permutation as
  *p*-value calculations and (if using run mode 3) per-group
  calculations are multithreaded. For very large datasets, running
  isopod on a cluster computing environment is recommended.

- `do_gene_level_comparisons` : Setting this to `TRUE` allows isopod to
  additionally calculate gene-level permutation *p*-values. A
  significant gene-level *p*-value indicates that an isoform is changing
  proportions between groups within the gene, but the specific isoform
  exhibiting DTU is not identified from this *p*-value. It is
  recommended to always calculate gene-level *p*-values, as this
  increases the capabilities of the `cutoff` (explained in [Permutation
  analysis](#permutation-analysis)) without significantly increasing
  runtime.

- `return_detailed_pvalue_tables` : Setting this to `TRUE` causes a few
  additional lists to be returned in the output. These lists contain the
  raw chi-squared p-value for each isoform (and gene if gene-level
  calculations are being run) for each permutation. This is useful for
  understanding how the permutation *p*-value was calculated for each
  isoform, since you can create a histogram of all *p*-values across all
  permutations for an isoform here. As a word of warning, this causes
  the object returned from `get_permutation_pvals()` to be quite large
  (easily 1-30GB depending on the dataset size, number of cell groups,
  and number of permutations), so beware if you don’t have access to
  much disk space or RAM. **We generally recommend keeping this option
  disabled**.

- `report_adjusted_pvalues` : Setting this to `FALSE` stops isopod from
  Benjamini-Hochberg adjusting the original chi-squared and permutation
  *p*-values. We recommend leaving this to the default `TRUE`, since raw
  *p*-values are always reported alongside the adjusted *p*-values
  anyway, but the option is included if you’re certain you’d like to
  adjust the *p*-values using a different adjustment algorithm.

- `cutoff` : A decimal between `0` and `1`, with a default value of
  `0.1`. See [Permutation analysis](#permutation-analysis) for details.

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

# DTU Method Details

Further details on the permutation analysis are provided below.

### Filtering the counts table

The filtering step will always apply the filtering conditions in the
following order:

1.  Genes with one associated transcript are removed (always enabled).
2.  Genes below the gene count threshold are removed.
3.  Genes without sufficient appearances across cells are removed.
4.  Transcripts with low counts are aggregated into a collapsed isoform.

### Permutation analysis

For each permutation, isopod tests whether the proportion of each
transcript within its gene changes significantly between the cell group
and all other cell groups by creating a per-isoform contingency table.
This contingency table is passed to a chi-squared calculation, and the
resulting p-value is recorded and stored for later use. Between
permutations, cells are shuffled randomly between all cell groups. Cell
groups therefore retain their size but have randomly-allocated cells
within. The final permutation p-value is a function of the number of
permutations where the calculated p-value was equal to or more
significant than the original chi-squared p-value. The more permutations
that were as or more significant than the original, the less significant
the p-value.

By default, a *p*-value cutoff is applied to the first permutation. What
this means is that isopod checks one or two *p*-values after the first
loop to determine whether permutation analysis should continue on the
transcript or gene. Specifically:

1.  If gene-level testing is being performed, genes with a first-loop
    chi-squared *p*-value equal to or above the cutoff will be filtered
    out, and permutation analysis will not be performed on isoforms
    belonging to these genes. If gene-level testing is not being
    performed, this step is ignored.
2.  Isoforms with a first-loop chi-squared *p*-value equal to or above
    this cutoff will be ignored in subsequent permutations. Their counts
    will remain to accurately calculate pseudo-bulked contingency tables
    for other isoforms in the gene, but permutations wil not be
    calculated for them.

This filtering of the isoforms and genes passed to subsequent
permutations after the first loop serves to reduce unnecessary
calculations, and often greatly speeds up permutations without
sacrificing DTU detections. Nevertheless, the cutoff can be disabled
through the argument `cutoff = 1` if required.

isopod

Explore Similar Repositories

Was this content helpful?

Recent Feedback

Download README