ManuscriptS4 - Network-to-network information transfer mapping. Assessing statistical significance by permuting FC topology between networks

Permute FC topology and see whether information transfer persists

Analysis for Supplementary Figure 4.

Master code for Ito et al., 2017

Takuya Ito (takuya.ito@rutgers.edu)



In [3]:

    
import sys
sys.path.append('utils/')
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import statsmodels.sandbox.stats.multicomp as mc
import multiprocessing as mp
%matplotlib inline
import os
os.environ['OMP_NUM_THREADS'] = str(1)
import warnings
warnings.filterwarnings('ignore')
from statsmodels.distributions.empirical_distribution import ECDF

import networkinformationtransfer as n2n


from matplotlib.colors import Normalize


class MidpointNormalize(Normalize):
    def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
        self.midpoint = midpoint
        Normalize.__init__(self, vmin, vmax, clip)

    def __call__(self, value, clip=None):
        # I'm ignoring masked values and all kinds of edge cases to make a
        # simple example...
        x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
        return np.ma.masked_array(np.interp(value, x, y))

0.0 Basic parameters



In [4]:

    
# Set basic parameters
datadir = '/projects2/ModalityControl2/data/'
runLength = 4648
subjNums = ['032', '033', '037', '038', '039', '045', 
            '013', '014', '016', '017', '018', '021', 
            '023', '024', '025', '026', '027', '031', 
            '035', '046', '042', '028', '048', '053', 
            '040', '049', '057', '062', '050', '030', '047', '034']

# Load in network array
networkdef = np.loadtxt(datadir + 'GlasserKKPartition/Glasser_KK_Partitions.csv', delimiter=',')


# Load in network keys (each network associated with a number in network array)
networkmappings = {'fpn':7, 'vis':1, 'smn':2, 'con':3, 'dmn':6, 'aud1':8, 'aud2':9, 'dan':11}
# Force aud2 key to be the same as aud1 (merging two auditory networks)
aud2_ind = np.where(networkdef==networkmappings['aud2'])[0]
networkdef[aud2_ind] = networkmappings['aud1']
# Redefine new network mappings with no aud1/aud2 distinction
networkmappings = {'fpn':7, 'vis':1, 'smn':2, 'con':3, 'dmn':6, 'aud':8, 'dan':11}

1.0 Run information transfer mapping procedure

1.1 First load in resting-state FC matrices for subjects



In [5]:

    
fcmat = {}
for subj in subjNums:
    fcmat[subj] = np.loadtxt(datadir + 'resultsMaster/MultRegConnRestFC_GlasserParcels/' + subj + '_multregconn_restfc.csv', delimiter=',')

1.2 Perform network-to-network information transfer mapping procedure using python module

1.2.1 First construct a wrapper to pass through a multiprocessing scheme



In [6]:

    
def informationTransferMappingWrapper((subj,fcmat)):
    """
    A wrapper so we can use multiprocessing to run subjects in parallel
    """
    out = n2n.networkToNetworkInformationTransferMapping(subj,fcmat,null=False)
    return out

1.2.2 Run using multiprocessing (parallel processing) to speed-up computation



In [7]:

    
inputs = []
for subj in subjNums:
    inputs.append((subj,fcmat[subj]))
    
pool = mp.Pool(processes=16)
results = pool.map_async(informationTransferMappingWrapper,inputs).get()
pool.close()
pool.join()

# Collect results
ruledims = ['logic','sensory','motor']
ite_matrix = {}
for ruledim in ruledims:
    ite_matrix[ruledim] = np.zeros((len(networkmappings),len(networkmappings),len(subjNums)))

scount = 0
for result in results:
    for ruledim in ruledims:
        ite_matrix[ruledim][:,:,scount] = result[ruledim]
        
    scount += 1

2.0 Calculate group average (True value)

Keep track of networks to matrix indices



In [8]:

    
# Create dictionary that reflects network ordering for matrix rows and columns
netkeys = {0:'vis',1:'smn',2:'con',3:'dmn',4:'fpn', 5:'aud', 6:'dan'}
num_networks=len(netkeys)



In [9]:

    
baseline = 0.0
avg_rho = {}
tstats = {}
pvals = {}
for ruledim in ruledims:
    avg_rho[ruledim] = np.zeros((num_networks,num_networks))
    tstats[ruledim] = np.zeros((num_networks,num_networks))
    for net1 in netkeys:
        
        for net2 in netkeys:
            # Skip if net1 and net2
            if net1==net2: 
                continue
                
            # Store results
            avg_rho[ruledim][net1,net2] = np.mean(ite_matrix[ruledim][net1,net2,:])
            tstats[ruledim][net1,net2] = stats.ttest_1samp(ite_matrix[ruledim][net1,net2,:],0)[0]

3.0 Load in null distribution for each network-to-network transfer

This was performed on a supercomputer due to computational expense. Code for these permutations can be found at: ./SupercomputerScripts/SFig4_FCPermutationTest



In [10]:

    
nPermutations = 1100
permdata = '/projects2/ModalityControl2/data/resultsMaster/ManuscriptS4_NetworkToNetworkITE_NullModel/PermutationTests_WithinNetPerm/'
cdf_val = {}
pmat = {}
for ruledim in ruledims:
    cdf_val[ruledim] = np.zeros((num_networks,num_networks))
    pmat[ruledim] = np.zeros((num_networks,num_networks))
    
    null_rhos = np.zeros((num_networks,num_networks,nPermutations))
    print 'Computing p-stats for', ruledim
    
    for i in range(nPermutations):
        null_rhos[:,:,i] = np.loadtxt(permdata + 'NetworkITE_NullModel_' + ruledim + '_Permutation' + str(i) + '.csv',
                                    delimiter=',')

    for net1 in netkeys:

        for net2 in netkeys:
            # Skip if net1 and net2
            if net1==net2: 
                continue

            # Store results
            # Compute ecdf for each net1 - net2 config
            ecdf = ECDF(null_rhos[net1,net2,:])
            # Now find where in the empirical distribution the true rho value lies
            cdf_val[ruledim][net1,net2] = ecdf(avg_rho[ruledim][net1,net2])
            # Since p-values are inverted... and we're only doing 1-sided tests
            pmat[ruledim][net1,net2] = 1.0 - ecdf(avg_rho[ruledim][net1,net2])









    



Computing p-stats for logic
Computing p-stats for sensory
Computing p-stats for motor

3.1 Perform multiple comparisons (using false discovery rate)



In [11]:

    
# Compute group stats
baseline = 0.0
triu_indices = np.triu_indices(len(networkmappings),k=1)
tril_indices = np.tril_indices(len(networkmappings),k=-1)

qmat = {}
for ruledim in ruledims:

    
    qmat[ruledim] = np.zeros((num_networks,num_networks))
    
    tmpq = []
    tmpq.extend(pmat[ruledim][triu_indices])
    tmpq.extend(pmat[ruledim][tril_indices])
    tmpq = mc.fdrcorrection0(tmpq)[1]
    qmat[ruledim][triu_indices] = tmpq[0:len(triu_indices[0])]
    qmat[ruledim][tril_indices] = tmpq[len(triu_indices[0]):]

3.2 Plot results



In [12]:

    
for ruledim in ruledims:
    plt.figure(figsize=(12,10))
    # First visualize unthresholded results
    plt.subplot(121)
    plt.title('NetworkToNetwork Information Transfer Mapping\n' + ruledim + ' domain', fontsize=14, y=1.04)
    mat = avg_rho[ruledim]
    np.fill_diagonal(mat,0)
    norm = MidpointNormalize(midpoint=0)
    plt.imshow(mat, origin='lower',norm=norm, vmin=0, cmap='seismic', interpolation='none')
    plt.xticks(netkeys.keys(),netkeys.values())
    plt.yticks(netkeys.keys(),netkeys.values())
    plt.ylabel('Source Network',fontsize=16)
    plt.xlabel('Target Network',fontsize=16)
    plt.colorbar(fraction=.046)
    
    # Next visualize thresholded results (after multiple comparisons)
    plt.subplot(122)
    plt.title('NetworkToNetwork Information Transfer Mapping\n' + ruledim + ' domain' , fontsize=14, y=1.04)
    mat = avg_rho[ruledim]
    thresh = qmat[ruledim] < 0.05
    # Threshold using q < 0.05
    mat = np.multiply(mat,thresh)
    np.fill_diagonal(mat,0)
    norm = MidpointNormalize(midpoint=0)
    plt.imshow(mat, origin='lower',norm=norm, vmin=0, cmap='seismic', interpolation='none')
    plt.xticks(netkeys.keys(),netkeys.values())
    plt.yticks(netkeys.keys(),netkeys.values())
    plt.ylabel('Source Network',fontsize=16)
    plt.xlabel('Target Network',fontsize=16)
    plt.colorbar(fraction=.046)
    plt.tight_layout()
#     plt.savefig('NetworkToNetworkITE_' + ruledim + '_PermutationTestResults.pdf')

Statistical testing with FWE permutation testing

4.0 Load in null distribution for each network-to-network transfer

4.1 Perform multiple comparisons (using FWE testing)



In [13]:

    
ruledims = ['logic', 'sensory', 'motor']
nPermutations = 1000
indices = np.ones((num_networks,num_networks))
np.fill_diagonal(indices,0)
flatten_ind = np.where(indices==1)

permdata = '/projects2/ModalityControl2/data/resultsMaster/ManuscriptS4_NetworkToNetworkITE_NullModel/PermutationTests_WithinNetPerm/'
pfwe = {}
for ruledim in ruledims:
    pfwe[ruledim] = np.zeros((num_networks,num_networks))
    
    print 'Computing p-stats for', ruledim
    maxTs = []
    for i in range(100,nPermutations+100):
        mat = np.loadtxt(permdata + 'NetworkITE_NullModel_' + ruledim + '_Permutation' + str(i) + '_Tstat.csv',
                                    delimiter=',')
        
        # Obtain max t-value 
        maxTs.append(np.max(mat[flatten_ind[0], flatten_ind[1]]))
    
    ecdf = ECDF(np.asarray(maxTs))

    for net1 in netkeys:

        for net2 in netkeys:
            # Skip if net1 and net2
            if net1==net2: 
                continue

            # Now find where in the empirical distribution the true rho value lies
            # Since p-values are inverted... and we're only doing 1-sided tests
            pfwe[ruledim][net1,net2] = 1.0 - ecdf(tstats[ruledim][net1,net2])









    



Computing p-stats for logic
Computing p-stats for sensory
Computing p-stats for motor

2.2 Plot results



In [15]:

    
for ruledim in ruledims:
    plt.figure(figsize=(12,10))
    # First visualize unthresholded results
    plt.subplot(121)
    plt.title('NetworkToNetwork Information Transfer Mapping\n' + ruledim + ' domain', fontsize=14, y=1.04)
    mat = avg_rho[ruledim]
    np.fill_diagonal(mat,0)
    norm = MidpointNormalize(midpoint=0)
    plt.imshow(mat, origin='lower',norm=norm, vmin=0, cmap='seismic', interpolation='none')
    plt.xticks(netkeys.keys(),netkeys.values())
    plt.yticks(netkeys.keys(),netkeys.values())
    plt.ylabel('Source Network',fontsize=16)
    plt.xlabel('Target Network',fontsize=16)
    plt.colorbar(fraction=.046)
    
    # Next visualize thresholded results (after multiple comparisons)
    plt.subplot(122)
    plt.title('NetworkToNetwork Information Transfer Mapping\n' + ruledim + ' domain' , fontsize=14, y=1.04)
    mat = avg_rho[ruledim]
    thresh = pfwe[ruledim] < 0.05
    # Threshold using q < 0.05
    mat = np.multiply(mat,thresh)
    np.fill_diagonal(mat,0)
    norm = MidpointNormalize(midpoint=0)
    plt.imshow(mat, origin='lower',norm=norm, vmin=0, cmap='seismic', interpolation='none')
    plt.xticks(netkeys.keys(),netkeys.values())
    plt.yticks(netkeys.keys(),netkeys.values())
    plt.ylabel('Source Network',fontsize=16)
    plt.xlabel('Target Network',fontsize=16)
    plt.colorbar(fraction=.046)
    plt.tight_layout()
    plt.savefig('NetworkToNetworkITE_' + ruledim + '_PermutationTestResults.pdf')