Translating brain anatomy and disease from mouse to human in latent gene expression space.

1. #!/usr/bin/env python
2. import argparse
3. import os
4. import numpy as np
5. import pandas as pd
6. import nibabel as nb
8. from alleninf.api import get_probes_from_genes,\
9. get_expression_values_from_probe_ids, get_mni_coordinates_from_wells
10. from alleninf.data import get_values_at_locations, combine_expression_values
11. from alleninf.analysis import fixed_effects, approximate_random_effects,\
12. bayesian_random_effects
15. def nifti_file(string):
16. if not os.path.exists(string):
17. msg = "%r does not exist" % string
18. raise argparse.ArgumentTypeError(msg)
19. try:
20. nii = nb.load(string)
21. except IOError as e:
22. raise argparse.ArgumentTypeError(str(e))
23. except:
24. msg = "%r is not a nifti file" % string
25. raise argparse.ArgumentTypeError(msg)
26. else:
27. if len(nii.shape) == 4 and nii.shape[3] > 1:
28. msg = "%r is four dimensional" % string
29. raise argparse.ArgumentTypeError(msg)
30. return string
33. def main():
34. parser = argparse.ArgumentParser(
35. description="Compare a statistical map with gene expression patterns from Allen Human Brain Atlas.")
36. parser.add_argument(
37. "stat_map", help="Unthresholded statistical map in the form of a 3D NIFTI file (.nii or .nii.gz) in MNI space.", type=nifti_file)
38. parser.add_argument("gene_name", help="Name of the gene you want to compare your map with. For list of all available genes see: "
39. "http://help.brain-map.org/download/attachments/2818165/HBA_ISH_GeneList.pdf?version=1&modificationDate=1348783035873.",
40. type=str)
41. parser.add_argument("--inference_method", help="Which model to use: fixed - fixed effects, approximate_random - approximate random effects (default), "
42. "bayesian_random - Bayesian hierarchical model (requires PyMC3).",
43. default="approximate_random")
44. parser.add_argument("--n_samples", help="(Bayesian hierarchical model) Number of samples for MCMC model estimation (default 2000).",
45. default=2000, type=int)
46. parser.add_argument("--n_burnin", help="(Bayesian hierarchical model) How many of the first samples to discard (default 500).",
47. default=500, type=float)
48. parser.add_argument("--probes_reduction_method", help="How to combine multiple probes: average (default) or pca - use first principal component (requires scikit-learn).",
49. default="average")
50. parser.add_argument("--mask", help="Explicit mask for the analysis in the form of a 3D NIFTI file (.nii or .nii.gz) in the same space and "
51. "dimensionality as the stat_map. If not specified an implicit mask (non zero and non NaN voxels) will be used.",
52. type=nifti_file)
53. parser.add_argument("--radius", help="Radius in mm of of the sphere used to average statistical values at the location of each probe (default: 4mm).",
54. default=4, type=float)
55. parser.add_argument("--probe_exclusion_keyword", help="If the probe name includes this string the probe will not be used.",
56. type=str)
58. args = parser.parse_args()
60. print "Fetching probe ids for gene %s" % args.gene_name
61. probes_dict = get_probes_from_genes(args.gene_name)
62. print "Found %s probes: %s" % (len(probes_dict), ", ".join(probes_dict.values()))
64. if args.probe_exclusion_keyword:
65. probes_dict = {probe_id: probe_name for (probe_id, probe_name) in probes_dict.iteritems() if not args.probe_exclusion_keyword in probe_name}
66. print "Probes after applying exclusion cryterion: %s" % (", ".join(probes_dict.values()))
68. print "Fetching expression values for probes %s" % (", ".join(probes_dict.values()))
69. expression_values, well_ids, donor_names = get_expression_values_from_probe_ids(
70. probes_dict.keys())
71. print "Found data from %s wells sampled across %s donors" % (len(well_ids), len(set(donor_names)))
73. print "Combining information from selected probes"
74. combined_expression_values = combine_expression_values(
75. expression_values, method=args.probes_reduction_method)
77. print "Translating locations of the wells to MNI space"
78. mni_coordinates = get_mni_coordinates_from_wells(well_ids)
80. print "Checking values of the provided NIFTI file at well locations"
81. nifti_values = get_values_at_locations(
82. args.stat_map, mni_coordinates, mask_file=args.mask, radius=args.radius, verbose=True)
84. # preparing the data frame
85. names = ["NIFTI values", "%s expression" % args.gene_name, "donor ID"]
86. data = pd.DataFrame(np.array(
87. [nifti_values, combined_expression_values, donor_names]).T, columns=names)
88. data = data.convert_objects(convert_numeric=True)
89. len_before = len(data)
90. data.dropna(axis=0, inplace=True)
91. nans = len_before - len(data)
92. if nans > 0:
93. print "%s wells fall outside of the mask" % nans
95. if args.inference_method == "fixed":
96. print "Performing fixed effect analysis"
97. fixed_effects(data, ["NIFTI values", "%s expression" % args.gene_name])
99. if args.inference_method == "approximate_random":
100. print "Performing approximate random effect analysis"
101. approximate_random_effects(
102. data, ["NIFTI values", "%s expression" % args.gene_name], "donor ID")
104. if args.inference_method == "bayesian_random":
105. print "Fitting Bayesian hierarchical model"
106. bayesian_random_effects(
107. data, ["NIFTI values", "%s expression" % args.gene_name], "donor ID", args.n_samples, args.n_burnin)
110. if __name__ == '__main__':
111. main()

scripts.py at commit bc6c8f4, under BSD-3-Clause · at the source