CanDI and DESeq2

Let’s say I want to look at changes in RNA expression across some cell lines in CCLE. DESeq2 is my preffered tool for doing differential expression analysis, unforutantely it’s written in R. CanDI makes it easy to format CCLE read counts data into the shape that DESeq2 expects.

import pandas as pd
import anndata as ad
import CanDI.candi as can

from CanDI.pipelines import diffexp

Object Instantiation

For this example I’m going to do differential expression analysis across male and female KRAS mutant cell lines. The cell below uses CanDI to generate the correct CellLineCluster objects for our purpose.

if type(can.data.mutations) != pd.DataFrame:
    can.data.load('mutations')

lung = can.Cancer("Lung Cancer", subtype = "NSCLC")
lung = can.CellLineCluster(lung.mutated("KRAS", variant = "Variant_Classification", item = "Missense_Mutation"))

lung_male = can.CellLineCluster(list(lung._info.loc[lung._info.sex == "Male",].index))
lung_female = can.CellLineCluster(list(lung._info.loc[lung._info.sex == "Female"].index))

Data Munging

The follow function takes two objects that we want to compare and automatically generates the counts and coldata matrices that DESeq2 needs to run. It’s typically a good idea to filter our genes/transcripts with consistently low counts prior to running DESeq2. This speeds up analysis and avoids issues related to read count scaling and multiple hypthothesis testing corrections. In this case we don’t care about different splicing of the same genes so I sum counts for duplicate indeces for all samples.

if type(can.data.rnaseq_reads) != pd.DataFrame:
    can.data.load('rnaseq_reads')

def make_counts_coldata(obj1, obj2, condition, factor1, factor2):
    
    counts1 = obj1.rnaseq_reads
    coldat1 = pd.Series(counts1.shape[1] * [factor1], index = counts1.columns, name = condition)
    
    counts2 = obj2.rnaseq_reads
    coldat2 = pd.Series(counts2.shape[1] * [factor2], index = counts2.columns, name = condition)
    
    #Concatenate Column Data
    coldat = pd.concat([coldat1, coldat2], axis = 0)
    #Concatenate read count data 
    counts_mat = pd.concat([counts1, counts2], axis = 1)
    #Sum duplicate indeces
    counts_mat = counts_mat.groupby(counts_mat.index).sum().astype(int)
    
    adata = ad.AnnData(
        counts_mat.T,
        obs = coldat.to_frame()
    )

    return adata

adata = make_counts_coldata(lung_male, lung_female, "sex", "male", "female")

Running pyDESeq2

results = diffexp.run_deseq(
    adata,
    design = 'sex',
    tested_level = 'male',
    ref_level = 'female'
)

Fitting size factors...
... done in 0.05 seconds.

Fitting dispersions...
... done in 1.39 seconds.

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

Fitting MAP dispersions...
... done in 1.42 seconds.

Fitting LFCs...
... done in 2.28 seconds.

Replacing 5676 outlier genes.

Fitting dispersions...
... done in 0.28 seconds.

Fitting MAP dispersions...
... done in 0.28 seconds.

Fitting LFCs...
... done in 0.36 seconds.

Running Wald tests...
... done in 1.64 seconds.

Log2 fold change & Wald test p-value: sex male vs female
              baseMean  log2FoldChange     lfcSE      stat    pvalue      padj
gene                                                                          
A1BG        200.696432        0.027944  0.671884  0.041591  0.966825  0.994527
A1BG-AS1    255.912479        0.202698  0.626175  0.323708  0.746159  0.950009
A1CF         12.057056       -1.367054  0.537503 -2.543342  0.010980  0.236614
A2M          18.685828        0.365012  0.619866  0.588857  0.555957  0.902738
A2M-AS1      49.283425       -0.266289  0.618898 -0.430264  0.667004  0.933027
...                ...             ...       ...       ...       ...       ...
ZYG11AP1      0.038949        0.057154  3.391294  0.016853  0.986554       NaN
ZYG11B     2200.470135       -0.200760  0.214134 -0.937546  0.348478  0.816844
ZYX       11155.922014        0.206356  0.377946  0.545994  0.585070  0.912754
ZZEF1      4400.173132       -0.483351  0.222239 -2.174911  0.029637  0.362384
ZZZ3       3532.326301       -0.212461  0.170959 -1.242756  0.213958  0.722301

[52443 rows x 6 columns]

Analyzing Results

Now we can read the results of the differential expression analysis back into our python enviroment and continue analysis as necessary. Looking at the genes with the lowest adjusted p-value we see that XIST is the most significant differentially expressed genes. XIST is a lncRNA responsible for X inactivation in females so this a good positive control for our analysis.

results.sort_values("padj").head()

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
gene
XIST	4041.971706	-7.145007	0.710899	-10.050660	9.125401e-24	2.692176e-19
RPS4Y1	6613.044427	8.723193	1.031580	8.456144	2.763634e-17	4.076636e-13
CEACAM5	10444.137420	-7.282857	0.892029	-8.164374	3.231066e-16	3.177430e-12
DDX3Y	1195.388509	7.642864	0.983916	7.767805	7.985781e-15	5.889913e-11
GJB1	58.891176	-5.623715	0.752347	-7.474896	7.726502e-14	4.558945e-10