Single Sample DNA + RNA#

This notebook is set up to work with single-sample h5 files, generated with v3 chemistry. This notebook will only work with Mosaic versions 3.17 and higher.

Objective: To showcase the minimum number of steps required for tertiary analysis of DNA (single-cell genotyping and CNV) and RNA analytes and explore different ways of visualizing the data.

Major questions answered:

Can we identify DNA clones based on genotypes (SNVs/Indels)?
Do we detect CNV events (e.g., copy number amplification, copy number loss)?
What RNA expression patterns can we see?
Do DNA clones (genotypes) correlate with specific RNA-defined clusters and/or CNV profiles?

Sections:

Setup
Data Structure
DNA Analysis
CNV Analysis
RNA Analysis
Combined Visualizations
Export and Save Data
Appendix

Not shown: All available methods and options - documented here

Setup#

Topics covered

Load dependencies
Load data

Load Dependencies#

NOTE: importing dependencies can sometimes take a couple of minutes.

import warnings
warnings.filterwarnings('ignore')

# Import mosaic libraries
import missionbio.mosaic as ms
from missionbio.mosaic.assay import _Assay
from missionbio.mosaic.constants import READS, LABEL, NORMALIZED_READS

# Import graph_objects from the plotly package to display figures when saving the notebook as an HTML
import plotly as px
import plotly.graph_objects as go

# Import additional packages for specific visuals
import matplotlib.pyplot as plt
import plotly.offline as pyo
pyo.init_notebook_mode()
import pandas as pd
import numpy as np
import seaborn as sns
from missionbio.plotting.multimap import MultiMap
from missionbio.plotting.heatmap import Heatmap
import plotly.graph_objects as go
from IPython.display import display, HTML
from sklearn.decomposition import PCA

# Import dependencies to create a new report
from pathlib import Path
from missionbio.config.config import Config
from missionbio.tertiary.config.prepare import DEFAULTS
from missionbio.tertiary.step.report import create_single_patient_report

from typing import Sequence, Optional, Tuple

pd.set_option("display.max_rows", 400)

# Note: when exporting the notebook as an HTML, plots that use the "go.Figure(fig)" command are saved

NOTE: The next two cells define variables and two functions to be used later on in the RNA analysis clustering section of the notebook.

# Define random state, umap metric and heatmap colorscales
RANDOM_STATE = 3
UMAP_METRIC = 'cosine'
HEATMAP_COLORSCALE = "Viridis"

# Cell is collapsed to hide length of code
# Build functions for elbow_max_distance and get number of PCs
def elbow_max_distance(x: Sequence, y: Sequence) -> int:
    """
    Find the elbow point of a decreasing curve using the maximum distance to chord method.
    
    Parameters
    ----------
    x : Sequence
    y : Sequence

    Returns
    -------
    idx : int
        Index of the elbow point in the original arrays.
    """
    if len(x) != len(y):
        raise ValueError("x and y must be sequences of equal length")

    # sort by decreasing y
    order = np.argsort(y)[::-1]
    x, y = x[order], y[order]

    # line through endpoints
    x0, y0 = x[0], y[0]
    x1, y1 = x[-1], y[-1]

    # calc perpendicular distances
    numerator = np.abs((y1 - y0) * (x - x0) - (x1 - x0) * (y - y0))
    denominator = np.hypot(y1 - y0, x1 - x0)
    d = numerator / denominator

    k = int(np.argmax(d))
    return order[k]


def get_n_PCs(
    assay: _Assay,
    layer: str = NORMALIZED_READS,
    plot_height: int = 700,
    plot_width: int = 700,
    random_state: Optional[int] = RANDOM_STATE,
) -> Tuple[int, go.Figure]:
    """
    Get the number of PCs for clustering using the elbow method. Return the number of PCs selected
    and scatterplot.

    Parameters
    ----------
    assay : _Assay
        missionbio.mosaic.assay._Assay with data for PCA analysis
    layer: str
        The _Assay layer to use for PCA
    plot_height : int = 700
        Height of plotly scatterplot
    plot_width : int = 700
        Width of plotly scatterplot
    random_state: int or None
        Optional random state for PCA

    Returns
    -------
    selected_PCs : int
        The selected number of PCs
    fig : go.Figure
        A scatterplot showing variance per PC ordered from greatest to least along the x, with the
        PC at the elbow highlighted
    """
    data = assay.get_attribute(NORMALIZED_READS, constraint="row").values
    components = min(data.shape)

    pca = PCA(n_components=components, random_state=RANDOM_STATE)
    pca.fit(data.T)
    pca_data = pca.components_.T

    elbow = elbow_max_distance(np.arange(len(pca.explained_variance_)), pca.explained_variance_)
    selected_PCs = elbow + 1 # elbow is an index, add one for number of selected PCs
    components_display = np.minimum(components, 50)

    scatter = go.Scatter(
        x = np.arange(components_display), 
        y = pca.explained_variance_ratio_[:components_display],
        mode="markers",
    )
    fig = go.Figure(data=[scatter])
    fig.add_vline(x=elbow, line_width=1.5, line_dash="dash", line_color="red", layer="below")
    fig.update_layout(
        xaxis=dict(title="Principal component", zeroline=False),
        yaxis=dict(title="Explained variance ratio", zeroline=False),
        height=plot_height,
        width=plot_width,
    )

    return selected_PCs, fig

# This code is optional, but will make the notebook cells/figures display across the entire width of the browser
from IPython.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

# Check version; this notebook is designed for Mosaic v3.17 or higher
print(ms.__version__)

# Any function's parameters and default values can be looked up via the 'help' function
# Here, the function is 'ms.load'
help(ms.load)

Load the data#

In the ms.load statement, there are several arguments:

raw : always set raw=False; if raw=True, ALL barcodes will be loaded (rather than cell-associated barcodes)
filter_variants : always set filter_variants=True unless an expected (target) variant cannot be found. Additional filtering options are included in the DNA section below. NOTE : if loading using filter_variants=False, the whitelist may need to set to None or removed from the command, otherwise it may not load variants in
single : always set single=False for multisample/merged h5 files
whitelist : The whitelist option loads any variant that is in the vcf.gz file (e.g. “chr1:179520511:C/G”). The whitelist argument takes many variant formats, they are detailed here
filter_cells : Always set filter_cells=False; if filter_cells=True, only cells called by the completeness algorithm are loaded

# Specify the h5 file to be used in this analysis: h5path = '/path/to/h5/file/test.h5'
# If working with Windows, an 'r' may need to be added before the path: h5path = r'/path/to/h5/file/test.h5'
h5path = '/Users/botero/Desktop/RNA notebooks/INR121.h5'

# Load the data
sample = ms.load(h5path, raw=False, filter_variants=True, single=True, whitelist = [], filter_cells=False) 

Data Structure#

DNA, CNV, and RNA are sub-classes of the Assay class. The information is stored in four categories, and the user can modify each of those:

1. metadata (add_metadata / del_metadata):

dictionary containing assay parameters and basic performance metrics

2. layers (add_layer / del_layer):

dictionary containing matrices of assay metrics
all the matrices have the shape (num barcodes) x (num ids)
for DNA assays, this includes AF, GQ, DP, etc. (per cell, per variant)
for CNV assays, this includes read counts (per cell, per amplicon)
for RNA assays, this includes read counts (per cell, per gene)

3. row_attrs (add_row_attr / del_row_attr):

dictionary containing ‘barcode’ as one of the keys (where the value is a list of all barcodes)
for all other keys, the values must be of the same length, i.e. match the number of barcodes
this is the attribute where ‘label’, ‘pca’, and ‘umap’ values are added

4. col_attrs (add_col_attr / del_col_attr):

dictionary containing ‘id’ as one of the keys (where the value is a list of all ids)
for DNA assays, ‘ids’ are variants; for CNV assays, ‘ids’ are genes; for RNA assays, ‘ids’ are amplicons
for all other keys, the values must be of the same length, i.e. match the number ids

# Summary of DNA assay 
print("\'sample.dna\':", sample.dna, '\n')
print("\'row_attrs\':", "\n\t", list(sample.dna.row_attrs.keys()), '\n')
print("\'col_attrs\':", "\n\t", list(sample.dna.col_attrs.keys()), '\n')
print("\'layers\':", "\n\t", list(sample.dna.layers.keys()), '\n')
print("\'metadata\':", "\n")
for i in list(sample.dna.metadata.keys()):
    print("\t", i, ": ", sample.dna.metadata[i], sep="")

# For DNA, ids are variants
# sample.dna.ids() is a shortcut for sample.dna.col_attrs['id']
sample.dna.ids()

# Summary of CNV assay
print("\'sample.cnv\':", sample.cnv, '\n')
print("\'row_attrs\':", "\n\t", list(sample.cnv.row_attrs.keys()), '\n')
print("\'col_attrs\':", "\n\t", list(sample.cnv.col_attrs.keys()), '\n')
print("\'layers\':", "\n\t", list(sample.cnv.layers.keys()), '\n')
print("\'metadata\':", "\n")
for i in list(sample.cnv.metadata.keys()):
    print("\t", i, ": ", sample.cnv.metadata[i], sep="")

# For CNV, ids are amplicons
# sample.cnv.ids() is a shortcut for sample.cnv.col_attrs['id']
sample.cnv.ids()

# Summary of RNA assay 
print("\'sample.rna\':", sample.rna, '\n')
print("\'row_attrs\':", "\n\t", list(sample.rna.row_attrs.keys()), '\n')
print("\'col_attrs\':", "\n\t", list(sample.rna.col_attrs.keys()), '\n')
print("\'layers\':", "\n\t", list(sample.rna.layers.keys()), '\n')
print("\'metadata\':", "\n")
for i in list(sample.rna.metadata.keys()):
    print("\t", i, ": ", sample.rna.metadata[i], sep="")

# For RNA, ids are genes
# sample.rna.ids() is a shortcut for sample.rna.col_attrs['id']
sample.rna.ids()

DNA Analysis#

Topics covered

Filtering of DNA variants
Subsetting dataset for variants of interest, including whitelisted variants
Addition of annotations to the variants
Manual variant selection and clone identification
Visualizations and customization options

Basic filtering#

For filtering DNA variants, there are two methods:
Method 1: Utilize the function sample.dna.filter_somatic_variants. This function automatically filters variants, eliminating germline and non-pathogenic variants from the data. More about the function can be seen here. This method is used by many of our tertiary analysis pipelines, such as Clonal Insights.
Method 2: Manually filter the variants using sample.dna.filter_variants function. This method allows for more flexibility with filters, but retains all germline and non-pathogenic variants in the data.

There are many options for filtering DNA variants. Use the help() function to understand the approach listed below.

Method 1#

# Filter germline and nonpathogenic variants
dna_vars = sample.dna.filter_somatic_variants()
sample.dna = sample.dna[:, dna_vars]
dna_vars = sample.dna.filter_variants()

Method 2#

# Manually adjust filters if needed by overwriting dna_vars
dna_vars = sample.dna.filter_variants(
    min_dp=10,
    min_gq=30,
    vaf_ref=5,
    vaf_hom=95,
    vaf_het=20,
    min_prct_cells=50,
    min_mut_prct_cells=1
)

# Check the number of filtered variants. When using the default filters, the number of 
# variants is likely smaller compared to the originally loaded variants due to the more 
# stringent filtering criteria (e.g., vaf_ref=5, vaf_hom=95, vaf_het=20).
print('Number of variants:', len(dna_vars))
print(dna_vars)

NOTE: Variants that are whitelisted during sample loading may be removed at this step, but can be added back in below. Whitelisted variants can be added to the list of final variants, or used exclusively in the code below.

# Ensure correct nomenclature, ie whitelist = ["chr13:28589657:T/G","chrX:39921424:G/A"]
# While the load statement accepts many variant formats, this whitelist does not
# If there are no whitelist variants, this can be left blank

white_list = []

# Combine whitelisted and filtered variants
final_vars = list(set(list(dna_vars) + white_list))

# To use only white_list variants:
# final_vars = white_list

#Check the length of the final list of variants
len(final_vars)

# Dimensionality of the original sample.dna dataframe
# First number = number of cells (rows); second number = number of variants (columns)
sample.dna.shape

# Before subsetting, verify that all the chosen variants are in the current sample.dna ids (should return True)
print(set(final_vars).issubset(set(sample.dna.ids())))

NOTE: If this returns False, the data will not subset properly. Please resolve before moving forward.

# Subsetting sample.dna (columns) based on reduced variant list.
sample.dna = sample.dna[:, final_vars]

# Check the shape of the final filtered DNA object, i.e. (number of barcodes/cells, number of ids/variants)
sample.dna.shape

Annotation Addition#

help(sample.dna.get_annotations)

# Fetch annotations using varsome
# Note: run this on a filtered DNA sample - too many variants (e.g., 100+) are not handled correctly by the method
ann = sample.dna.get_annotations() 

Variant selection and Subclone identification#

In this section of the notebook, all variants remaining in the dataset will populate in a variant table. This table is interactive, variants can be selected, and rows can be sorted by ascending/descending values. The variant name can be clicked on and will navigate to the variant’s varsome page in the default browser.

Variants selected in this table will populate in a subclone table below.
The variants in the subclone table can be highlighted and assessed for Read Depth, Genotype Quality and Variant Allele Frequency.
Subclones can be renamed by clicking on the pencil icon.
ADO score: The ADO score threshold can be adjusted, but by default is set to 1. Any clones with an ADO score higher this will be moved into the “ADO subclones” column. We recommend moving any clones with a score of .8 or higher into this column. The algorithm used to generate the score is detailed here.
Min clone size: The Min Clone Size can also be adjusted, but by default is set to 1. Any clones that represent less than 1% of the sample will be moved into the “Small Subclones” column.
Cells with missing genotype information across any of the selected variants will be moved into the “Missing GT” column.

More on this workflow can be found here.

# Run the variant table workflow to select variants and begin clone identification
wfv = ms.workflows.VariantSubcloneTable(sample)
wfv.run()
# The width can be adjusted if needed for long variant names
# wfv.run(width=3000)

NOTE: When continuing onto the next piece of the notebook, the workflow will stop working.

# Subsample the variants to only the variants selected from the workflow
variants = wfv.selected_variants
sample.dna = sample.dna[:,variants]

# Add annotation to the id names
sample.dna.set_annotated_ids(add_clinvar=False)

# Annotations are now added to the variants
sample.dna.ids()

# Use sample.dna.reset_ids() to get the original ids

# Plot heatmap using NGT_FILTERED
fig = sample.dna.heatmap(attribute='NGT_FILTERED')
go.Figure(fig)

# Clone removal
# Remove barcodes from the missing GT, small subclones, ADO subclones or FP labels
clones = ['missing', 'small', 'ADO', 'FP']
for c in clones:
    cells = sample.dna.barcodes({c})
    if len(cells) > 0:
        sample.dna = sample.dna.drop(cells)
set(sample.dna.get_labels())

# Redraw heatmap
fig = sample.dna.heatmap(attribute='NGT_FILTERED')
go.Figure(fig)

# Evaluate new total number of cells after the above filtering
sample.dna.shape

Adjusting subclone colors or heatmap colors#

To change the color palette used for the subclones, use set_palette. For this provide a list of ALL subclone names, pointing to the hexcode/color to change that subclone to. Note: this function also works with the protein and CNV assays. See below for an example of this.

To change the color palette used for any of the assays/layers, use ms.Config.Colorscale. Then list the assay and layer, and assign the plotly colorscale to use for those graphics. Visualize all available colorscales by using plotly.colors.sequential.swatches_continuous().show(). To reset all color palettes back to their defaults, use ms.Config.Colorscle.reset().

# Example of changing the colors assigned to each clone
sample.dna.set_palette({'Clone 1': '#800080', 'Clone 2': '#FF69B4'})

# Change colorScale to change the heatmap colors
# This will change the color of the dna.ngt layer to be viridis
ms.Config.Colorscale.Dna.NGT = 'viridis'
fig = sample.dna.heatmap('NGT')
fig.show()

# Run the following to reset the colorscale
#ms.Config.Colorscale.Dna.reset("NGT")
# Or run this to reset any/all modifications:
#ms.Config.Colorscale.reset()
# Colorscale can be changed for DNA, CNV and Protein assays

To return back to the original population of cells, reset the data using: sample.reset("dna") This command ‘sample.reset’ works on all assays, including CNV and protein.

CNV Analysis#

Topics covered

Amplicon and cell filtering and ploidy estimation
Visualization of ploidy across subclones present

# Get gene names for amplicons
sample.cnv.get_gene_names();

CNV workflow#

This workflow will normalize all reads and filter amplicons/cells based on the settings set at the beginning of the workflow:

Amplicon completeness: refers to the minimum percentage of cells that must have reads greater than or equal to the minimum read depth set. By default this is set to 50.
Amplicon read depth: refers to the minimum read depth for each amplicon-barcode combination to not be considered missing. By default this is set to 10.
Mean cell read depth: refers the minimum mean read depth for a cell to be included in the analysis, otherwise the cell will be removed. By default this is set to 40.
Diploid clone in DNA: refers to which subclone is used as the true diploid population. All ploidy estimates will be calculated in relation to this diploid population. We recommend setting this to the ‘WT’ population or earliest clone present.

NOTE: If large copy number events are expected, Amplicon Completeness and Amplicon Read Depth may need to be reduced.

Once the above filters are set, the visualizations can be changed.

Plot: Can be changed from Heatmap positions, to Heatmap genes, Line-plot positions, Line-plot genes.
Clone for line plot: If one of the line-plot visualizations is selected, only one clone can be shown at a time. This determines which one is plotted.
X-axis features: To plot a subset of the data (chromosomes or genes), select which chromosomes/genes should be plotted with this function. Chromosomes can be selected for ‘positions’ type plots, and Genes can be selected for ‘genes’ type plots.

More on this workflow can be found here.

# CNV workflow to filter, normalize and estimate ploidy
wfc = ms.workflows.CopyNumber(sample)
wfc.run()

NOTE: When continuing onto the next piece of the notebook, the workflow will stop working.

# Heatmap with the features ordered by the default amplicon order
# To plot a subset of chromosome, the chromosomes can be put in list format in the features argument
fig = sample.cnv.heatmap('ploidy', features='positions') #features=['7', '17', '20']

# Optionally, restrict the range of ploidy values based on observed/expected CNV events (commented out)
#fig.layout.coloraxis.cmax = 4
#fig.layout.coloraxis.cmin = 0

# Optionally, change the size of the figure:
#fig.layout.width = 1600
#fig.layout.height = 1500

go.Figure(fig)

# Heatmap with the features grouped by the genes
# To plot just a subset of genes, put them in list format for features
fig = sample.cnv.heatmap('ploidy', features='genes', convolve=1) #features=["ASXL1", "EZH2",'TP53']

# Optionally, update the separating lines to be black
#for shape in fig.layout.shapes:
#    shape.line.color = '#000000'

go.Figure(fig)

# Show heatmap with convolve and subclustering turned off
bars = sample.cnv.clustered_barcodes('ploidy', subcluster=False)

# This is useful to create "convolved" heatmaps which are easier to interpret
# With the subclustering off and convolve=20, the noise will be reduced and real signals will be easier to visualize
fig = sample.cnv.heatmap('ploidy', bars_order=bars, convolve=20)
fig.layout.width = 900
fig.show()

RNA Analysis#

Topics covered

Data normalization
Preparing for data clustering
Clustering
Differential expression

Normalization#

# To get more information on the available parameters in rna.normalize_reads
help(sample.rna.normalize_reads)

# Normalize RNA reads
sample.rna.normalize_reads(
    correct_background = False, 
    negative_control_genes = ["BFP", "EGFP", "RFP"], 
    bg_method = "max", 
    min_reads_per_cell = 5,  
    min_cells_detec_gene = 0.05, 
    max_fraction = 0.05, 
    exclude_highly_expressed = False
)

Preparing for clustering#

Select PCs for clustering and UMAP#

The number of PCs used has a large influence on the RNA clusters that are called. Using fewer PCs selects the genes with the largest effect on variation, but may miss some important sources of variance. Using more PCs incorporates more genes into the clustering, but may result in splitting biologically relevant groups.

The below plot shows the number of PCs selected by the knee algorithm. This is how the tertiary module selects PCs, however, this can be changed for custom analysis.

# Determine the number of PCs to use for PCA
n_PCs, fig = get_n_PCs(sample.rna)

print(f"{n_PCs} PCs selected by knee method")
fig.show()

# Run the PCA and then the UMAP
# This step may take some time to run
sample.rna.run_pca(
    attribute=NORMALIZED_READS,
    components=n_PCs, # manually tune number of PCs if desired
    random_state=RANDOM_STATE,
)

sample.rna.run_umap(
    attribute='pca', 
    n_neighbors=50, # UMAP n neighbors may also be tuned to adjust plot
    metric="cosine",
    min_dist=0, 
    random_state=RANDOM_STATE,
)

Clustering#

The tertiary analysis module uses Leiden clustering since it was found to produce the most consistently well-connected and biologically meaningful clusters. However, other clustering methods can be used.

In the clustering, UMAP can be used if desired in place of PCA. Additionally, the n_neighbors parameter can be used to modify the clustering.

sample.rna.cluster(
    'pca', 
    method="leiden", 
    n_neighbors=30, 
    random_state=RANDOM_STATE
)

fig = sample.rna.scatterplot('umap', colorby='label')
fig.show(renderer='png')

# RNA UMAP colored by one gene
gene = "B2M"
fig = sample.rna.scatterplot('umap', colorby=NORMALIZED_READS, features=[gene])
fig.update_layout(title=gene)

# Heatmap of RNA normalized reads
fig = sample.rna.heatmap(attribute=NORMALIZED_READS)
fig.update_layout(width=1200, coloraxis=dict(colorscale=HEATMAP_COLORSCALE))
fig.show(renderer='png')

# Generate a filtered RNA heatmap
# Plots genes detected in at least 5 cells with a normalized reads greater than 0
norm_reads = sample.rna.get_attribute(NORMALIZED_READS)
cells_detected = (norm_reads > 0).sum(axis=0)
min_cells_detected = 5
detected_genes = cells_detected[cells_detected >= min_cells_detected].index.tolist()
fig = sample.rna[:, detected_genes].heatmap(attribute=NORMALIZED_READS)
fig.update_layout(width=1200, coloraxis=dict(colorscale=HEATMAP_COLORSCALE))

Differential expression#

Compare the expression of genes between clusters.

# Get the list of the RNA markers
# This is a list of marker genes for each RNA cluster using DE analysis
cm = sample.rna.get_markers() 
# The number of markers can be changed
# cm = sample.rna.get_markers(top_n_markers=10) 

# Get the marker genes
np.unique(cm.top_markers['group'])
marker_genes = list(pd.unique(cm.top_markers['gene']))
marker_genes

# Heatmap of RNA normalized reads, only showing marker genes
fig = sample.rna[:, marker_genes].heatmap(attribute=NORMALIZED_READS)
fig.update_layout(width=1200, coloraxis=dict(colorscale=HEATMAP_COLORSCALE))
fig.show(renderer='png')

# UMAP colored by the expression of each marker gene
fig = sample.rna.scatterplot(attribute='umap',
                           colorby='normalized_counts',
                           features=marker_genes)
go.Figure(fig)

# Dotplot by RNA cluster
cm.dotplot()

# Name the RNA clusters
# Combine clusters by assigning identical names
# Just an example! Need to fill this in based on knowledge of the sample and markers
sample.rna.rename_labels(
    {
        '1': 'Cell type A', 
        '2': 'Cell type B', 
        '3': 'Cell type C',
        '4': 'Cell type D',
        '5': 'Cell type E'
    }
)

# RNA violinplot for one gene by RNA cluster
gene = "B2M"
fig = sample.rna.violinplot(NORMALIZED_READS, splitby="label", features=[gene])
fig.show()

# Marker gene violin plot of one cluster versus all other clusters
fig = cm.marker_violin_plot('1')
fig.update_layout(width=900)
fig.show()

Combined Visualizations#

In this section the dateset will be subset to only retain barcodes with remaining DNA, CNV and RNA data. Then this data can be plotted together using sample.heatmap(), sample.signaturemap() and other functions.

Overlaps in the dataset across different assays can be further quantified and visualized, using sample.crosstabmap().

# This will return the barcodes common to all assays in the sample.
sample.common_barcodes()

# Use that to filter the sample so that only the common barcodes are present in all assays
sample = sample[sample.common_barcodes()]

# Check dimensionality for each assay; the number of cells (first number) should be the same in each data set
print(sample.dna.shape,
sample.cnv.shape,
sample.rna.shape)

# Add DNA clone information to RNA assay
sample.rna.add_row_attr('clone', sample.dna.get_labels())
cm_dna = sample.rna.get_markers(split_by="clone")

# DNA, CNV, RNA combined heatmap
fig = sample.heatmap(
    clusterby=('dna', 'cnv', 'rna'),  # The first assay is used for the labels
    attributes=['AF', 'ploidy', 'normalized_counts'],
    features=[None, 'genes', None],  # If None, then clustered_ids is used
    bars_order=None,  # The order of the barcodes
    widths=None,
    order=('dna', 'cnv', 'rna')  # The order in which the heatmaps should be drawn
)
fig.show()

# DNA and RNA combined heatmap
fig = sample.heatmap(
    clusterby=('dna', 'rna'),  # The first assay is used for the labels
    attributes=['AF', 'normalized_counts'],
    features=[None, None],  # If None, then clustered_ids is used
    bars_order=None,  # The order of the barcodes
    widths=None,
    order=('dna', 'rna')  # The order in which the heatmaps should be drawn
)
fig.show()

# DNA and RNA signaturemap
# Plot a combined signature heatmap, showing DNA, CNV and RNA signatures for all subclones
# To only plot 2 of the analytes 
# Adjust the code below by removing the arguments for the analyte that should be removed 
sample.signaturemap(
    clusterby=('dna', 'rna'),
    attributes=('NGT', 'normalized_counts'),
    features=[None, None],
    signature_kind=['median', 'median'],
    widths=None,
    order=('dna', 'rna')  # The order in which the heatmaps should be drawn
)

# RNA UMAP colored by DNA clone
fig = sample.rna.scatterplot('umap', colorby='clone')
fig.show(renderer='png')

# RNA heatmap by DNA subclone
fig = sample.rna.heatmap(attribute='normalized_counts',splitby=sample.dna.row_attrs['label'])
go.Figure(fig)

# RNA dotplot by DNA clone
cm_dna.dotplot()

# Plot the RNA UMAP and color with ploidy for the defined gene
gene_ploidy = sample.cnv.signature("ploidy", splitby=None, features="gene_name")
tp53_ploidy = gene_ploidy["TP53"].values
sample.rna.scatterplot("umap", colorby=tp53_ploidy)

# Plot the RNA ridge plots split by DNA clone
sample.rna.ridgeplot(attribute='normalized_counts',
                         features=['CD3G', 'CD3D', 'ZAP70', 'ATM'],
                         splitby=sample.dna.row_attrs['label'])

# Plot the RNA ridge plots split by RNA cluster
sample.rna.ridgeplot(attribute='normalized_counts',
                         features=['CD3G', 'CD3D', 'ZAP70', 'ATM'],
                         splitby=sample.rna.row_attrs['label'])

Quantify overlaps between DNA clones and RNA clusters#

# Use crosstab() to calculate the overlap between DNA subclones and RNA Clusters
# This will generate the dataframe
# What percentage of each RNA cluster is defined by each Subclone
sample.dna.crosstab(sample.rna.get_labels(), normalize='columns')

# This will generate the heatmap of the dataframe generated above
# What percentage of each RNA cluster is defined by each Subclone
sample.dna.crosstabmap(sample.rna.get_labels(), normalize='columns').show()

Export and Save Data#

In this section the data can export a filtered .h5 file, which will contain all new labels/layers, and contain only the filtered barcodes/cells remaining. All data can be exported (row attributes, column attributes and layers) for every assay (DNA, CNV and RNA) into easily parsable .csv tables. Additionally, an updated .html report can be generated based on the analysis performed in this notebook.

# Save new h5 file that includes only the final, cleaned dataset
ms.save(sample, 'FilteredData.h5')

# Export data into csv formats
ms.to_zip(sample, 'filename')

# Generate an updated .html report
config = Config.read({}, DEFAULTS)

create_single_patient_report(
    sample,
    location=Path("/Users/botero/Desktop/test_report.html"),
    config=config,
)

Appendix#

SPARC imputation#

This section of the notebook is OPTIONAL

SPARC (Single-cell Phylogeny And Reconstruction of Clones) is a Mission Bio developed tool to infer the evolutionary history and clonal structure of a population, integrating information from SNVs, Loss of Heterozygosity (LoH), and Copy Number Variations (CNVs).

More about this tool can be found here.

# Import SPARC from missionbio.demultiplex.phylogeny.sparc import SPARC # Reset variant ids in the DNA assay sample.dna.reset_ids()

# Define the set of variants you want SPARC to run off of # Varaints must be in the format: chromosome:position:ref/alt, Ex: "chr4:55599436:T/C" # If you want to use all variants identified from the VariantSubcloneWorkflow you can use the list "variants", Ex: vars = variants vars = variants #vars = [] sparc = SPARC() result = sparc.predict(sample.dna, variants=vars)

# Order the results and produce a phylogeny tree ordered_result = result.serialize_clones() sample.dna.set_labels(ordered_result.labels.values) fig = ordered_result.draw() # Draws the phylogenetic tree fig.show()

# Plot a new heatmap with new DNA subclone labels fig = sample.dna.heatmap(attribute='NGT_FILTERED') fig.show()

Single Sample DNA + RNA

Contents

Single Sample DNA + RNA#

Setup#

Load Dependencies#

Load the data#

Data Structure#

DNA Analysis#

Basic filtering#

Method 1#

Method 2#

Annotation Addition#

Variant selection and Subclone identification#

Adjusting subclone colors or heatmap colors#

CNV Analysis#

CNV workflow#

RNA Analysis#

Normalization#

Preparing for clustering#

Select PCs for clustering and UMAP#

Clustering#

Differential expression#

Combined Visualizations#

Quantify overlaps between DNA clones and RNA clusters#

Export and Save Data#

Appendix#

SPARC imputation#