scMitoMut 1.0.0

Install package:

`BiocManager::install("scMitoMut")`

Load packages:

```
library(scMitoMut)
library(data.table)
library(ggplot2)
library(rhdf5)
```

Key functions used in the scMitoMut package:

`parse_mgatk()`

: Parses the mgatk output and saves the result in an H5 file.`open_h5_file()`

: Opens the H5 file and returns a “mtmutObj” object.`subset_cell()`

: Subsets the cells in the mtmutObj object.`run_model_fit()`

: Runs the model fitting and saves the results in the H5 file.`mut_filter()`

: Filters the mutations based on specific criteria.`plot_heatmap()`

: Plots a heatmap with different types of data, including p-values, allele frequencies, or binary mutation status.`export_df()`

,`export_af()`

,`export_pval()`

, and`export_binary()`

: Export the mutation data in various formats.

- Somatic mutation: a mutation that occurs in a cell or a cell population, but not in the germline.
- AF or VAF: allele frequency or variant allele frequency. It is the ratio of the number of reads supporting the variant allele to the total number of reads.
- p-value: the probability of observing a test statistic at least as extreme as the one that was observed, assuming that the null hypothesis is correct.

**IMPORTANT**: In this vignette, I will use the term “mutation” to refer to the somatic mutation. For the somatic mutation, I use the dominant allele detected as the reference allele. And if the reference allele frequency is significant (FDR < 0.05) lowe, I will call the **locus** a mutation.

Single-cell genomics technology provides a powerful tool for understanding cellular heterogeneity and diversity within complex tissues.

Mitochondrial DNA (mtDNA) is small, without chromatin, and has multiple copies within a cell. Those features help us to get good mtDNA sequencing depth in single-cell ATAC seq data, which facilitates somatic mutation detections.

Together with chromatin openness data, mitochondrial somatic mutation can be used as a lineage label for helping us understand cellular differentiation, oncogenesis and development.

In this vignette, the scMitoMut package is used to identify and visualize mtDNA single nucleic somatic mutations.

In the following analysis, scMitoMut was used to analyze the allele count data, which is the output of mgatk. The only few loci have been selected for the demonstration purpose to reduce the file size and the run time.

The full demostration can be found XXX, which shows the analysis begins with mgatk and CellRanger output.

We load the allele count table with the `parse_table`

function. The allele count table consists with following columns:

`loc`

: the locus name`cell_barcode`

: the cell barcode of the cell`fwd_depth`

: the forward read count of the allele`rev_depth`

: the reverse read count of the allele`alt`

: the allele name`coverage`

: the total read count of the locus`ref`

: the reference allele name

Instead of using the table above as input, the output from the `mgatk`

also can be read in directly using the `parse_mgatk`

function.

Using the `parse_table`

function or `parse_mgatk`

function, the allele count data are transformed into an `H5`

file.
The `H5`

file works as a database, which does not occupy the memory, and data can be accessed by querying.
It helps for better memory usage and faster loading.

The process may take some minutes. The return value is the `H5`

file path.

```
## Load the allele count table
f <- system.file("extdata", "mini_dataset.tsv.gz", package = "scMitoMut")
f_h5_tmp <- tempfile(fileext = ".h5")
f_h5 <- parse_table(f, h5_file = f_h5_tmp)
```

`f_h5`

`## [1] "/tmp/RtmpeS06VE/file75c1919d2aa7e.h5"`

Once we have the `H5`

file, we can use `open_h5_file`

function to load it and get an object named “mtmutObj”.

**Detail**: On this step, the `mtmutObj`

has 6 slots
- `h5f`

is the `H5`

file handle
- `mut_table`

is the allele count table
- `loc_list`

is a list of available loci
- `loc_selected`

is the selected loci
- `cell_list`

is a list of available cell ids
- `cell_selected`

is the selected cell ids

```
## Open the h5 file as a scMitoMut object
x <- open_h5_file(f_h5)
str(x)
```

```
## List of 9
## $ file : chr "/tmp/RtmpeS06VE/file75c1919d2aa7e.h5"
## $ h5f :Formal class 'H5IdComponent' [package "rhdf5"] with 2 slots
## .. ..@ ID : chr "72057594037927938"
## .. ..@ native: logi FALSE
## $ mut_table :Formal class 'H5IdComponent' [package "rhdf5"] with 2 slots
## .. ..@ ID : chr "144115188075855873"
## .. ..@ native: logi FALSE
## $ loc_list : chr [1:16(1d)] "chrM.200" "chrM.204" "chrM.310" "chrM.824" ...
## $ loc_selected : chr [1:16(1d)] "chrM.200" "chrM.204" "chrM.310" "chrM.824" ...
## $ cell_list : chr [1:1359(1d)] "AAACGAAAGAACCCGA-1" "AAACGAAAGTACCTCA-1" "AAACGAACAGAAAGAG-1" "AAACGAAGTGGTCGAA-1" ...
## $ cell_selected: chr [1:1359(1d)] "AAACGAAAGAACCCGA-1" "AAACGAAAGTACCTCA-1" "AAACGAACAGAAAGAG-1" "AAACGAAGTGGTCGAA-1" ...
## $ loc_pass : NULL
## $ loc_filter :List of 5
## ..$ min_cell : num 1
## ..$ model : chr "bb"
## ..$ p_threshold : num 0.05
## ..$ alt_count_threshold: num 0
## ..$ p_adj_method : chr "fdr"
## - attr(*, "class")= chr "mtmutObj"
```

```
## Show what's in the h5 file
h5ls(x$h5f, recursive = FALSE)
```

```
## group name otype dclass dim
## 0 / cell_list H5I_DATASET STRING 1359
## 1 / cell_selected H5I_DATASET STRING 1359
## 2 / loc_list H5I_DATASET STRING 16
## 3 / loc_selected H5I_DATASET STRING 16
## 4 / mut_table H5I_GROUP
```

We are only interested in annotated good-quality cells.

So we will select the cells with annotation, which are good quality cells.

```
f <- system.file("extdata", "mini_dataset_cell_ann.csv", package = "scMitoMut")
cell_ann <- read.csv(f, row.names = 1)
## Subset the cells, the cell id can be found by colnames() for the Seurat object
x <- subset_cell(x, rownames(cell_ann))
```

After subsetting the cells, the `cell_selected`

slot will be updated.
Only the selected cells will be used in the following p-value calculation.

`head(x$cell_selected)`

```
## [1] "AAACGAAAGAACCCGA-1" "AAACGAAAGTACCTCA-1" "AAACGAACAGAAAGAG-1"
## [4] "AAACGAAGTGGTCGAA-1" "AAACGAATCAATCGTG-1" "AAACGAATCCCACGGA-1"
```

Similarly, we can select loci by using the `subset_locus`

function. It saves time when we only focus on a few loci.

Assuming an mtDNA locus’s majority base has constant allele frequency in cells, the allele frequency mutation was caused by sequencing error or sampling. Based on that assumption, we build a null hypothesis that there are no mutations for that locus in the cell. Then we fit the allele frequency distribution and calculate the probability of observing allele frequency for a specific locus in a cell. If the probability is small, we can reject the null hypothesis and conclude that there is a mutation for that locus in the cell.

To calculate the probability value (p-value), we run `run_calling`

function, which has 2 arguments:
- `mtmutObj`

is the `scMitoMut`

object
- `mc.cores`

is the number of CPU threads to be used

The process will take some time.
The output will be stored in the `pval`

group of the `H5`

file.
The result is stored in the hard drive, instead of in memory.
We don’t need to re-run the mutation calling when loading the `H5`

file next time.

The mutation calling is based on beta-binomial distribution. The mutation p-value is the probability that with the null hypothesis: there are no mutations for that locus in the cell.

**Detail**: For a specific locus, we calculate the p-value using the following steps.
1. Defining the wild-type allele as the allele with the highest average AF among cells.
2. Fitting a 2 components binomial-mixture model to the AF of the wild-type allele for all cells.
Using the model, we define the wild-type cells if it has a probability >= 0.001 to be the wild type.
3. Using those wild-type cells, we fit the beta-binomial model.
4. At last, based on the model, we calculate the p-value of observing the AF of the wild-type allele in the cell.

```
## Run the model fitting
run_model_fit(x, mc.cores = 1)
```

`## chrM.200`

`## chrM.204`

`## chrM.310`

`## chrM.824`

`## chrM.1000`

`## chrM.1001`

`## chrM.1227`

`## chrM.2285`

`## chrM.6081`

`## chrM.9429`

`## chrM.9728`

`## chrM.9804`

`## chrM.9840`

`## chrM.12889`

`## chrM.16093`

`## chrM.16147`

```
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 1126843 60.2 2115388 113 2115388 113.0
## Vcells 2168869 16.6 8388608 64 5681889 43.4
```

```
## The p-value is kept in the pval group of H5 file
h5ls(x$h5f, recursive = FALSE)
```

```
## group name otype dclass dim
## 0 / cell_list H5I_DATASET STRING 1359
## 1 / cell_selected H5I_DATASET STRING 1359
## 2 / loc_list H5I_DATASET STRING 16
## 3 / loc_selected H5I_DATASET STRING 16
## 4 / model_par_bb H5I_DATASET COMPOUND 16
## 5 / model_par_bm H5I_DATASET COMPOUND 16
## 6 / mut_table H5I_GROUP
## 7 / pval H5I_GROUP
```

Then we will filter the mutations by using the `mut_filter`

function with the following criteria:
- The mutation has at least 5 cells mutant.
- The FDR-adjusted p-value is less than 0.05.

The output is a `data.frame`

with 2 columns
- `loc`

is the locus
- `mut_cell_n`

is the cell number

We can see that there are 12 loci after filtering.

**Detail**: The `mut_filter`

function has 4 arguments:
- `mtmutObj`

is the `mtmutObj`

object
- `min_cell`

is the minimum number of mutant cells
- `p_adj_method`

is the method used to adjust the p-value.
- `p_threshold`

is the adjusted p-value threshold

```
## Filter mutation
x <- filter_loc(
mtmutObj = x,
min_cell = 2,
model = "bb",
p_threshold = 0.01,
p_adj_method = "fdr"
)
x$loc_pass
```

```
## [1] "chrM.200" "chrM.204" "chrM.310" "chrM.824" "chrM.1227"
## [6] "chrM.2285" "chrM.6081" "chrM.9429" "chrM.9728" "chrM.9804"
## [11] "chrM.9840" "chrM.12889" "chrM.16093" "chrM.16147"
```

We will visualize the mutation by heatmap using the `plot_heatmap`

function.
It can draw a heatmap of p-value, AF, or binarized mutation status.
Its input is the `mtmutObj`

object.
It will independently apply all the filters we used in the `mut_filter`

function, and select the cells and loci that pass the filter criteria.
In all kinds of figures, the mutation status will be calculated, and the loci and cells are ordered by the mutation status.

**Detail**: The `plot_heatmap`

arguments.
- `mtmutObj`

is the `scMitoMut`

object
- `pos_list`

is the list of loci
- `cell_ann`

is the cell annotation
- `ann_colors`

is the color of the cell annotation
- `type`

is the type of the heatmap which can be `p`

, `af`

, or `binary`

- `p_adj_method`

is the method used to adjust the p-value
- `p_threshold`

is the adjusted p-value threshold to determine if a cell has mutation when selecting the cells and loci
- `min_cell_n`

is the minimum number of cells that have mutation when selecting the cells and loci
- `p_binary`

is the adjusted p-value threshold to get the binary mutation status
- `percent_interp`

is the percentage overlap threshold between mutations, to determine if two mutations are correlated for interpolating the mutation status
- `n_interp`

is the number of overlapped cells to determine if two mutations are correlated for interpolating

The interpolation is based on the assumption that the mutation are unique, it is rare to have two mutation in the same population. Therefore, when two mutations are correlated, one of them is likely a subclone of the other one.

The binary heatmap will show the mutation status of each cell for each locus. The red means mutant, the blue means no mutant, and the grey means the missing value.

For the mutation status, if two mutations are correlated, we set the high-frequency locus will get the mutant status from the low-frequency one. The correlation is determined by the percentage and number of cells overlapping between the two mutations.

```
## Prepare the color for cell annotation
colors <- c(
"Cancer Epi" = "#f28482",
Blood = "#f6bd60"
)
ann_colors <- list("SeuratCellTypes" = colors)
```

```
## binary heatmap
plot_heatmap(x,
cell_ann = cell_ann, ann_colors = ann_colors, type = "binary",
percent_interp = 0.2, n_interp = 3
)
```

Also we can turn off the interpolation by setting `percent_interp = 1`

.

```
## binary heatmap
plot_heatmap(x,
cell_ann = cell_ann, ann_colors = ann_colors, type = "binary",
percent_interp = 1, n_interp = 3
)
```

The p-value heatmap shows the adjusted p-value of each cell for each locus. The order of the cells and loci is determined by the mutation status.

```
## p value heatmap
plot_heatmap(x,
cell_ann = cell_ann, ann_colors = ann_colors, type = "p",
percent_interp = 0.2, n_interp = 3
)
```