shell:
    'mkdir -p {output.probtrack_dir} && singularity exec -e --nv {params.container} probtrackx2 --samples={params.bedpost_merged} --mask={input.mask} --seed={input.seed_res} --targetmasks={input.target_txt} --seedref={input.seed_res} --nsamples={params.nsamples} --dir={output.probtrack_dir} {params.probtrack_opts} -V 1  &> {log}'

shell:
    'c3d {input.label_vol} -replace 110 0 -o {output.label_vol_removepir110} &> {log}'

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import nibabel as nib
from functools import reduce
import os


#Set the name of the seed you want to create the Gephi input for
seeds = snakemake.params.seeds_py

#Set the k parcellation you want to create the Gephi input for
#this_k = 3
cluster_range = range(2,snakemake.params.max_k+1)

#Set the path to the location of the directory containing the connmapGroup directories for each seed and each k 
#(i.e. connmapClusters dir from diffparc located in /scratch/nchris5/HCP_template_piriform_segmentations/diffparc/results)
path = snakemake.params.connmapClusters_dir


#Read in the Glasser_2016_Table with the 'SectionsBOLD' column manually included, which indicates which of the 22 overarching regions each of the 180 regions belongs to; 
#and the 'ColourRGB' column manually included, which indicates the RGB for that specific SectionsBOLD
GlasserTable = pd.read_csv(snakemake.params.GlasserTable_22regions_colours, header=0, usecols=range(1,8))
GlasserTable.index = range(1,snakemake.params.num_targets+1)
print('Original GlasserTable : \n', GlasserTable)
#Sort by the section values from 1-->22
GlasserTable_sorted_PrimarySection = GlasserTable.sort_values(by='SectionsBOLD', ascending=True)
#Change the actual index to 1-->179
GlasserTable_sorted_PrimarySection.index = range(1,snakemake.params.num_targets+1)
#Insert a column that describes the sorted dataframe new index from 1-->180
GlasserTable_sorted_PrimarySection.insert(loc=1, column='Sorted_PrimarySection_Index', value=range(1,snakemake.params.num_targets+1))
#Create the list of target ID's with leading zeros so all ID's have 3 digits (i.e. 001, 010, 100), as for gephi to read in proper order
GlasserTable_sorted_PrimarySection_Index_ThreeDigits = [str(indID).zfill(3) for indID in GlasserTable_sorted_PrimarySection.Sorted_PrimarySection_Index]
print('GlasserTable_sorted_PrimarySection : \n', GlasserTable_sorted_PrimarySection)

print('\n ------ Creating Nodes df ------ \n ')
#Includes the Id, Label, color, and SectionsBOLD for each of the 179 or 358 Glasser targets (001-179 or 001-358)
Targets = pd.DataFrame({'Id': GlasserTable_sorted_PrimarySection_Index_ThreeDigits, 'Label': GlasserTable_sorted_PrimarySection.AreaName, 'color': GlasserTable_sorted_PrimarySection.ColourRGB, 'Modularity Class': GlasserTable_sorted_PrimarySection.SectionsBOLD}, index=range(1,snakemake.params.num_targets+1))
print('Targetsdf : \n', Targets)

for seed in seeds:
    for i,this_k in enumerate(cluster_range): #For 0,2 1,3 2,4
        #Set the RBG colours for each of the clusters within this_k to the same colour that is present in the clustering output in itksnap (i.e. diffparc spectral_clustering output)
        if this_k==2:
            this_k_ColourRGB = ['255,0,0','0,255,0'] #Cluster01=Red, #Cluster02=Green
        elif this_k==3:
            this_k_ColourRGB = ['255,0,0','0,255,0','0,0,255'] #Cluster01=Red, #Cluster02=Green, #Cluster03=Blue
        elif this_k==4:
            this_k_ColourRGB = ['255,0,0','0,255,0','0,0,255', '255,255,0'] #Cluster01=Red, #Cluster02=Green, #Cluster03=Blue #Yellow
        #Create the Id's and labels for the individual clusters  for and get the basename for the individual clusters (e.g. Cluster01, Cluster02, Cluster03), assigning Id values above 179 or 358 (since 179 or 358 targets) e.g. 180=Cluster01, 181=Cluster02, 182=Cluster03
        IdlistClusters = list()
        labellist = list()
        for i,this_label in enumerate(range(1,this_k+1)):
            IdlistClusters.append(snakemake.params.num_targets+1+i)
            labellist.append('Cluster{this_label:02d}'.format(this_label=this_label))
        #Set the Modularity Class for each individual cluster to be distinct (i.e. Cluster01=23, Cluster02=24, Cluster03=25)
        ModularityClassClusters = list(range(23,this_k+23))
        #Create the dataframe for the clusters to add to the list of targets
        individualclusterstoadd = pd.DataFrame({'Id': IdlistClusters, 'Label': labellist, 'color': this_k_ColourRGB, 'Modularity Class': ModularityClassClusters}, index=range(snakemake.params.num_targets,snakemake.params.num_targets+this_k))
        print('\n individualclusterstoadd : \n', individualclusterstoadd)

        #Append the cluster labels dataframe to the bottom of the targets list dataframe to yield the final Nodesdf
        Nodesdf = Targets.append(individualclusterstoadd)
        print('\n Nodesdf : \n ', Nodesdf)




        print('\n ------ Creating Edges df ------ \n ')
        print(print('------ Creating Edges df SourceId and SourceName columns ------ \n '))

        #Multiply the targets df by number of clusters (i.e. first 180 will represent cluster 1, next 180 will represent cluster 2, etc.)
        EdgesSourceId = np.zeros((snakemake.params.num_targets,this_k), dtype=int)
        EdgesSourceNames = list()
        for i,this_label in enumerate(range(1,this_k+1)):
            EdgesSourceId[:,i] = this_label+snakemake.params.num_targets
            tempSourceNames = list()
            for j in range(snakemake.params.num_targets):
                tempSourceNames.append('Cluster{this_label:02d}'.format(this_label=this_label))
            EdgesSourceNames.append(tempSourceNames)
        #Use functools reduce to remove iterative lists into a single list i.e. [179xCluster01], [179xCluster02], [179xCluster03] to [179xCluster01, 179xCluster02, 179xCluster03]
        EdgesSourceNames = reduce(lambda x,y: x+y, EdgesSourceNames)
        #Reshape the edges source col current with shape (179,numClusters) to 179*numClusters with order 'F' (i.e. 179x'181' --> 179x'182' --> 179x'183')
        EdgesSourceId = np.reshape(EdgesSourceId, (snakemake.params.num_targets*this_label), order='F')

        TargetsMultiplied = pd.concat([Targets]*this_k, ignore_index=True)
        print('TargetsMultiplied Shape = ', TargetsMultiplied.shape)
        print('TargetsMultiplied = \n', TargetsMultiplied)

        print('\n ------ Creating Edges df Columns:   Source, SLabel, Target, TLabel, ConnectivityScore, ConnectivityScoreNormalized, Rank_in_Cluster, Weight ------ \n ')

        #Create list of paths to individual clusternpz files for this_k
        connmapClustersDirList = list()
        thiskString = str(this_k)
        #Only select the directories that have the specified seed and also end with k=this_k
        for connmapClustersDir in os.listdir(path):
            if ((connmapClustersDir.find(seed) != -1) and (connmapClustersDir.endswith(thiskString))):
                connmapClustersDirList.append(connmapClustersDir)
                newpath = path + '/' + connmapClustersDirList[0]
                individualClusters_thisk_npzList = list()
                for individualClusters_thisk_npz in os.listdir(newpath):
                    if individualClusters_thisk_npz.endswith('.npz'):
                        fullpath = newpath + '/' + individualClusters_thisk_npz
                        individualClusters_thisk_npzList.append(fullpath)
                print(individualClusters_thisk_npzList)

        #Initiate an empty dataframe that will be appended to include the values from each cluster
        df_ALLCLUSTERS_InsertRank = pd.DataFrame(columns = ['targetnames', 'ConnectivityScore', 'ConnectivityScoreNormalized', 'SortedRank'])

        #Create a list of ranks from 1(lowest connectivity) to 180(highest connectivity)
        for i,thisk_thisCluster_npz in enumerate(individualClusters_thisk_npzList):
            data_thisCluster = np.load(thisk_thisCluster_npz)
            conngroup_thisCluster = data_thisCluster['conn_group']

            connAvg_across_subs_thisCluster = conngroup_thisCluster.mean(axis=0) #Shape of input is (#sub,#voxels,#targets)
            targetAvg_across_voxels_thisCluster = connAvg_across_subs_thisCluster.mean(axis=0) #Shape of input is (#voxels,#targets)
            sum_targetAvg_across_voxels_thisCluster = np.sum(targetAvg_across_voxels_thisCluster) #Sum all the targets to get sum of conn scores
            quotients = [number / sum_targetAvg_across_voxels_thisCluster for number in targetAvg_across_voxels_thisCluster]

            df_thisCluster = pd.DataFrame({'targetnames': GlasserTable.AreaName, 'ConnectivityScore': targetAvg_across_voxels_thisCluster, 'ConnectivityScoreNormalized': quotients}, index=range(1,snakemake.params.num_targets+1))
            #Insert a column the SectionsBold column from the GlasserTable
            df_thisCluster.insert(loc=1, column='SectionsBOLD', value=GlasserTable.SectionsBOLD)
            print('\n df_thisCluster_InsertedSectionsBold k-%s and cluster-%s : \n %s' % (this_k, i+1, df_thisCluster))
            df_thisCluster_sorted_PrimarySection =  df_thisCluster.sort_values(by='SectionsBOLD', ascending=True)
            df_thisCluster_sorted_PrimarySection.index = range(1,snakemake.params.num_targets+1)
            #Insert a column that describes the sorted dataframe new index from 1-->179
            df_thisCluster_sorted_PrimarySection.insert(loc=1, column='Sorted_PrimarySection_Index', value=range(1,snakemake.params.num_targets+1))
            print('\n df_thisCluster_sorted_PrimarySection k-%s and cluster-%s : \n %s' % (this_k, i+1, df_thisCluster_sorted_PrimarySection))

            #Assign weight based on ConnectivityScore
            df_sorted_thisCluster = df_thisCluster_sorted_PrimarySection.sort_values(by='ConnectivityScoreNormalized', ascending=False)
            print('\n df_sorted_thisCluster k-%s and cluster-%s : \n %s' % (this_k, i+1, df_sorted_thisCluster))

            df_sorted_thisCluster_InsertRank = df_sorted_thisCluster
            df_sorted_thisCluster_InsertRank['SortedRank'] = pd.Series(reversed(range(1,snakemake.params.num_targets+1)),index=(df_sorted_thisCluster.index))

            df_thisCluster_sorted_PrimarySection_InsertRank = df_thisCluster_sorted_PrimarySection
            df_thisCluster_sorted_PrimarySection_InsertRank['SortedRank'] = pd.Series(reversed(range(1,snakemake.params.num_targets+1)),index=(df_sorted_thisCluster.index))
            print('\n df_sorted_thisCluster k-%s and cluster-%s with Rank : \n %s' % (this_k, i+1, df_sorted_thisCluster_InsertRank))
            print('\n df_thisCluster k-%s and cluster -%s with Rank : \n %s' % (this_k, i+1, df_thisCluster_sorted_PrimarySection_InsertRank))

            df_ALLCLUSTERS_InsertRank = df_ALLCLUSTERS_InsertRank.append(df_thisCluster_sorted_PrimarySection_InsertRank, ignore_index=True)
            print('CURRENT df_ALLCLUSTERS_InsertRank Shape = ', df_ALLCLUSTERS_InsertRank.shape)
        print('FINAL df_ALLCLUSTERS_InsertRank Shape = ', df_ALLCLUSTERS_InsertRank.shape)
        print('FINAL df_ALLCLUSTERS_InsertRank : \n', df_ALLCLUSTERS_InsertRank)

        #Add in the edge weight column
        #Multiply the SortedRank*ConnectivityScoreNormalized to get the Weight of each edge, round to integer, then add 1 to ensure that lowest edge weight value is > 0
        df_ALLCLUSTERS_InsertRank['Weight'] = df_ALLCLUSTERS_InsertRank['SortedRank'] * df_ALLCLUSTERS_InsertRank['ConnectivityScoreNormalized']
        df_ALLCLUSTERS_InsertRank['Weight'] = df_ALLCLUSTERS_InsertRank['Weight'].apply(lambda x: round(x, 0) + 1)

        Edgesdf = pd.DataFrame({'Source': EdgesSourceId, 'SLabel': EdgesSourceNames, 'Target': TargetsMultiplied.Id, 'TLabel': TargetsMultiplied.Label, 'ConnectivityScore': df_ALLCLUSTERS_InsertRank.ConnectivityScore, 'ConnectivityScoreNormalized': df_ALLCLUSTERS_InsertRank.ConnectivityScoreNormalized, 'Rank_in_Cluster': df_ALLCLUSTERS_InsertRank.SortedRank, 'Weight': df_ALLCLUSTERS_InsertRank.Weight}, index=range(0,this_k*snakemake.params.num_targets))
        print('\n Shape of this Edges is %s and Edgesdf = \n %s  ' % (Edgesdf.shape, Edgesdf))


        print('\n ------ Saving Nodesdf and Edgees df to csv ------ \n ')

        if not os.path.exists('results/gephi_input_NodesandEdges_tables_' + seed):
            os.mkdir('results/gephi_input_NodesandEdges_tables_' + seed)

        #Create savepaths
        savepathNodesdf = 'results/gephi_input_NodesandEdges_tables_' + seed + '/GephiNodes_' + seed + '_k-' + str(this_k) + '_Sorted_Overarching22Categories_WithRGBColours.csv'
        savepathEdgesdf = 'results/gephi_input_NodesandEdges_tables_' + seed + '/GephiEdges_' + seed + '_k-' + str(this_k) + '_Sorted_Overarching22Categories_NormalizedWeight.csv'

        print('savepathNodesdf = ', savepathNodesdf)
        print('savepathEdgesdf = ', savepathEdgesdf)

        #Write to Nodesdf and Edgesdf to CSV files
        Nodesdf.to_csv(savepathNodesdf, index = False, header=True)
        Edgesdf.to_csv(savepathEdgesdf, index = False, header=True)

import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from scipy import stats

targetnames = pd.read_table('results/diffparc/sub-100610/target_images_from-hcp7Tsubj_rpiriform_overlap_removed.txt', names=['TargetNames'], header=None)
targetnamesList = targetnames.TargetNames.to_list()
targetnamesList = [ line.replace('results/diffparc/sub-100610/lh_rh_targets_from-hcp7Tsubj_rpiriform_space-native_resampled_dwi_resolution_pir-overlap-removed/','') for line in targetnamesList ] #Remove that string from each line
targetnamesList = [ line.replace('_ROI.nii.gz','') for line in targetnamesList ] #Remove this string from each
cluster_range = range(2,snakemake.params.max_k+1)
for i,this_k in enumerate(cluster_range): #For 0,2 1,3 2,4
    label = 1
    thiskdirslist = list()
    thiskString = str(this_k)
    for j,thiskdir in enumerate(snakemake.input.connmap_Clusters_npz_dir):
        if thiskdir.endswith(thiskString):
            thiskdirslist.append(thiskdir)
    newClusterGroupDir = '{snakeout}_k-{thisk}'.format(snakeout=snakemake.params.connmap_Clusters_group_npz_dir_base,thisk=this_k)
    if not os.path.exists(newClusterGroupDir):
        os.mkdir(newClusterGroupDir)
    fig = plt.figure(figsize=[6.4*this_k,4.8])
    while label <= this_k:
        thislabel_namesList = list()
        labelString = str(label) + '_connmap.npz'
        print('labelString = ',labelString)
        for l,eachthiskdir in enumerate(thiskdirslist):
            thislabel_names = [f for f in os.listdir(path=eachthiskdir) if f.endswith(labelString)]
            thislabel_namesString = ''.join(thislabel_names)
            thislabel_fullpath = eachthiskdir + '/' + thislabel_namesString
            thislabel_namesList.append(thislabel_fullpath)
        thisk_thislabelonesubj = thislabel_namesList[0]
        print('thisk_thislabelonesubj = ',thisk_thislabelonesubj)
        data_thisk_thislabelonesubj = np.load(thisk_thislabelonesubj)
        affine = data_thisk_thislabelonesubj['affine']
        mask = data_thisk_thislabelonesubj['mask']
        conn_shape = data_thisk_thislabelonesubj['conn'].shape 
        nsubjects = len(thislabel_namesList)
        print('nsubjects = ',nsubjects)
        conn_group = np.zeros([nsubjects,conn_shape[0],conn_shape[1]])
        savestring = '{newClusterGroupD}/cluster{k_i:02d}_connmap_group.npz'.format(newClusterGroupD=newClusterGroupDir,k_i=label)
        for sub,npz in enumerate(thislabel_namesList):
            data_thisk_thislabel_eachsub = np.load(npz)
            conn_group[sub,:,:] = data_thisk_thislabel_eachsub['conn']
        np.savez(savestring, conn_group=conn_group,mask=mask,affine=affine)
        connmatrix_for_each_target_avggroup = conn_group.mean(axis=0) #Average across rows (subjects) to get connmap matrix average across all subjects
        print('connmatrix_for_each_target_avggroup shape = ',connmatrix_for_each_target_avggroup.shape)
        avg_conn_for_each_target_avggroup = connmatrix_for_each_target_avggroup.mean(axis=0) #Average across rows (voxels in piriform, indexed) to get connmap matrix average across all subjects
        sum_avg_conn_for_each_target_avggroup = np.sum(avg_conn_for_each_target_avggroup)
        quotients = [number / sum_avg_conn_for_each_target_avggroup for number in avg_conn_for_each_target_avggroup]
        #Create dataframe that includes columns target names and the values of the targets
        df = pd.DataFrame({'targetnames': targetnamesList, 'ConnectivityScore': avg_conn_for_each_target_avggroup, 'ConnectivityScoreNormalized': quotients}, index=range(0,snakemake.params.num_targets))
        print("Targets prior to sorting based on avg conn for k-%s and label %s = %s" % (this_k, label, df))
        dfsorted = df.sort_values(by='ConnectivityScoreNormalized', ascending=False) #Sort the target regions descending
        print("Targets sorted based on avg conn for k-%s and label %s = %s" % (this_k, label, dfsorted))
        #Plot the 10 targets with the highest connectivity to the piriform
        x = np.linspace(1,10,10) #start,stop,number of ticks
        y = np.linspace(0,1,11) #Makes it incerement by 0.1
        ax = fig.add_subplot(1,this_k,label)
        ax.set_title("k-%s Cluster-0%s" % (this_k, label))
        ax.set_xticks(x)
        ax.set_yticks(y)
        ax.set_ylim([0,0.5])
        ax.set_xticklabels(dfsorted.targetnames[0:10], rotation='vertical')
        ax.set_xlabel('Target Name')
        ax.set_ylabel('Average Connectivity to Target')
        ax.bar(x,dfsorted.ConnectivityScoreNormalized[0:10]) #Gives the target regions with the top 10 highest average connectivites from the group averaged piriform seed
        label +=1
    figuresavestring = '{newClusterGroupDir}/k-{thisk}_Individual_Clusters_Top10connectivity_to_targets.png'.format(newClusterGroupDir=newClusterGroupDir,thisk=this_k)
    print('figuresavestring is : ',figuresavestring)
    plt.savefig(figuresavestring, bbox_inches = "tight")

import sklearn
import numpy as np
import nibabel as nib

# Define a function for saving niftis 
def save_label_nii (labels,mask,affine,out_nifti):
    labels_vol = np.zeros(mask.shape)
    labels_vol[mask > 0] = labels+1 #add a 1 so label 0 is diff from bgnd
    labels_nib = nib.Nifti1Image(labels_vol,affine)
    nib.save(labels_nib,out_nifti)

data = np.load(snakemake.input.connmap_group_npz)
cluster_range = range(2,snakemake.params.max_k+1)
out_nii_list = snakemake.output

conn_group = data['conn_group']
mask = data['mask']
affine = data['affine']

# Concat subjects
conn_group_m = np.moveaxis(conn_group,0,2)
conn_concat = conn_group_m.reshape([conn_group_m.shape[0],conn_group_m.shape[1]*conn_group_m.shape[2]])

# Run spectral clustering and save output nifti
for i,k in enumerate(cluster_range):
    from sklearn.cluster import SpectralClustering
    clustering = SpectralClustering(n_clusters=k, assign_labels="discretize",random_state=0,affinity='cosine').fit(conn_concat)
    print(f'i={i}, k={k},saving {out_nii_list[i]}')
    save_label_nii(clustering.labels_,mask,affine,out_nii_list[i])