# Analysis

Swan has several analysis options to use.

* [Differential expression tests](#de)
* [Isoform switching / differential isoform expression](#is)
* [Combining metadata columns](#multi_gb)
* [Exon skipping and intron retention](#es_ir)
* [More differential expression](#more_de)

```python
import swan_vis as swan

sg = swan.read('data/swan.p')
```

```
Read in graph from data/swan.p
```

## Differential expression tests <a href="#de" id="de"></a>

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](https://github.com/owkin/PyDESeq2) or Scanpy's [`rank_genes_groups`](https://scanpy.readthedocs.io/en/stable/generated/scanpy.tl.rank_genes_groups.html) test, which both support the AnnData format for simple compatibility with Swan's AnnData expression representation.

### Using PyDESeq2 <a href="#de" id="de"></a>

Please read the [PyDESeq2 documentation](https://pydeseq2.readthedocs.io/en/latest/) 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.

```python
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
```

|                   | baseMean | log2FoldChange | lfcSE    | stat      | pvalue   | padj     |
| ----------------- | -------- | -------------- | -------- | --------- | -------- | -------- |
| tid               |          |                |          |           |          |          |
| TALONT000296400   | 3.497574 | 2.416436       | 1.472616 | 1.640914  | 0.100815 | 0.193840 |
| ENST00000591581.1 | 0.162115 | 0.624956       | 5.002765 | 0.124922  | 0.900585 | NaN      |
| ENST00000546893.5 | 8.173445 | -0.053334      | 0.793179 | -0.067241 | 0.946390 | 0.968527 |
| ENST00000537289.1 | 5.140216 | -1.248175      | 0.978927 | -1.275044 | 0.202294 | 0.328286 |
| ENST00000258382.9 | 2.012104 | -0.011819      | 1.493866 | -0.007912 | 0.993687 | NaN      |
| ...               | ...      | ...            | ...      | ...       | ...      | ...      |
| ENST00000506914.1 | 1.022692 | 0.947588       | 2.272813 | 0.416923  | 0.676735 | NaN      |
| ENST00000571080.1 | 0.473263 | -2.824409      | 3.426473 | -0.824291 | 0.409774 | NaN      |
| ENST00000378615.7 | 1.334873 | -1.199810      | 1.790781 | -0.669993 | 0.502862 | NaN      |
| ENST00000409586.7 | 0.338615 | 1.464364       | 4.023896 | 0.363917  | 0.715920 | NaN      |
| ENST00000370278.3 | 1.328549 | -1.201246      | 1.962358 | -0.612144 | 0.540443 | NaN      |

34814 rows × 6 columns

```python
df.head()
```

|                   | baseMean | log2FoldChange | lfcSE    | stat      | pvalue   | padj     |
| ----------------- | -------- | -------------- | -------- | --------- | -------- | -------- |
| tid               |          |                |          |           |          |          |
| TALONT000296400   | 3.497574 | 2.416436       | 1.472616 | 1.640914  | 0.100815 | 0.193840 |
| ENST00000591581.1 | 0.162115 | 0.624956       | 5.002765 | 0.124922  | 0.900585 | NaN      |
| ENST00000546893.5 | 8.173445 | -0.053334      | 0.793179 | -0.067241 | 0.946390 | 0.968527 |
| ENST00000537289.1 | 5.140216 | -1.248175      | 0.978927 | -1.275044 | 0.202294 | 0.328286 |
| ENST00000258382.9 | 2.012104 | -0.011819      | 1.493866 | -0.007912 | 0.993687 | NaN      |

## Isoform switching / Differential isoform expression testing <a href="#is" id="is"></a>

Isoform switching / differential isoform expression (DIE) testing is implemented according to the strategy in [Joglekar et. al., 2021](https://www.nature.com/articles/s41467-020-20343-5). 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.

```python
# look at valid metadata options
sg.adata.obs
```

|          | cell\_line | replicate | dataset  | total\_counts |
| -------- | ---------- | --------- | -------- | ------------- |
| dataset  |            |           |          |               |
| hepg2\_1 | hepg2      | 1         | hepg2\_1 | 499647.0      |
| hepg2\_2 | hepg2      | 2         | hepg2\_2 | 848447.0      |
| hffc6\_1 | hffc6      | 1         | hffc6\_1 | 761493.0      |
| hffc6\_2 | hffc6      | 2         | hffc6\_2 | 787967.0      |
| hffc6\_3 | hffc6      | 3         | hffc6\_3 | 614921.0      |

```python
# 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](https://www.nature.com/articles/s41467-020-20343-5).

```python
die_table.head()
```

|   | gid             | p\_val   | dpi       | pos\_iso\_1       | pos\_iso\_2       | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1       | neg\_iso\_2       | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname   |
| - | --------------- | -------- | --------- | ----------------- | ----------------- | ---------------- | ---------------- | ----------------- | ----------------- | ---------------- | ---------------- | ----------- | ------- |
| 0 | ENSG00000130175 | 0.206667 | 39.264992 | TALONT000296399   | NaN               | 39.264992        | NaN              | TALONT000296400   | ENST00000589838.5 | -20.116056       | -10.638298       | 0.469420    | PRKCSH  |
| 1 | ENSG00000130202 | 0.000367 | 11.893251 | ENST00000252485.8 | ENST00000252483.9 | 9.684967         | 2.208281         | TALONT000406668   | ENST00000591581.1 | -11.473083       | -0.420168        | 0.003560    | NECTIN2 |
| 2 | ENSG00000111371 | 0.680435 | 9.401713  | ENST00000398637.9 | ENST00000439706.5 | 8.547012         | 0.854701         | ENST00000546893.5 | NaN               | -9.401711        | NaN              | 0.886243    | SLC38A1 |
| 3 | ENSG00000181924 | 0.028195 | 9.452934  | ENST00000537289.1 | ENST00000545127.1 | 7.298603         | 2.154328         | ENST00000355693.4 | NaN               | -9.452934        | NaN              | 0.122619    | COA4    |
| 4 | ENSG00000163468 | 0.255788 | 0.680048  | ENST00000295688.7 | TALONT000476055   | 0.366966         | 0.277693         | ENST00000368259.6 | ENST00000489870.1 | -0.568503        | -0.111545        | 0.525066    | CCT3    |

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:

```python
# 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)
```

|   | gid             | p\_val   | dpi       | pos\_iso\_1       | pos\_iso\_2       | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1     | neg\_iso\_2       | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname   |
| - | --------------- | -------- | --------- | ----------------- | ----------------- | ---------------- | ---------------- | --------------- | ----------------- | ---------------- | ---------------- | ----------- | ------- |
| 0 | ENSG00000130175 | 0.206667 | 39.264992 | TALONT000296399   | NaN               | 39.264992        | NaN              | TALONT000296400 | ENST00000589838.5 | -20.116056       | -10.638298       | 0.46942     | PRKCSH  |
| 1 | ENSG00000130202 | 0.000367 | 11.893251 | ENST00000252485.8 | ENST00000252483.9 | 9.684967         | 2.208281         | TALONT000406668 | ENST00000591581.1 | -11.473083       | -0.420168        | 0.00356     | NECTIN2 |

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

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

|    | gid             | p\_val        | dpi       | pos\_iso\_1       | pos\_iso\_2       | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1       | neg\_iso\_2       | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val   | gname   |
| -- | --------------- | ------------- | --------- | ----------------- | ----------------- | ---------------- | ---------------- | ----------------- | ----------------- | ---------------- | ---------------- | ------------- | ------- |
| 1  | ENSG00000130202 | 3.674234e-04  | 11.893251 | ENST00000252485.8 | ENST00000252483.9 | 9.684967         | 2.208281         | TALONT000406668   | ENST00000591581.1 | -11.473083       | -0.420168        | 3.560035e-03  | NECTIN2 |
| 8  | ENSG00000025156 | 1.857411e-03  | 35.714287 | ENST00000452194.5 | ENST00000465214.2 | 25.714287        | 10.000000        | ENST00000368455.8 | NaN               | -35.714287       | NaN              | 1.405048e-02  | HSF2    |
| 20 | ENSG00000105254 | 4.779408e-38  | 25.419458 | ENST00000585910.5 | TALONT000366329   | 20.394853        | 4.702971         | ENST00000221855.7 | ENST00000589996.5 | -25.016304       | -0.403154        | 7.343219e-36  | TBCB    |
| 23 | ENSG00000105379 | 1.802782e-307 | 81.136333 | ENST00000354232.8 | NaN               | 81.136333        | NaN              | ENST00000309244.8 | ENST00000596253.1 | -80.857935       | -0.278394        | 1.551113e-304 | ETFB    |
| 28 | ENSG00000148180 | 0.000000e+00  | 89.319039 | ENST00000373818.8 | TALONT000419680   | 85.245850        | 4.073189         | TALONT000418752   | ENST00000373808.8 | -68.603180       | -7.768958        | 0.000000e+00  | GSN     |

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`.

```python
# 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]
```

|   | gid             | p\_val   | dpi       | pos\_iso\_1        | pos\_iso\_2        | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1        | neg\_iso\_2        | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname   |
| - | --------------- | -------- | --------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------------ | ------------------ | ---------------- | ---------------- | ----------- | ------- |
| 0 | ENSG00000000419 | 0.249561 | 10.002187 | ENSG00000000419\_2 | ENSG00000000419\_1 | 9.906052         | 0.096133         | ENSG00000000419\_3 | ENSG00000000419\_4 | -7.218704        | -2.783483        | 0.536906    | DPM1    |
| 1 | ENSG00000001461 | 0.852914 | 9.093739  | ENSG00000001461\_5 | NaN                | 9.093739         | NaN              | ENSG00000001461\_1 | ENSG00000001461\_3 | -5.543442        | -1.775148        | 1.000000    | NIPAL3  |
| 2 | ENSG00000001497 | 0.891496 | 3.846153  | ENSG00000001497\_1 | NaN                | 3.846146         | NaN              | ENSG00000001497\_2 | NaN                | -3.846153        | NaN              | 1.000000    | LAS1L   |
| 3 | ENSG00000001630 | 0.286866 | 2.580168  | ENSG00000001630\_3 | NaN                | 2.580168         | NaN              | ENSG00000001630\_2 | ENSG00000001630\_1 | -1.549232        | -1.030928        | 0.582636    | CYP51A1 |
| 4 | ENSG00000002330 | 0.184635 | 12.817431 | ENSG00000002330\_2 | ENSG00000002330\_1 | 11.868484        | 0.948944         | ENSG00000002330\_4 | ENSG00000002330\_3 | -12.603584       | -0.213847        | 0.454053    | BAD     |

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

```python
# 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)
```

|   | gid             | p\_val   | dpi       | pos\_iso\_1        | pos\_iso\_2        | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1        | neg\_iso\_2        | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname  |
| - | --------------- | -------- | --------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------------ | ------------------ | ---------------- | ---------------- | ----------- | ------ |
| 0 | ENSG00000000419 | 0.249561 | 10.002187 | ENSG00000000419\_2 | ENSG00000000419\_1 | 9.906052         | 0.096133         | ENSG00000000419\_3 | ENSG00000000419\_4 | -7.218704        | -2.783483        | 0.536906    | DPM1   |
| 1 | ENSG00000001461 | 0.852914 | 9.093739  | ENSG00000001461\_5 | NaN                | 9.093739         | NaN              | ENSG00000001461\_1 | ENSG00000001461\_3 | -5.543442        | -1.775148        | 1.000000    | NIPAL3 |

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

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

|    | gid             | p\_val       | dpi       | pos\_iso\_1        | pos\_iso\_2        | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1        | neg\_iso\_2        | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val  | gname   |
| -- | --------------- | ------------ | --------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------ | ------- |
| 10 | ENSG00000003402 | 7.338098e-15 | 84.848486 | ENSG00000003402\_3 | ENSG00000003402\_7 | 77.272728        | 7.575758         | ENSG00000003402\_1 | ENSG00000003402\_5 | -31.717173       | -22.222223       | 4.390279e-13 | CFLAR   |
| 11 | ENSG00000003436 | 3.812087e-06 | 34.072281 | ENSG00000003436\_1 | ENSG00000003436\_3 | 28.073771        | 4.564083         | ENSG00000003436\_4 | NaN                | -34.072281       | NaN              | 7.064168e-05 | TFPI    |
| 16 | ENSG00000004487 | 1.135407e-03 | 75.714287 | ENSG00000004487\_1 | NaN                | 75.714287        | NaN              | ENSG00000004487\_2 | NaN                | -75.714285       | NaN              | 1.020405e-02 | KDM1A   |
| 33 | ENSG00000006282 | 4.544048e-03 | 19.819595 | ENSG00000006282\_2 | ENSG00000006282\_1 | 13.679245        | 6.140350         | ENSG00000006282\_3 | NaN                | -19.819595       | NaN              | 3.236475e-02 | SPATA20 |
| 43 | ENSG00000007376 | 1.256379e-04 | 27.898551 | ENSG00000007376\_3 | ENSG00000007376\_1 | 23.454107        | 4.444445         | ENSG00000007376\_5 | ENSG00000007376\_7 | -14.130434       | -7.512077        | 1.525134e-03 | RPUSD1  |

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()`.

```python
# 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]
```

|   | gid             | p\_val   | dpi       | pos\_iso\_1        | pos\_iso\_2        | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1        | neg\_iso\_2        | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname   |
| - | --------------- | -------- | --------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------------ | ------------------ | ---------------- | ---------------- | ----------- | ------- |
| 0 | ENSG00000000419 | 1.000000 | 0.096133  | ENSG00000000419\_2 | NaN                | 0.096133         | NaN              | ENSG00000000419\_1 | NaN                | -0.096130        | NaN              | 1.000000    | DPM1    |
| 1 | ENSG00000001461 | 0.852914 | 9.093739  | ENSG00000001461\_3 | NaN                | 9.093739         | NaN              | ENSG00000001461\_4 | ENSG00000001461\_1 | -5.543442        | -1.775148        | 1.000000    | NIPAL3  |
| 2 | ENSG00000001630 | 0.286866 | 2.580168  | ENSG00000001630\_2 | NaN                | 2.580168         | NaN              | ENSG00000001630\_1 | ENSG00000001630\_3 | -1.549232        | -1.030928        | 0.618077    | CYP51A1 |
| 3 | ENSG00000002330 | 0.184635 | 12.817431 | ENSG00000002330\_1 | ENSG00000002330\_4 | 11.868484        | 0.948944         | ENSG00000002330\_2 | ENSG00000002330\_3 | -12.603584       | -0.213847        | 0.480860    | BAD     |
| 4 | ENSG00000002549 | 0.679148 | 0.694543  | ENSG00000002549\_2 | NaN                | 0.694542         | NaN              | ENSG00000002549\_1 | NaN                | -0.694543        | NaN              | 0.926889    | LAP3    |

```python
# 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)
```

|   | gid             | p\_val   | dpi      | pos\_iso\_1        | pos\_iso\_2 | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1        | neg\_iso\_2        | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname  |
| - | --------------- | -------- | -------- | ------------------ | ----------- | ---------------- | ---------------- | ------------------ | ------------------ | ---------------- | ---------------- | ----------- | ------ |
| 0 | ENSG00000000419 | 1.000000 | 0.096133 | ENSG00000000419\_2 | NaN         | 0.096133         | NaN              | ENSG00000000419\_1 | NaN                | -0.096130        | NaN              | 1.0         | DPM1   |
| 1 | ENSG00000001461 | 0.852914 | 9.093739 | ENSG00000001461\_3 | NaN         | 9.093739         | NaN              | ENSG00000001461\_4 | ENSG00000001461\_1 | -5.543442        | -1.775148        | 1.0         | NIPAL3 |

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

|    | gid             | p\_val       | dpi       | pos\_iso\_1        | pos\_iso\_2        | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1        | neg\_iso\_2        | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val  | gname   |
| -- | --------------- | ------------ | --------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------------ | ------------------ | ---------------- | ---------------- | ------------ | ------- |
| 9  | ENSG00000003402 | 7.338098e-15 | 84.848486 | ENSG00000003402\_5 | ENSG00000003402\_6 | 77.272728        | 7.575758         | ENSG00000003402\_9 | ENSG00000003402\_7 | -31.717173       | -22.222223       | 5.296535e-13 | CFLAR   |
| 10 | ENSG00000003436 | 2.772096e-07 | 32.637854 | ENSG00000003436\_1 | ENSG00000003436\_4 | 22.082711        | 10.555142        | ENSG00000003436\_2 | ENSG00000003436\_5 | -30.435921       | -1.818182        | 7.003007e-06 | TFPI    |
| 14 | ENSG00000004487 | 1.135407e-03 | 75.714287 | ENSG00000004487\_2 | NaN                | 75.714287        | NaN              | ENSG00000004487\_1 | NaN                | -75.714285       | NaN              | 1.141621e-02 | KDM1A   |
| 32 | ENSG00000006282 | 6.858641e-03 | 19.819595 | ENSG00000006282\_2 | NaN                | 19.819595        | NaN              | ENSG00000006282\_1 | NaN                | -19.819595       | NaN              | 4.915014e-02 | SPATA20 |
| 60 | ENSG00000010278 | 3.191221e-22 | 39.024387 | ENSG00000010278\_2 | NaN                | 39.024387        | NaN              | ENSG00000010278\_3 | ENSG00000010278\_1 | -38.605976       | -0.418410        | 3.486194e-20 | CD9     |

Finally, for transcriptomes generated with [Cerberus](https://github.com/mortazavilab/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.

```python
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
```

|   | gid               | p\_val       | dpi       | pos\_iso\_1          | pos\_iso\_2          | pos\_iso\_1\_dpi | pos\_iso\_2\_dpi | neg\_iso\_1          | neg\_iso\_2 | neg\_iso\_1\_dpi | neg\_iso\_2\_dpi | adj\_p\_val | gname             |
| - | ----------------- | ------------ | --------- | -------------------- | -------------------- | ---------------- | ---------------- | -------------------- | ----------- | ---------------- | ---------------- | ----------- | ----------------- |
| 0 | ENCODEMG000055991 | 4.100695e-01 | 6.071428  | ENCODEMG000055991\_2 | ENCODEMG000055991\_3 | 3.571428         | 2.500000         | ENCODEMG000055991\_1 | NaN         | -6.071426        | NaN              | 0.707992    | ENCODEMG000055991 |
| 1 | ENCODEMG000055998 | 6.190162e-01 | 9.722223  | ENCODEMG000055998\_2 | NaN                  | 9.722223         | NaN              | ENCODEMG000055998\_1 | NaN         | -9.722218        | NaN              | 0.836200    | ENCODEMG000055998 |
| 2 | ENCODEMG000056718 | 8.081718e-07 | 6.289159  | NaN                  | NaN                  | 2.329770         | 1.953835         | NaN                  | NaN         | -3.339438        | -2.949721        | 0.000020    | ENCODEMG000056718 |
| 3 | ENCODEMG000056804 | 2.107047e-03 | 27.863049 | ENCODEMG000056804\_1 | NaN                  | 27.863047        | NaN              | ENCODEMG000056804\_2 | NaN         | -27.863049       | NaN              | 0.021510    | ENCODEMG000056804 |
| 4 | ENCODEMG000063411 | 7.699701e-01 | 12.500000 | ENCODEMG000063411\_1 | NaN                  | 12.500000        | NaN              | ENCODEMG000063411\_2 | NaN         | -12.500000       | NaN              | 0.919808    | ENCODEMG000063411 |

## Combining metadata columns <a href="#multi_gb" id="multi_gb"></a>

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?

```python
sg.adata.obs.head()
```

|          | cell\_line | replicate | dataset  | total\_counts |
| -------- | ---------- | --------- | -------- | ------------- |
| dataset  |            |           |          |               |
| hepg2\_1 | hepg2      | 1         | hepg2\_1 | 499647.0      |
| hepg2\_2 | hepg2      | 2         | hepg2\_2 | 848447.0      |
| hffc6\_1 | hffc6      | 1         | hffc6\_1 | 761493.0      |
| hffc6\_2 | hffc6      | 2         | hffc6\_2 | 787967.0      |
| hffc6\_3 | hffc6      | 3         | hffc6\_3 | 614921.0      |

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:

```python
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:

```python
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 <a href="#es_ir" id="es_ir"></a>

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

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

```python
# 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.
```

```python
es_df.head()
```

|   | gid                | tid             | edge\_id |
| - | ------------------ | --------------- | -------- |
| 0 | ENSG00000157916.19 | TALONT000218256 | 952616   |
| 0 | ENSG00000122406.13 | TALONT000425229 | 952716   |
| 0 | ENSG00000224093.5  | TALONT000434035 | 952720   |
| 0 | ENSG00000224093.5  | TALONT000434035 | 952720   |
| 0 | ENSG00000224093.5  | TALONT000434035 | 952720   |

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

```python
# 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][A
Testing each novel edge for intron retention:   0%|          | 0/1186 [00:00<?, ?it/s][A

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

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.
