-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathGeneExpression_on_CorticalSurface.py
115 lines (90 loc) · 3.61 KB
/
GeneExpression_on_CorticalSurface.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#========================================================================
# Visualizing the expression pattern on the cortical surface
def generate_FreeSurferLUT(labels,data,mapname,filename):
"""
This function generates a colour look-up table to be used with FreeSurfer.
It needs the values to be plotted and the labels corresponding to their anatomical location.
It further accepts an input for the colour map to be used for plotting - see matplotlib documentation for accepted names.
written by Joe Bathelt, PhD
MRC Cognition & Brain Sciences Unit
"""
import matplotlib as mpl
import matplotlib.cm as cm
import pandas as pd
# Mapping the data to a colour map
norm = mpl.colors.Normalize(vmin=np.min(data), vmax=np.max(data))
cmap = eval('cm.' + mapname)
m = cm.ScalarMappable(norm=norm, cmap=cmap)
colormap = m.to_rgba(data)
colormap[:,0:3] = np.asarray(colormap[:,0:3]*255)
colormap = np.array(colormap).astype('int')
# Creating a dataframe with the labels and colour values
df = pd.DataFrame(colormap,columns=['R','G','B','A'])
df['Label'] = labels
df = df[['Label','R','G','B','A']]
# Writing the colour table
df.to_csv(filename,sep=' ')
return colormap
def generate_ROI_file(FreeSurfer_ROI_file):
"""
This script generates a dictionary of ROIs found in the FreeSurfer parcellation filename
"""
from nipype.interfaces.freesurfer import MRIConvert
mc = MRIConvert()
mc.inputs.in_file = FreeSurfer_ROI_file
mc.inputs.out_type = 'niigz'
mc.run()
import nipype.interfaces.cmtk as cmtk
rg = cmtk.ROIGen()
rg.inputs.aparc_aseg_file = FreeSurfer_ROI_file.split('.')[0] + '_out.nii.gz'
rg.inputs.use_freesurfer_LUT = True
out_file = rg.run()
return out_file
def get_parcellation_labels(label_dict):
import networkx as nx
import pandas as pd
numbers = list()
labels = list()
Labels = nx.read_gpickle(label_dict)
for label in Labels:
numbers.append(label)
labels.append(Labels[label]['labels'])
return pd.Series(labels,index=numbers,name='label')
def FreeSurfer_LUT_for_gene(gene,label_file):
"""
This script generates a colour look-up table for a expression data as generated by the French & Paus 2015 scripts.
written by Joe Bathelt, PhD
MRC Cognition & Brain Sciences Unit
"""
# Loading the expression information
import pandas as pd
import numpy as np
from scipy.stats.mstats import zscore
lh_expression = pd.read_csv(gene + '_expression_lh.csv')
lh_expression = lh_expression[gene].values
rh_expression = pd.read_csv(gene + '_expression_rh.csv')
rh_expression = rh_expression[gene].values
# Transforming expression values to z-scores
zscores = zscore(np.hstack([lh_expression,rh_expression]))
lh_expression = zscores[0:34]
rh_expression = zscores[34:]
# Removing non-cortical ROIs
import re
counter = 0
remove = list()
labels = get_parcellation_labels(generate_ROI_file(FreeSurfer_ROI_file)).values
for label in labels:
if not re.search('ctx',label):
remove.append(counter)
counter += 1
labels = np.delete(labels,remove)
lh_colormap = generate_FreeSurferLUT(labels[0:34],lh_expression,'seismic','lh-' + gene + '-ExpressionColorLUT.txt')
rh_colormap = generate_FreeSurferLUT(labels[34:],rh_expression,'seismic','rh-' + gene + '-ExpressionColorLUT.txt')
"""
## Plotting the expression values on the surface
# Generate the colour look-up table for the gene:
FreeSurfer_LUT_for_gene(gene)
# Call tksurfer and use the GUI to replace the look-up table with the one generated above:
from subprocess import call
call('tksurfer fsaverage rh inflated -annotation aparc',shell=True)
"""