Single Cell Comparative Network Reconstruction with scCNR

In this Notebook we use scCNR to reconstruct the same signaling network of the scMRA_Tutorial comparatively between two cell populations. Next to the wildtype populations that we analyzed in the scMRA_Tutorial we here have as well a simulated BRAF cell population. ScCNR tries to not only identify the network topology, but also identify the main differences in signaling between the wildtype and BRAF mutant population.

import scmra
import pandas as pd
import numpy as np
import seaborn as sns; sns.set_theme()
import matplotlib.pyplot as plt
from IPython.display import display, HTML
import networkx as nx

#read in data
dataWT = pd.read_csv("../data/scData_wildtype.tsv", sep='\t')
dataBRAF = pd.read_csv("../data/scData_BRAFmut.tsv", sep='\t')

Data Preparation

Similar to scMRA we need to calculate deviations of single cells to the population mean. Since we now consider several populations of single cells, we additionally need to assign group labels to map a cellID to its population.

#add a SC-ID
dataWT.index = ["CTR."+str(i) for i in range(dataWT.shape[0])]
dataBRAF.index = ["BRAFmut."+str(i) for i in range(dataBRAF.shape[0])]

#calculate input: deviatipon of protein from cell population mean
populations = [dataWT, dataBRAF]
#later we also have to access the populations specific rloc values by these labels
labels = ["Ctr", "BRAFmut"]
dataScaled, groupAnnot = scmra.prepare_data(populations, labels)

#get active protein deviations (frmo those we infer interaction strengths)
rglob = dataScaled.loc[[n for n in list(dataScaled.index) if ("Active" in n)]]
rglob.index = [i.replace("Active", "") for i in rglob.index]

#get total protein deviations: 'natural perturbations'
rtot = dataScaled.loc[[n for n in list(dataScaled.index) if ("Tot" in n)]]
rtot.index = [i.replace("Tot", "") for i in rtot.index]

display(rglob.head)

<bound method NDFrame.head of            CTR.0     CTR.1     CTR.2     CTR.3     CTR.4     CTR.5     CTR.6  \
Sos    -0.051155  0.085232  0.008321  0.073530  0.037336  0.077462  0.092302   
Ras    -0.076395  0.054095 -0.007728  0.025145  0.028869  0.058157  0.088900   
Raf1   -0.103680  0.069777  0.053183  0.041268  0.066881  0.064294  0.005852   
Mek    -0.033101 -0.152314 -0.050093 -0.128995  0.069491 -0.148370  0.166537   
Erk    -0.106257 -0.177075  0.030958 -0.113891  0.051923 -0.161771  0.092711   
P90Rsk -0.129504 -0.117604 -0.007665 -0.074405 -0.126830 -0.019424  0.120107   
PI3K    0.142959  0.021876 -0.143707  0.023346  0.029151  0.198616  0.008005   
Akt     0.116549 -0.098978 -0.143954 -0.084988  0.046098  0.104219  0.110371   
C3G     0.007717 -0.082415 -0.006071 -0.138128  0.056060 -0.083842 -0.024886   
Rap1   -0.005878 -0.069144  0.016994 -0.137141  0.027891 -0.114642  0.005731   
bRaf    0.072555 -0.084724 -0.030218 -0.107014  0.048830 -0.116599  0.097450   

           CTR.7     CTR.8     CTR.9  ...  BRAFmut.240  BRAFmut.241  \
Sos    -0.008603 -0.222532  0.126810  ...    -0.146622    -0.093920   
Ras     0.004750 -0.215548  0.117134  ...    -0.129316    -0.127258   
Raf1    0.066933 -0.259652  0.135834  ...    -0.123249    -0.081300   
Mek     0.058568  0.107128 -0.142194  ...     0.111068    -0.166231   
Erk     0.024951  0.221935 -0.127559  ...     0.102615     0.071023   
P90Rsk  0.015984  0.237644 -0.188018  ...     0.238576    -0.028478   
PI3K    0.087303 -0.074450 -0.020944  ...     0.052805     0.006293   
Akt    -0.176313  0.044074 -0.082037  ...     0.036224    -0.175243   
C3G     0.005292  0.115988  0.113998  ...     0.052111     0.107812   
Rap1    0.009414  0.123415  0.086880  ...     0.055672     0.117495   
bRaf    0.006757  0.086650  0.085719  ...     0.031120     0.183854   

        BRAFmut.242  BRAFmut.243  BRAFmut.244  BRAFmut.245  BRAFmut.246  \
Sos        0.041602     0.083655     0.090291    -0.179314     0.160474   
Ras        0.041264     0.118002     0.074339    -0.171745     0.168231   
Raf1       0.017540     0.185428     0.098197    -0.200488     0.163796   
Mek       -0.066659    -0.113333    -0.065334     0.011382     0.044878   
Erk        0.223207    -0.092159    -0.069101    -0.063405    -0.034418   
P90Rsk     0.026048    -0.060250    -0.011676     0.124937     0.005208   
PI3K       0.030364    -0.048787    -0.154532    -0.080240     0.048716   
Akt        0.104959    -0.075421    -0.116958     0.152433    -0.011654   
C3G        0.138902     0.092160    -0.031306    -0.015758    -0.005661   
Rap1       0.164348     0.098692    -0.030684     0.013867    -0.006914   
bRaf      -0.019847     0.222383    -0.149625    -0.095272     0.202139   

        BRAFmut.247  BRAFmut.248  BRAFmut.249  
Sos        0.366688     0.029405    -0.023816  
Ras        0.407325     0.011223    -0.040518  
Raf1       0.496338    -0.056441    -0.009805  
Mek       -0.043485    -0.000050     0.031329  
Erk       -0.059881    -0.025874     0.281526  
P90Rsk    -0.180146     0.017860     0.090534  
PI3K      -0.188240    -0.041037    -0.089324  
Akt        0.003887     0.010586    -0.002628  
C3G        0.099245    -0.083501    -0.032298  
Rap1       0.110646    -0.056336    -0.041837  
bRaf       0.052489    -0.135774     0.252340  

[11 rows x 500 columns]>

Define Network Sparcity

scCNR penalizes the number of non-equal edges in the network(additional to the number of edges in the network). Similar to the estimation of a suitable η-parameter, the θ-penalty on non-equal edges can also be estimated (see scMRA_Tutorial). Here we continue with parameters that we estimated for the method evaluation in the paper.

Reconstruct the network

The result of the scCNR ananlysis return not only one weighted adjacency matrix of the inferred network, but one matrix for each cell population. Depending on the θ-penalty we force edges to be similar between networks, and would ideally retain 'rloc' matrices, that are fairly similar and differ in the main interactions that explain hte main signaling differences between the populations. Here we display the differences between the wildtype and BRAFmut population and see that most differences are detected downstream of BRAF.

%%capture

#make MRA simulation
eta = 0.02
theta = 0.1
scd = scmra.ScData(rglob=rglob, rtot = rtot, group_annot= groupAnnot)
p = scmra.ScCnrProblem(scd, eta=eta, theta=theta) 
p.cpx.solve()
result = scmra.ScCnrResult(p)

#display difference between both adjacency matrices
display(result.rloc["Ctr"] - result.rloc["BRAFmut"])
differenceRloc = result.rloc["Ctr"] - result.rloc["BRAFmut"]

#diagonal rloc entries are per definition -1, replace those with 0 to plot network
A = differenceRloc.transpose()
A[A == -1] = 0
G = nx.from_pandas_adjacency(A, create_using=nx.DiGraph)

pos=nx.shell_layout(G)
nx.draw_networkx(G, with_labels=True, node_size = 750, pos = pos, node_color = "yellow")
#create edge labels
edgeLabels = {}
edgeColor = {}
for row_label, row_data in A.iterrows():
    for col_label, value in row_data.items():
        key = (col_label, row_label)
        if(np.abs(value) > 0):
            edgeLabels[key] = round(value,1)

nx.draw_networkx_edge_labels(
    G, pos,
    edge_labels=edgeLabels,
    font_color="black"
)
plt.show()

	Sos	Ras	Raf1	Mek	Erk	P90Rsk	PI3K	Akt	C3G	Rap1	bRaf
Sos	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
Ras	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
Raf1	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
Mek	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.791585
Erk	0.0	0.0	0.0	0.965511	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
P90Rsk	0.0	0.0	0.0	0.000000	0.701235	0.0	0.0	0.0	0.0	0.0	0.000000
PI3K	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
Akt	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
C3G	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
Rap1	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000
bRaf	0.0	0.0	0.0	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.000000