Analysis

Swan has several analysis options to use.

import swan_vis as swan

sg = swan.read('data/swan.p')
Read in graph from data/swan.p

Differential expression tests

Swan's old differential gene and transcript expression tests using diffxpy have now been deprecated as it seems that the library is unsupported. I recommend that users interested in running differential gene or transcript expression test either use PyDESeq2 or Scanpy's rank_genes_groups test, which both support the AnnData format for simple compatibility with Swan's AnnData expression representation.

Using PyDESeq2

Please read the PyDESeq2 documentation for details on how to use one of the SwanGraph AnnData objects to obtain differential expression results. Below is an example on how to find differentially expressed transcripts between cell lines.

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats
import numpy as np

adata = sg.adata.copy()

# PyDESeq2 currently doesn't support column names with underscores, so change that
adata.obs.rename({'cell_line': 'cellline'}, axis=1, inplace=True)
obs_col = 'cellline'

threads = 8

# densify matrix
adata.X = np.array(adata.X.todense())

# run test
dds = DeseqDataSet(adata=adata,
               design_factors=obs_col,
               n_cpus=threads,
               refit_cooks=True)
dds.deseq2()
stat_res = DeseqStats(dds,
                  n_cpus=threads)
stat_res.summary()

df = stat_res.results_df
Fitting size factors...
... done in 0.00 seconds.

Fitting dispersions...
... done in 34.09 seconds.

Fitting dispersion trend curve...
... done in 12.01 seconds.

Fitting MAP dispersions...
... done in 13.97 seconds.

Fitting LFCs...
... done in 4.86 seconds.

Refitting 0 outliers.

Running Wald tests...
... done in 2.36 seconds.

Log2 fold change & Wald test p-value: cellline hffc6 vs hepg2

34814 rows × 6 columns

df.head()

Isoform switching / Differential isoform expression testing

Isoform switching / differential isoform expression (DIE) testing is implemented according to the strategy in Joglekar et. al., 2021. DIE can roughly be described as finding statistically significant changes in isoform expression between two conditions along with a change in percent isoform usage per gene.

Pairwise comparisons can be set up using different columns in the metadata that was added to the SwanGraph with the obs_col and obs_conditions arguments.

# look at valid metadata options
sg.adata.obs
# find genes that exhibit DIE between HFFc6 and HepG2
obs_col = 'cell_line'
obs_conditions = ['hepg2', 'hffc6']
die_table, die_results = sg.die_gene_test(obs_col=obs_col,
                                          obs_conditions=obs_conditions,
                                          verbose=True)
Testing for DIE for each gene: 100%|██████████| 14684/14684 [03:50<00:00, 123.69it/s]

The resultant table contains an entry for each gene with the p value (p_val), adjusted p value (adj_p_val), and change in percent isoform usage for the top two isoforms (dpi), as well as the identities of the top isoforms involved in the switch. Exact details on these calculations can be found in Joglekar et. al., 2021.

die_table.head()

Differential isoform expression testing results are stored automatically in sg.adata.uns, and will be saved to the SwanGraph if you save it. You can regenerate this key and access the summary table by running the following code:

# die_iso - isoform level differential isoform expression test results
uns_key = swan.make_uns_key('die',
                            obs_col=obs_col,
                            obs_conditions=obs_conditions)
test = sg.adata.uns[uns_key]
test.head(2)

Swan comes with an easy way to filter your DIE test results based on adjusted p value and dpi thresholds.

test = sg.get_die_genes(obs_col=obs_col, obs_conditions=obs_conditions,
                       p=0.05, dpi=10)
test.head()

Swan also now automatically tracks transcription start site (TSS) and transcription end site (TES) usage, and can find genes that exhibit DIE on the basis of their starts or ends. To do this, use the kind argument to die_gene_test.

# find genes that exhibit DIE for TSSs between HFFc6 and HepG2
die_table, die_results = sg.die_gene_test(kind='tss',
                                          obs_col=obs_col,
                                          obs_conditions=obs_conditions,
                                          verbose=True)
die_table.head()
Testing for DIE for each gene: 100%|█████████▉| 14674/14684 [02:18<00:00, 162.60it/s]

To access the die results on the tss level, use die_kind='tss' as input to make_uns_key().

# die_iso - TSS level differential isoform expression test results
uns_key = swan.make_uns_key('die',
                            obs_col=obs_col,
                            obs_conditions=obs_conditions,
                            die_kind='tss')
test = sg.adata.uns[uns_key]
test.head(2)

And provide the kind='tss' option to get_die_genes() when trying to filter your test results.

test = sg.get_die_genes(kind='tss', obs_col=obs_col,
                        obs_conditions=obs_conditions,
                        p=0.05, dpi=10)
test.head()

For TESs, use kind='tes' as input to die_genes_test(), die_kind='tes' to make_uns_key(), and kind='tes' to get_die_genes().

# find genes that exhibit DIE for TESs between HFFc6 and HepG2
die_table, die_results = sg.die_gene_test(kind='tes',
                                          obs_col='cell_line',
                                          obs_conditions=['hepg2', 'hffc6'])
die_table.head()
Testing for DIE for each gene: 100%|██████████| 14684/14684 [02:19<00:00, 105.51it/s]
# die_iso - TES level differential isoform expression test results
uns_key = swan.make_uns_key('die',
                            obs_col=obs_col,
                            obs_conditions=obs_conditions,
                            die_kind='tes')
test = sg.adata.uns[uns_key]
test.head(2)
test = sg.get_die_genes(kind='tes', obs_col=obs_col,
                        obs_conditions=obs_conditions,
                        p=0.05, dpi=10)
test.head()

Finally, for transcriptomes generated with Cerberus, Swan automatically tracks intron chain usage, and you can perform intron chain switching analysis with kind='ic' as input to the die_gene_test() function.

sg_brain = swan.read('swan_modelad.p')

# find genes that exhibit DIE for ICs between genotypes
die_table, die_results = sg_brain.die_gene_test(kind='ic',
                                                obs_col='genotype',
                                                obs_conditions=['b6n', '5xfad'])
die_table.head()
Read in graph from swan_modelad.p

Combining metadata columns

What if none of the metadata columns you have summarize the comparisons you want to make? What if I want to find differentially expressed genes or transcripts, or find isoform-switching genes between hffc6 replicate 3 and hepg2 replicate 1?

sg.adata.obs.head()

Let's ignore for a moment the fact that the dataset column does effectively capture both replicate as well as cell line metadata. This may not always be the case with more complex datasets. Swan has a function to concatenate columns together and add them as an additional column to the metadata tables. Use the following code to generate a new column that concatenates as many preexisting metadata columns as you wish:

col_name = sg.add_multi_groupby(['cell_line', 'replicate'])

print(col_name)
print(sg.adata.obs.head())
cell_line_replicate
        cell_line replicate  dataset  total_counts cell_line_replicate
dataset                                                               
hepg2_1     hepg2         1  hepg2_1      499647.0             hepg2_1
hepg2_2     hepg2         2  hepg2_2      848447.0             hepg2_2
hffc6_1     hffc6         1  hffc6_1      761493.0             hffc6_1
hffc6_2     hffc6         2  hffc6_2      787967.0             hffc6_2
hffc6_3     hffc6         3  hffc6_3      614921.0             hffc6_3

The added column in col_name can then be used as the obs_col input to die_gene_test() as follows:

obs_col = col_name
obs_conditions = ['hffc6_3', 'hepg2_1']


die_table, die_results = sg.die_gene_test(kind='iso',
                                          obs_col=obs_col,
                                          obs_conditions=obs_conditions)

Exon skipping and intron retention

Swan can detect novel (unannotated) exon skipping and intron retention events.

To obtain a dataframe of novel exon skipping events, run the following code:

# returns a DataFrame of genes, transcripts, and specific edges in
# the SwanGraph with novel exon skipping events
es_df = sg.find_es_genes(verbose=True)
Testing each novel edge for exon skipping:   0%|          | 0/855 [00:00<?, ?it/s]

Analyzing 855 intronic edges for ES


Testing each novel edge for exon skipping: 100%|██████████| 855/855 [1:26:06<00:00,  6.13s/it]

Found 529 novel es events in 149 transcripts.
es_df.head()

To obtain a list of genes containing novel intron retention events, run the following code:

# returns a DataFrame of genes, transcripts, and specific edges in
# the SwanGraph with novel intron retaining events
ir_df = sg.find_ir_genes(verbose=True)
  0%|          | 0/1186 [00:00<?, ?it/s]
Testing each novel edge for intron retention:   0%|          | 0/1186 [00:00<?, ?it/s]

Analyzing 1186 exonic edges for IR
Testing each novel edge for intron retention: 100%|██████████| 1186/1186 [1:59:16<00:00,  5.97s/it]

Found 35 novel ir events in 27 transcripts.

You can pass gene IDs from es_df into gen_report() or plot_graph() to visualize where these exon-skipping events are in gene reports or gene summary graphs respectively.

PreviousGetting startedNextVisualization

Last updated