OSCR

Long-term editing of brain circuits using an engineered electrical synapse.

Code ↔ Paper

13 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 13 matches · 4 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Automated FETCH output processing pipeline ↔ FETCH.py, lines 270–279 · score 0.75 · optimal bandwidth, kernel density, cross validation, fit, FETCH
  2. [2] § Methods › Quantifying MD single-unit responses to direct IL activation ↔ LinCx_Opto_Stim_mkcg.m, the whole file · a weak match · score 0.69 · Laser Power, light pulses, intervals, optical, width, stimulation
  3. [3] § Methods › Neural data acquisition and analysis in IL→MD circuit-edited mice ↔ OptoLinCx_Analysis_v2.ipynb, lines 232–298 · score 0.67 · pre stimulus dominant, outliers, voltage, window, peaks, mice
  4. [4] § Methods › Quantifying MD single-unit responses to direct IL activation ↔ LinCx_Opto_Stim.m, the whole file · a weak match · score 0.66 · Laser Power, light pulses, intervals, optical, stimulation, analog
  5. [5] § Methods › Automated FETCH output processing pipeline ↔ FETCH.py, lines 446–497 · score 0.64 · FETCH score, Q2, Q4, Q1, Q3, cells
  6. [6] § Cx34.7(M1)–Cx35(M1) potentiates a long-range circuit ↔ OptoLinCx_Analysis_v2.ipynb, lines 1499–1573 · score 0.62 · medial dorsal thalamus, infralimbic cortex, Cx34.7, MD, Cx35, IL
  7. [7] § Methods › Automated FETCH output processing pipeline ↔ FETCH.py, lines 500–554 · score 0.62 · mApple, FITC, PE, SSC, fcs, FSC
  8. [8] § Methods › Neural data acquisition and analysis in IL→MD circuit-edited mice ↔ OptoLinCx_Analysis.ipynb, lines 194–216 · score 0.60 · pre stimulus dominant, window, peaks, filtered, channel
  9. [9] § Methods › Neural data acquisition and analysis in IL→MD circuit-edited mice ↔ LinCx_Opto_Stim_mkcg.m, the whole file · a weak match · score 0.57 · baseline recording, intensities, pulse, intervals, laser, stimulated
  10. [10] § Cx34.7(M1)–Cx35(M1) potentiates a long-range circuit ↔ OptoLinCx_Analysis_v2.ipynb, lines 1354–1497 · score 0.54 · MD channel, IL channel, Cx34.7, box, width, GFP
  11. [11] § Methods › Automated FETCH output processing pipeline ↔ FETCH.py, lines 282–326 · score 0.53 · kernel density, Gaussian, bandwidth, fitted, gate, FETCH
  12. [12] § Methods › Characterizing gap junction biophysical properties using Xenopus oocytes ↔ OptoLinCx_Analysis_v2.ipynb, lines 1057–1067 · score 0.51 · Cx34.7M1, Cx35M1
  13. [13] § Methods › Neural data acquisition and analysis in IL→MD circuit-edited mice ↔ LinCx_Opto_Stim.m, the whole file · a weak match · score 0.51 · baseline recording, pulse, intervals, laser, stimulated, analog

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 841 lines · 40 KB · no license · 4 matches

  1. #e.g. $python FETCH.py -f example -p my_cool_project -s FC114_A2_A02_002.fcs FC114_C1_C01_025.fcs
  3. import argparse
  4. import matplotlib
  5. import matplotlib.pyplot as plt
  6. from matplotlib.pyplot import figure
  7. import scipy
  8. import warnings
  9. import random
  11. def _centered(arr, newsize):
  12. # Return the center newsize portion of the array.
  13. newsize = np.asarray(newsize)
  14. currsize = np.array(arr.shape)
  15. startind = (currsize - newsize) // 2
  16. endind = startind + newsize
  17. myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
  18. return arr[tuple(myslice)]
  21. scipy.signal.signaltools._centered = _centered
  23. import flowkit as fk
  24. import numpy as np
  25. from numpy.linalg import norm
  26. import scipy.stats as st
  27. from os.path import join
  28. import os
  29. #these are still libraries
  31. from sklearn.neighbors import KernelDensity
  32. from sklearn.model_selection import GridSearchCV
  33. from numpy.linalg import eig, inv
  34. from matplotlib.patches import Ellipse
  35. from matplotlib.path import Path
  36. from sklearn.model_selection import GridSearchCV
  37. from matplotlib.offsetbox import AnchoredText
  38. from matplotlib import path
  39. import pandas as pd
  40. import seaborn as sns
  42. plt.switch_backend('agg')
  43. plt.ioff()
  44. #plot formatting
  46. #def = define function; in this case... allows argument input into terminal..
  47. def parse_arguments():
  48. ap = argparse.ArgumentParser(description="Parse arguments")
  50. ap.add_argument(
  51. "-f", "--folder",
  52. default="",
  53. help="The path to a directory containing your .fcs files",
  54. type=str )
  56. ap.add_argument(
  57. "-p", "--project",
  58. default="Untitled_project",
  59. help="The project name",
  60. type=str )
  62. ap.add_argument(
  63. "-s", "--skip_renaming",
  64. default='',
  65. help="This script expects fcs files to be formatted as 'FC114_A1_A01_001.fcs'; if just a few files are not formatted like that, pass them as arguments in this function to avoid errors; otherwise edit summarize function too parse your filenames correctly", nargs='*')
  67. ap.add_argument(
  68. "-l", "--legacy_analysis",
  69. default="False",
  70. help ="An optional argument that sets a more narrow first gate and quantile cutoff in the last gate for replicability of past analysis",
  71. type=str )
  73. ap.add_argument(
  74. "-n", "--negative_control",
  75. default="None",
  76. help ="An optional argument: negative control file name w/o fcs; draws the third gate against a known negative control instead of automatic gating",
  77. type=str )
  79. ap.add_argument(
  80. "-pn", "--predefined_negative_control",
  81. default="None",
  82. help ="An optional argument: can specify x axis and y axis precise values to draw the third gate in the format e.g.: -pn 100,240",
  83. type=str )
  85. ap.add_argument(
  86. "-g", "--log_gate",
  87. default="None",
  88. help ="An optional argument: one uses log transform on the first gate for the unorthodox FETCH template",
  89. type=str )
  91. return ap.parse_args()
  93. #this function writes gatefiles that can be opened in flowjo specifically
  94. def gate_writer(vertices1, vertices2, boundaries, filename, channame1=None, channame2=None, fluorophore1=None, fluorophore2=None):
  95. gate_text = ['<?xml version="1.0" encoding="UTF-8"?>',
  96. '<gating:Gating-ML',
  97. ' xmlns:gating="http://www.isac-net.org/std/Gating-ML/v2.0/gating"',
  98. ' xmlns:data-type="http://www.isac-net.org/std/Gating-ML/v2.0/datatypes">']
  99. polygon1 = [' <gating:PolygonGate gating:id="Polygon1">',
  100. ' <gating:dimension gating:compensation-ref="uncompensated">',
  101. ' <data-type:fcs-dimension data-type:name="SSC-A" />',
  102. ' </gating:dimension>',
  103. ' <gating:dimension gating:compensation-ref="uncompensated">',
  104. ' <data-type:fcs-dimension data-type:name="FSC-A" />',
  105. ' </gating:dimension>']
  106. for vertex in vertices1:
  107. vertex = [' <gating:vertex>',
  108. ' <gating:coordinate data-type:value="' + str(vertex[0]) + '" />',
  109. ' <gating:coordinate data-type:value="' + str(vertex[1]) + '" />',
  110. ' </gating:vertex>']
  111. polygon1 = polygon1 + vertex
  112. polygon1 = polygon1 + [' </gating:PolygonGate>']
  113. if vertices2 is not None:
  114. polygon2 = [' <gating:PolygonGate gating:id="Polygon2">',
  115. ' <gating:dimension gating:compensation-ref="uncompensated">',
  116. ' <data-type:fcs-dimension data-type:name="FSC-A" />',
  117. ' </gating:dimension>',
  118. ' <gating:dimension gating:compensation-ref="uncompensated">',
  119. ' <data-type:fcs-dimension data-type:name="FSC-H" />',
  120. ' </gating:dimension>']
  121. for vertex in vertices2:
  122. vertex = [' <gating:vertex>',
  123. ' <gating:coordinate data-type:value="' + str(vertex[0]) + '" />',
  124. ' <gating:coordinate data-type:value="' + str(vertex[1]) + '" />',
  125. ' </gating:vertex>']
  126. polygon2 = polygon2 + vertex
  127. polygon2 = polygon2 + [' </gating:PolygonGate>']
  128. and_gate = [' <gating:BooleanGate gating:id="And1">',
  129. ' <data-type:custom_info>',
  130. ' Only keep results satisfying both Polygon gates',
  131. ' </data-type:custom_info>',
  132. ' <gating:and>',
  133. ' <gating:gateReference gating:ref="Polygon1" />',
  134. ' <gating:gateReference gating:ref="Polygon2" />',
  135. ' </gating:and>',
  136. ' </gating:BooleanGate>']
  137. if boundaries is not None:
  138. quadrants = [' <gating:QuadrantGate gating:id="Quadrant1" gating:parent_id="And1">',
  139. ' <gating:divider gating:id="A" gating:compensation-ref="uncompensated">',
  140. ' <data-type:fcs-dimension data-type:name="' + channame1 + '" />',
  141. ' <gating:value>' + str(boundaries[1]) + '</gating:value>',
  142. ' </gating:divider>',
  143. ' <gating:divider gating:id="B" gating:compensation-ref="uncompensated">',
  144. ' <data-type:fcs-dimension data-type:name="' + channame2 + '" />',
  145. ' <gating:value>' + str(boundaries[0]) + '</gating:value>',
  146. ' </gating:divider>',
  147. ' <gating:Quadrant gating:id="Untransfected">',
  148. ' <gating:position gating:divider_ref="A" gating:location="' + str(boundaries[1] - 1) + '" />',
  149. ' <gating:position gating:divider_ref="B" gating:location="' + str(boundaries[0] - 1) + '" />',
  150. ' </gating:Quadrant>',
  151. ' <gating:Quadrant gating:id="Double-Positive">',
  152. ' <gating:position gating:divider_ref="A" gating:location="' + str(boundaries[1] + 1) + '" />',
  153. ' <gating:position gating:divider_ref="B" gating:location="' + str(boundaries[0] + 1) + '" />',
  154. ' </gating:Quadrant>',
  155. ' <gating:Quadrant gating:id="fluorophorea">',
  156. ' <gating:position gating:divider_ref="A" gating:location="' + str(boundaries[1] + 1) + '" />',
  157. ' <gating:position gating:divider_ref="B" gating:location="' + str(boundaries[0] - 1) + '" />',
  158. ' </gating:Quadrant>',
  159. ' <gating:Quadrant gating:id="fluorophoreb">',
  160. ' <gating:position gating:divider_ref="A" gating:location="' + str(boundaries[1] - 1) + '" />',
  161. ' <gating:position gating:divider_ref="B" gating:location="' + str(boundaries[0] + 1) + '" />',
  162. ' </gating:Quadrant>',
  163. ' </gating:QuadrantGate>']
  164. second_and_gate = [' <gating:BooleanGate gating:id="And2Untransfected">',
  165. ' <data-type:custom_info>',
  166. ' Only keep results satisfying both Polygon gates and Untransfected',
  167. ' </data-type:custom_info>',
  168. ' <gating:and>',
  169. ' <gating:gateReference gating:ref="And1" />',
  170. ' <gating:gateReference gating:ref="Untransfected" />',
  171. ' </gating:and>',
  172. ' </gating:BooleanGate>']
  173. third_and_gate = [' <gating:BooleanGate gating:id="And3DoublePositive">',
  174. ' <data-type:custom_info>',
  175. ' Only keep results satisfying both Polygon gates and Double-Positive',
  176. ' </data-type:custom_info>',
  177. ' <gating:and>',
  178. ' <gating:gateReference gating:ref="And1" />',
  179. ' <gating:gateReference gating:ref="Double-Positive" />',
  180. ' </gating:and>',
  181. ' </gating:BooleanGate>']
  182. fourth_and_gate = [' <gating:BooleanGate gating:id="And4fluorophorea">',
  183. ' <data-type:custom_info>',
  184. ' Only keep results satisfying both Polygon gates and fluorophorea',
  185. ' </data-type:custom_info>',
  186. ' <gating:and>',
  187. ' <gating:gateReference gating:ref="And1" />',
  188. ' <gating:gateReference gating:ref="fluorophorea" />',
  189. ' </gating:and>',
  190. ' </gating:BooleanGate>']
  191. fifth_and_gate = [' <gating:BooleanGate gating:id="And5fluorophoreb">',
  192. ' <data-type:custom_info>',
  193. ' Only keep results satisfying both Polygon gates and fluorophoreb',
  194. ' </data-type:custom_info>',
  195. ' <gating:and>',
  196. ' <gating:gateReference gating:ref="And1" />',
  197. ' <gating:gateReference gating:ref="fluorophoreb" />',
  198. ' </gating:and>',
  199. ' </gating:BooleanGate>']
  200. gate_text = gate_text + polygon1 + polygon2 + and_gate + quadrants + \
  201. second_and_gate + third_and_gate + fourth_and_gate + fifth_and_gate + ['</gating:Gating-ML>']
  202. else:
  203. gate_text = gate_text + polygon1 + polygon2 + and_gate + ['</gating:Gating-ML>']
  204. else:
  205. gate_text = gate_text + polygon1 + ['</gating:Gating-ML>']
  206. with open(filename, "w") as f:
  207. for line in gate_text:
  208. if line not in ['\n', '\r\n']:
  209. f.write("%s\n" % line)
  211. #selection for the first gate
  212. def fitEllipse(x,y):
  213. x = x[:,np.newaxis]
  214. y = y[:,np.newaxis]
  215. D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
  216. S = np.dot(D.T,D)
  217. C = np.zeros([6,6])
  218. C[0,2] = C[2,0] = 2; C[1,1] = -1
  219. E, V = eig(np.dot(inv(S), C))
  220. n = np.argmax(np.abs(E))
  221. a = V[:,n]
  222. b,c,d,f,g,a=a[1]/2., a[2], a[3]/2., a[4]/2., a[5], a[0]
  223. num=b*b-a*c
  224. cx=(c*d-b*f)/num
  225. cy=(a*f-b*d)/num
  227. angle=0.5*np.arctan(2*b/(a-c))*180/np.pi
  228. up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
  229. down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
  230. down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
  231. a=np.sqrt(abs(up/down1))
  232. b=np.sqrt(abs(up/down2))
  233. ell=Ellipse((cx,cy),a*2.,b*2.,angle=angle)
  234. ell_coord=ell.get_verts()
  235. return [ell_coord, cx, cy, a*2, b*2, angle]
  237. #kde = kernel density estimation; matching density to color
  238. def make_kde(points):
  239. x = points[:, 0]
  240. y = points[:, 1]
  241. # Define the borders
  242. deltaX = (max(x) - min(x))/1000
  243. deltaY = (max(y) - min(y))/1000
  244. xmin = min(x) - deltaX
  245. xmax = max(x) + deltaX
  246. ymin = min(y) - deltaY
  247. ymax = max(y) + deltaY
  248. # Create meshgrid
  249. xx, yy = np.mgrid[xmin:xmax:50j, ymin:ymax:50j]
  250. positions = np.vstack([xx.ravel(), yy.ravel()])
  251. values = np.vstack([x, y])
  252. kernel = st.gaussian_kde(values)
  253. f = np.reshape(kernel(positions).T, xx.shape)
  254. fig = plt.figure(figsize=(16,16))
  255. ax = fig.gca()
  256. cset1 = ax.contour(xx, yy, f, levels=50, colors='k')
  257. figure_centre = [(xmin + xmax)/2, (ymin + ymax)/2]
  258. plt.cla()
  259. plt.clf()
  260. return [cset1.allsegs, figure_centre, x, y]
  263. #similar to above
  264. def getKernelDensityEstimation(values, x, bandwidth = 0.2, kernel = 'gaussian'):
  265. model = KernelDensity(kernel = kernel, bandwidth=bandwidth)
  266. model.fit(values[:, np.newaxis])
  267. log_density = model.score_samples(x[:, np.newaxis])
  268. return np.exp(log_density)
  270. #a helper function for kde color matching
  271. def bestBandwidth(data, minBandwidth = 0.1, maxBandwidth = 2, nb_bandwidths = 30, cv = 30):
  272. """
  273. Run a cross validation grid search to identify the optimal bandwidth for the kernel density
  274. estimation.
  275. """
  276. model = GridSearchCV(KernelDensity(),
  277. {'bandwidth': np.linspace(minBandwidth, maxBandwidth, nb_bandwidths)}, cv=cv, n_jobs=-1)
  278. model.fit(data)
  279. return model.best_params_['bandwidth']
  282. #possibly for defining the last gate
  283. def z(samplename, Z, a, vertices1, vertices2, dest, sample, fluorophore1, fluorophore2, channame1, channame2, neg_cntrl, leg_g1):
  284. fluorophores = fluorophore1 + '_' + fluorophore2
  285. filename = join(dest, 'gates.xml')
  286. d = Z[a]
  288. #Take a pseudorandom subsample of d if it is > 10000 points:
  289. # if d.shape[0] > 10000:
  290. # random.seed(42)
  291. # print(d.shape)
  292. # rand_idx = random.sample(range(d.shape[0]), k = 10000)
  293. # d = d[rand_idx, :]
  295. if len(d) == 0:
  296. warnings.warn("No sample for last gate")
  297. return [samplename, 0, None, 0, [0, 0], [0, 0, 0, 0]]
  298. new_Z = np.array([d[:, 0] + abs(min(d[:, 0])) + 1, d[:, 1] + abs(min(d[:, 1])) + 1])
  299. new_Z= new_Z.T
  300. Z_log = np.log(new_Z)
  301. try:
  302. [alls, figure_centre, x, y] = make_kde(Z_log)
  303. except ValueError:
  304. warnings.warn("Can't make kde")
  305. return [samplename, 0, None, 0, [0, 0], [0, 0, 0, 0]]
  307. max_area = 0
  308. best_top_point = None
  309. best_right_point = None
  311. xy = np.vstack([x,y])
  312. cv_bandwidth = bestBandwidth(xy.T)
  313. kde_model = KernelDensity(kernel='gaussian', bandwidth=cv_bandwidth).fit(xy.T)
  314. kde = np.exp(kde_model.score_samples(xy.T))
  315. idx = kde.argsort()
  316. candidates = []
  318. #These are variables for the negative control gating
  319. nc_top_point = 0
  320. nc_right_point = 0
  322. if neg_cntrl[0] != "None" and samplename.rsplit('.fcs')[0] != neg_cntrl[0]:
  323. best_top_point, best_right_point = neg_cntrl[1]
  324. else:
  325. for j in range(len(alls)):
  326. for ii, seg in enumerate(alls[j]):
  328. #To find the best points for negative control, identify the rightmost and topmost points on a contour
  329. if neg_cntrl[0] != "None" and samplename.rsplit('.fcs')[0] == neg_cntrl[0]:
  330. if seg[0][0] == seg[-1][0] and seg[0][1] == seg[-1][1]:
  331. top_point_log = max(seg[:,1])
  332. right_point_log = max(seg[:,0])
  333. top_point = np.exp(top_point_log) - abs(min(d[:, 1])) - 1
  334. right_point = np.exp(right_point_log) - abs(min(d[:, 0])) - 1
  335. if top_point > nc_top_point:
  336. nc_top_point = top_point
  337. if right_point > nc_right_point:
  338. nc_right_point = right_point
  339. else:
  340. #The following applies to FETCH, isn't relevant for negative control
  341. p = Path(seg) # make a polygon
  342. grid = p.contains_points(xy.T)
  343. mean_kde = np.mean(kde[grid])
  344. top_point_log = max(seg[:,1])
  345. right_point_log = max(seg[:,0])
  346. top_point = np.exp(top_point_log) - abs(min(d[:, 1])) - 1
  347. right_point = np.exp(right_point_log) - abs(min(d[:, 0])) - 1
  348. if leg_g1:
  349. transfected_cells_x = right_point >= 500
  350. transfected_cells_y = top_point >= 500
  351. else:
  352. transfected_cells_x = right_point >= 1100
  353. transfected_cells_y = top_point >= 1100
  354. plt.plot(seg[:,0], seg[:,1], '.-')
  355. if transfected_cells_x or transfected_cells_y:
  356. continue
  357. area = 0.5*np.abs(np.dot(seg[:,0],np.roll(seg[:,1],1))-np.dot(seg[:,1],np.roll(seg[:,0],1)))
  358. if area > max_area:
  359. candidates.append([area, kde[grid], top_point, right_point])
  360. if neg_cntrl[0] == "None":
  361. largest_cand = [0]
  362. if len(candidates) == 0:
  363. warnings.warn("Pipeline error on this file")
  364. return [samplename, 0, None, 0, [0, 0], [0, 0, 0, 0]]
  365. for cand in candidates:
  366. if cand[0] > largest_cand[0]:
  367. largest_cand = cand
  368. h_plt = plt.hist(largest_cand[1], bins=100)
  369. bin_count = h_plt[0]
  370. cutoff = h_plt[1]
  371. #this if/else is not currently used, but could be used to make the last gate more stringent
  372. if leg_g1:
  373. quantile_cutoff_val = 0.60
  374. else:
  375. quantile_cutoff_val = 0.30
  376. quantile_cutoff = np.quantile(cutoff, quantile_cutoff_val)
  377. for cand in candidates:
  378. if np.mean(cand[1]) < quantile_cutoff:
  379. continue
  380. if cand[0] > max_area:
  381. max_area = cand[0]
  382. best_top_point = cand[2]
  383. best_right_point = cand[3]
  385. #Plot kde lines on the log-transformed data for debugging:
  386. plt.figure(num=None, figsize=(16, 16), dpi=80, facecolor='w', edgecolor='k')
  387. for j in range(len(alls)):
  388. for ii, seg in enumerate(alls[j]):
  389. plt.plot(seg[:,0], seg[:,1], '.-')
  390. plt.scatter(Z_log[:, 0], Z_log[:, 1], s=12.5)
  391. if best_top_point != None:
  392. best_log_top = np.log(best_top_point + abs(min(d[:, 1])) + 1)
  393. best_log_right = np.log(best_right_point + abs(min(d[:, 0])) + 1)
  394. plt.plot([min(Z_log[:, 0]), max(Z_log[:, 0])], [best_log_top, best_log_top], c='black')
  395. plt.plot([best_log_right, best_log_right], [min(Z_log[:, 1]), max(Z_log[:, 1])], c='black')
  396. print(channame1)
  397. print(channame2)
  398. plt.savefig(join(dest, channame1 + '_' + channame2 + '_debug_third_gate.pdf'), format='pdf', bbox_inches='tight')
  399. plt.cla()
  400. plt.clf()
  402. boundaries = [best_right_point, best_top_point]
  403. if len(sample.channels) == 7:
  404. gname = join(dest, fluorophores + '_gates.xml')
  405. else:
  406. gname = join(dest, 'gates.xml')
  407. if neg_cntrl[0] != "None":
  408. if neg_cntrl[0] == samplename.rsplit('.fcs')[0]:
  409. best_top_point = nc_top_point
  410. best_right_point = nc_right_point
  411. boundaries = [best_right_point, best_top_point]
  412. else:
  413. best_right_point, best_top_point = neg_cntrl[1]
  414. boundaries = [best_right_point, best_top_point]
  415. samplename = samplename + "_" + channame1 + "_" + channame2
  417. gate_writer(vertices1, vertices2, boundaries, gname, channame1, channame2, fluorophore1, fluorophore2)
  418. g_strat = fk.parse_gating_xml(gname)
  419. gs_results = g_strat.gate_sample(sample)
  420. e = gs_results.get_gate_membership('And3DoublePositive')
  421. fig = figure(num=None, figsize=(16, 16), dpi=80, facecolor='w', edgecolor='k')
  423. xy = d
  424. x = d[:, 0]
  425. y = d[:, 1]
  426. x_sorted, y_sorted, kde_sorted = x[idx], y[idx], kde[idx]
  427. # parameters of the main output plot
  428. plt.scatter(x, y, c=kde, cmap = 'turbo', s=15)
  429. plt.yscale('symlog', linthresh=1000)
  430. plt.xscale('symlog', linthresh=1000)
  431. # plt.scatter(new_Z[:, 0], new_Z[:, 1], c=kde, cmap = 'turbo', s=15)
  432. # plt.yscale('log')
  433. # plt.xscale('log')
  434. #these are gate lines; min, max are the range of point values; best points define position of the gate
  435. plt.plot([min(x), max(x)], [best_top_point, best_top_point], c='black')
  436. plt.plot([best_right_point, best_right_point], [min(y), max(y)], c='black')
  437. #df is a table format for... parsed.. data..
  438. df = gs_results.report
  439. df = df.reset_index()
  440. #numbers of each individual quadrant
  441. double_positives = list(df.loc[df['gate_name'] == 'And3DoublePositive']['count'])[0]
  442. green = list(df.loc[df['gate_name'] == 'And5fluorophoreb']['count'])[0]
  443. red = list(df.loc[df['gate_name'] == 'And4fluorophorea']['count'])[0]
  444. untransfected = list(df.loc[df['gate_name'] == 'And2Untransfected']['count'])[0]
  446. #for each sample, for each pair of colors in it, export a dataframe with fluorescence values for each cell in it
  447. color_df = pd.DataFrame(columns = [channame2, channame1])
  448. color_df[channame2] = x
  449. color_df[channame1] = y
  450. color_df.to_csv(join(dest, channame1 + "_" + channame2 + ".csv"))
  452. #this is a contingency for blank samples
  453. if neg_cntrl[0] == "None" and double_positives + green + red == 0:
  454. print("Only untransfected cells found")
  455. return [samplename, 0, None, 0, boundaries, [untransfected, red, green, double_positives]]
  456. #defining the FETCH score
  457. try:
  458. FETCH_score = double_positives/(double_positives + green + red)
  459. except ZeroDivisionError:
  460. FETCH_score = 0
  462. #another contingency
  463. if neg_cntrl[0] == "None" and (FETCH_score > 0.90 or untransfected/(double_positives + green + red + untransfected) > 0.90):
  464. print("FETCH score unreasonably high -- something went wrong")
  465. return [samplename, 0, None, double_positives + green + red + untransfected, boundaries, [untransfected, red, green, double_positives]]
  466. try:
  467. r_g = red/green
  468. except ZeroDivisionError:
  469. print("Can't calculate the proportion of red to green: there is no green cells")
  470. r_g = 0
  471. ax = plt.gca()
  472. minor = matplotlib.ticker.LogLocator(base = 10.0, subs = np.arange(1.0, 10.0) * 0.1, numticks = 10)
  473. ax.yaxis.set_minor_locator(minor)
  474. ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
  475. ax.xaxis.set_minor_locator(minor)
  476. ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
  477. ax.tick_params(which='minor', length=10, width=2)
  478. ax.tick_params(which='major', length=20, width=3)
  479. plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=14)
  480. plt.setp(ax.get_yticklabels(), fontsize=14)
  481. #text boxes in the output plot
  482. txt1 = AnchoredText('Q1\n' + str(round(100*red/(double_positives + green + red + untransfected), 1)), loc="upper left", pad=0.4, borderpad=0, prop={"fontsize":14})
  483. txt2 = AnchoredText('Q2\n' + str(round(100*double_positives/(double_positives + green + red + untransfected), 1)), loc="upper right", pad=0.4, borderpad=0, prop={"fontsize":14})
  484. txt3 = AnchoredText('Q3\n' + str(round(100*green/(double_positives + green + red + untransfected), 1)), loc="lower right", pad=0.4, borderpad=0, prop={"fontsize":14})
  485. txt4 = AnchoredText('Q4\n' + str(round(100*untransfected/(double_positives + green + red + untransfected), 1)), loc="lower left", pad=0.4, borderpad=0, prop={"fontsize":14})
  486. #this puts texts boxes onto the plot
  487. ax.add_artist(txt1)
  488. ax.add_artist(txt2)
  489. ax.add_artist(txt3)
  490. ax.add_artist(txt4)
  491. ax.set_title("FETCH Score: " + str(FETCH_score))
  492. #look into the bbox, bounding box
  493. plt.savefig(join(dest, channame1 + '_' + channame2 + '_double_positive_final.pdf'), format='pdf', bbox_inches='tight')
  494. #clear axes and figure to plot next
  495. plt.cla()
  496. plt.clf()
  497. return [samplename, FETCH_score, r_g, double_positives + green + red + untransfected, boundaries, [untransfected, red, green, double_positives]]
  500. #Above, the helper files were added, Below here, the actual processing, central functions are listed
  501. def FETCH_analysis(inputlist):
  502. plt.close('all')
  503. fcs_path, samplename, dest, leg_g1, neg_cntrl, log_gate = inputlist
  504. #make a directly for an FCS file and use flowkit to parse that, to get variable called sample
  505. os.mkdir(dest)
  506. plt.grid(visible=None)
  507. sample = fk.Sample(fcs_path)
  510. fsc_a_loc = sample.channels.loc[sample.channels['pnn'].str.contains('FSC-A')].index[0]
  511. ssc_a_loc = sample.channels.loc[sample.channels['pnn'].str.contains('SSC-A')].index[0]
  512. fsc_h_loc = sample.channels.loc[sample.channels['pnn'].str.contains('FSC-H')].index[0]
  513. remaining_rows = [r for r in range(sample.channels.shape[0]) if r not in [fsc_a_loc, ssc_a_loc, fsc_h_loc]]
  514. other_chans = list(sample.channels.iloc[remaining_rows]['pnn'])
  515. other_chans = [chn for chn in other_chans if 'Time' not in chn]
  516. fluor_chan_n = len(other_chans)
  518. arr1 = sample.get_channel_events(fsc_a_loc, source='raw', subsample=False) #FSC-A
  519. arr2 = sample.get_channel_events(ssc_a_loc, source='raw', subsample=False) #SSC-A
  520. arr3 = sample.get_channel_events(fsc_h_loc, source='raw', subsample=False) #FSC-H
  522. if 'FITC' in other_chans or '1-A' in other_chans: #always have green along x axis
  523. if 'FITC' in other_chans:
  524. grn_ch_name = 'FITC'
  525. else:
  526. grn_ch_name = '1-A'
  527. x_ax_index = [chn for chn in other_chans if grn_ch_name in chn][0]
  528. x_chan_name = list(sample.channels.loc[sample.channels['pnn'].str.contains(grn_ch_name)]['pns'])[0]
  529. green_loc = sample.channels.loc[sample.channels['pnn'].str.contains(grn_ch_name)].index[0]
  530. other_chans = [chn for chn in other_chans if grn_ch_name not in chn]
  531. arr4 = sample.get_channel_events(green_loc, source='raw', subsample=False) #1-A or FITC-A(Emerald)
  533. else:
  534. x_ax_index = other_chans[0]
  535. x_chan_name = list(sample.channels.loc[sample.channels['pnn'].str.contains(other_chans[0])]['pns'])[0]
  536. first_loc = sample.channels.loc[sample.channels['pnn'].str.contains(other_chans[0])].index[0]
  537. other_chans = [chn for chn in other_chans if other_chans[0] not in chn]
  538. arr4 = sample.get_channel_events(first_loc, source='raw', subsample=False) #1-A or FITC-A(Emerald)
  540. if x_chan_name == '':
  541. x_chan_name = x_ax_index
  542. if len(other_chans) == 0: #If only got one fluorescent channel, use it as x and FSC-A as y
  543. y_chan_name = 'FSC-A'
  544. y_ax_index = 'FSC-A'
  545. Z = np.stack((arr4, arr1), axis=1)
  546. else:
  547. second_loc = sample.channels.loc[sample.channels['pnn'].str.contains(other_chans[0])].index[0]
  548. y_ax_index = other_chans[0]
  549. y_chan_name = list(sample.channels.loc[sample.channels['pnn'].str.contains(other_chans[0])]['pns'])[0]
  550. if y_chan_name == '':
  551. y_chan_name = y_ax_index
  552. arr5 = sample.get_channel_events(second_loc, source='raw', subsample=False) #5-A(RFP670), PE-Texas Red-A(mCherry), PE-A (mApple), or any other color
  554. Z = np.stack((arr4, arr5), axis=1)
  556. if fluor_chan_n == 3: # have 3 fluorescent channels
  557. third_loc = sample.channels.loc[sample.channels['pnn'].str.contains(other_chans[1])].index[0]
  558. z_chan_name = list(sample.channels.loc[sample.channels['pnn'].str.contains(other_chans[1])]['pns'])[0]
  559. z_ax_index = other_chans[1]
  560. arr6 = sample.get_channel_events(third_loc, source='raw', subsample=False)
  561. Z_ea = np.stack((arr4, arr5), axis=1)
  562. Z_er = np.stack((arr4, arr6), axis=1)
  563. Z_ar = np.stack((arr5, arr6), axis=1)
  564. elif fluor_chan_n > 3:
  565. raise Exception("Something is wrong with your channel number")
  566. #arr = array, plot.. x = 1st gate and y = 2nd gate
  568. if log_gate:
  569. alter_X = np.array([arr2 + abs(min(arr2)) + 1, arr1 + abs(min(arr1)) + 1])
  570. alter_X = alter_X.T
  571. X = np.log(alter_X)
  572. else:
  573. X = np.stack((arr2, arr1), axis=1)
  574. Y = np.stack((arr1, arr3), axis=1)
  575. #this loop goes through the contours of the first gate
  576. [alls, figure_centre, x, y] = make_kde(X)
  577. max_area = 0
  578. best_seg = None
  579. point_num = None
  580. plt.figure(num=None, figsize=(16, 16), dpi=80, facecolor='w', edgecolor='k')
  581. seg_list = []
  582. min_points = 4000
  583. for j in range(len(alls)):
  584. for ii, seg in enumerate(alls[j]):
  585. non_single_cells_x = min(seg[:,0]) <= 25000
  586. non_single_cells_y = min(seg[:,1]) <= 25000
  587. out_of_bounds_x = max(seg[:,0]) > (max(x) - 10000)
  588. out_of_bounds_y = max(seg[:,1]) > (max(y) - 10000)
  589. plt.plot(seg[:,0], seg[:,1], '.-')
  590. area = 0.5*np.abs(np.dot(seg[:,0],np.roll(seg[:,1],1))-np.dot(seg[:,1],np.roll(seg[:,0],1)))
  591. p = path.Path(seg)
  592. mask = p.contains_points(X)
  593. num_points = X[mask].shape[0]
  594. if num_points >= min_points:
  595. seg_list.append([j, area, ii])
  596. if non_single_cells_x or non_single_cells_y or out_of_bounds_x or out_of_bounds_y or len(seg[:,0])<10:
  597. continue
  598. if area > max_area:
  599. max_area = area
  600. best_seg = seg
  601. point_num = num_points
  602. if best_seg is None or point_num < min_points:
  603. if len(seg_list) == 0:
  604. print('No cells in the first gate')
  605. return [samplename, 0, None, 0]
  606. newbest = None
  607. smallest_area = np.inf
  608. for item in seg_list:
  609. if item[1] < smallest_area:
  610. smallest_area = item[1]
  611. best_seg = alls[item[0]][item[2]]
  612. #fitting an elipse to our identified best fit contour
  613. if log_gate:
  614. return_X = np.array([np.exp(X[:, 0]), np.exp(X[:, 1])])
  615. return_X = return_X.T
  616. X = np.array([return_X[:, 0] - abs(min(arr2)) - 1, return_X[:, 1] - abs(min(arr1)) - 1])
  617. X = X.T
  618. return_best_seg = np.array([np.exp(best_seg[:, 0]), np.exp(best_seg[:, 1])])
  619. return_best_seg = return_best_seg.T
  620. best_seg = np.array([return_best_seg[:, 0] - abs(min(arr2)) - 1, return_best_seg[:, 1] - abs(min(arr1)) - 1])
  621. best_seg = best_seg.T
  622. ell_coord, el_cx, el_cy, el_w, el_h, el_angle = fitEllipse(best_seg[:,0],best_seg[:,1])
  623. vertices1 = np.round(ell_coord, 0)
  624. if not leg_g1: #Keep the definition of vertices1 only if replicating old data
  625. # Step 1: Filter points whose y values are within 1000 of the target value
  626. filtered_points = X[np.abs(X[:, 1] - min(vertices1[:, 1])) <= 1000]
  627. thresh = 1000
  628. while np.shape(filtered_points)[0] == 0:
  629. thresh += 500
  630. filtered_points = X[np.abs(X[:, 1] - min(vertices1[:, 1])) <= thresh]
  632. # # Step 2: Sort the filtered points based on x values
  633. sorted_points = filtered_points[np.argsort(filtered_points[:, 0])]
  635. # # Step 3: Select the point with the smallest x value from the sorted array
  636. left_low = sorted_points[0]
  638. # # # Step 4: Filter points whose y values are at least 10 more than the left_low's
  639. # filtered_points2 = X[(X[:, 1] > (left_low[1] + 10)) & (X[:, 0] > left_low[0])]
  641. # # # Step5: get the left_mid point to get the slope of the left bound
  642. # left_mid = filtered_points2[np.argsort(filtered_points2[:, 0])][0]
  643. # # left_mid = [left_mid[1], left_mid[0]]
  644. # slope = (left_mid[1] - left_low[1]) / (left_mid[0] - left_low[0])
  645. # y_intercept = left_low[1] - slope * left_low[0]
  646. # top_left = [(max(X[:, 1]) - 10 - y_intercept) / slope, max(X[:, 1]) -10]
  647. top_left = [left_low[0], max(X[:, 1]) -10]
  648. right_top = [max(X[:, 0]) -10, max(X[:, 1]) -10]
  650. # #Step 6: get the line parallel to the ellipse's angle and perpendicular to its second principal component to define left bound
  651. rad90 = np.radians(90)
  652. if el_angle < 0:
  653. el_angle = el_angle + 90
  654. ell_side_point = [el_cx + el_h/2*np.sin(np.radians(el_angle))/np.sin(rad90), el_cy - el_h/2*np.sin(np.radians(90 - el_angle))/np.sin(rad90)]
  655. slope_right = np.tan(np.radians(el_angle))
  656. y_intercept_right = ell_side_point[1] - slope_right * ell_side_point[0]
  657. mid_right = [max(X[:, 0]) -10, slope_right*(max(X[:, 0]) -10) + y_intercept_right]
  658. low_right = [(left_low[1] - y_intercept_right)/slope_right, left_low[1]]
  659. vertices1 = np.round([left_low, top_left, right_top, mid_right, low_right], 0)
  661. #generates the .xml file
  662. filename = join(dest, 'gates.xml')
  663. gate_writer(vertices1, None, None, filename)
  664. g_strat = fk.parse_gating_xml(filename)
  665. gs_results = g_strat.gate_sample(sample)
  666. #gets the indices of selected cells to move into gate 2
  667. a = gs_results.get_gate_membership('Polygon1')
  668. b = Y[a]
  669. #len= length; a contingency
  670. if len(b) == 0:
  671. print('No cells in the second gate')
  672. return [samplename, 0, None, 0]
  673. plt.scatter(X[:, 0], X[:, 1], c=a, s=12.5)
  674. ax = plt.gca()
  675. ax.set_xlabel('SSC-A', fontsize=36)
  676. ax.set_ylabel('FSC-A', fontsize=36)
  677. ax.set_title("First Gate", fontsize=30)
  678. plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=36)
  679. plt.setp(ax.get_yticklabels(), fontsize=36)
  680. plt.savefig(join(dest, 'first_gate_KDE.pdf'), format='pdf', bbox_inches='tight')
  681. plt.cla()
  682. plt.clf()
  683. coefficients = np.polyfit(b[:, 0], b[:, 1], 1)
  684. poly = np.poly1d(coefficients)
  685. new_x = np.linspace(min(b[:, 0]), max(b[:, 0]), 2)
  686. new_y = poly(new_x)
  687. norms = []
  688. p1 = np.array([new_x[0], new_y[0]])
  689. p2 = np.array([new_x[1], new_y[1]])
  690. for point in b:
  691. d = norm(np.cross(p2-p1, p1-point))/norm(p2-p1)
  692. norms.append(d)
  693. std = np.array(norms).std()
  694. mask2 = []
  695. for i in range(len(b)):
  696. if norms[i] > 4*std:
  697. mask2.append(False)
  698. else:
  699. mask2.append(True)
  700. coefficients2 = np.polyfit(b[:, 0], b[:, 1], 1)
  701. poly2 = np.poly1d(coefficients)
  702. new_x2 = np.linspace(min(b[:, 0]), max(b[:, 0]), 2)
  703. new_y2 = poly(new_x)
  704. factor = 4*std
  705. plt.figure(num=None, figsize=(16, 16), dpi=80, facecolor='w', edgecolor='k')
  706. plt.scatter(b[:, 0], b[:, 1], c=mask2, s=12.5)
  707. plt.plot(new_x.tolist(), new_y.tolist(), marker = "o", c='red')
  708. plt.plot(new_x, [new_y[0] - factor, new_y[1] - factor], marker = "o", c='black')
  709. plt.plot(new_x, [new_y[0] + factor, new_y[1] + factor], marker = "o", c='black')
  710. ax = plt.gca()
  711. ax.set_xlabel('FSC-A', fontsize=36)
  712. ax.set_ylabel('FSC-H', fontsize=36)
  713. ax.set_title("Second Gate", fontsize=36)
  714. plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=36)
  715. plt.setp(ax.get_yticklabels(), fontsize=36)
  716. plt.savefig(join(dest, 'secondgate.pdf'), format='pdf', bbox_inches='tight')
  717. plt.cla()
  718. plt.clf()
  719. vertices2 = np.array([[new_x[0], new_y[0] - factor],
  720. [new_x[1], new_y[1] - factor],
  721. [new_x[1], new_y[1] + factor],
  722. [new_x[0], new_y[0] + factor]])
  723. if np.isnan(vertices1).any() or np.isnan(vertices2).any():
  724. return [samplename, 0, None, 0]
  725. gate_writer(vertices1, vertices2, None, filename)
  726. g_strat = fk.parse_gating_xml(filename)
  727. gs_results = g_strat.gate_sample(sample)
  728. a = gs_results.get_gate_membership('And1')
  729. c = Y[a]
  730. fluorescent_chan_names = list(sample.channels['pns'].unique())
  731. if fluor_chan_n == 1 or fluor_chan_n == 2:
  732. return z(samplename, Z, a, vertices1, vertices2, dest, sample, y_chan_name, x_chan_name, y_ax_index, x_ax_index, [neg_cntrl[0], neg_cntrl[1][0]], leg_g1)
  733. elif len(sample.channels) == 7:
  734. first = z(samplename, Z_ea, a, vertices1, vertices2, dest, sample, y_chan_name, x_chan_name, y_ax_index, x_ax_index, [neg_cntrl[0], neg_cntrl[1][0]], leg_g1)
  735. second = z(samplename, Z_er, a, vertices1, vertices2, dest, sample, z_chan_name, x_chan_name, z_ax_index, x_ax_index, [neg_cntrl[0], neg_cntrl[1][1]], leg_g1)
  736. third = z(samplename, Z_ar, a, vertices1, vertices2, dest, sample, z_chan_name, y_chan_name, z_ax_index, y_ax_index, [neg_cntrl[0], neg_cntrl[1][2]], leg_g1)
  737. return [first, second, third]
  740. def summarize(outputs, fcs_folder, project_name, skip_renaming):
  741. dataf = pd.DataFrame(outputs)
  742. print(dataf)
  743. dataf = dataf.rename({0: "File", 1: "FETCH score", 2 : "r_g", 3:"n_tot", 4: "3rd_gate_coord", 5:"raw_counts"}, axis='columns')
  744. dataf['Dubious?'] = [False for i in range(dataf.shape[0])]
  745. dataf['FETCH score'] = dataf['FETCH score'].astype(float)
  746. dataf['r_g'] = dataf['r_g'].astype(float)
  747. dataf['n_tot'] = dataf['n_tot'].astype(int)
  748. dataf['Dubious?'] = dataf.apply(lambda row : 'yes' if ((row['r_g'] >= 2) or (row['r_g'] <= 0.5) or (row['n_tot'] < 500) or np.isnan(row['r_g'])) else 'no',
  749. axis=1)
  750. dataf = dataf.drop(['r_g'], axis=1)
  751. try:
  752. dataf["File"] = dataf["File"].apply(lambda x: x.rsplit('_')[2] + '_' + x.rsplit('_')[0] + '_' + x.rsplit('_')[1] + '_' + x.rsplit('_')[3] if x not in skip_renaming else x)
  753. except IndexError:
  754. pass
  755. dataf = dataf.sort_values(by='File')
  756. dataf['Numbername'] = [i for i in range(dataf.shape[0])]
  757. sns.set(font_scale=2)
  758. figure(num=None, figsize=(32, 16), dpi=80, facecolor='w', edgecolor='k')
  759. colors = [(0, 0, 0) if dataf["Dubious?"].iloc[i] == 'no' else (1, 0, 0) for i in range(dataf["File"].unique().shape[0])]
  761. sns.set_style('ticks')
  762. sns.catplot(x='File', y='FETCH score', palette=colors, capsize=.2, kind="point", ci="sd", data=dataf, height=10, aspect=2)
  763. g = sns.swarmplot(x='File', y='FETCH score', data=dataf, color="purple", size=5)
  764. ax = plt.gca()
  765. ax.set_xlabel('FETCH_id')
  766. ax.set_ylabel('FETCH Score')
  767. plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=8)
  768. plt.setp(ax.get_yticklabels(), fontsize=26)
  769. ax.set(facecolor = "white")
  770. ax.set_title(project_name, fontsize=26)
  771. dataf['FETCH score'] = dataf['FETCH score'].fillna(0)
  772. ax.set(ylim=(0, max(dataf['FETCH score'])+0.1))
  773. plt.savefig(join(fcs_folder, project_name + ".pdf"), format='pdf', bbox_inches='tight')
  774. plt.cla()
  775. plt.clf()
  776. plt.close()
  778. dataf = dataf.set_index("File")
  779. dataf.to_csv(join(fcs_folder, project_name + ".csv"))
  782. #identify which folder contains our fcs files, etc.
  783. def main(args):
  784. fcs_folder = args.folder
  785. project_name = args.project
  786. leg_g1 = not(args.legacy_analysis == 'False') #for replicability of past analysis, add an optional argument that sets a more narrow first gate
  787. skip_renaming = args.skip_renaming
  788. if skip_renaming == '':
  789. skip_renaming = []
  790. negative_control = args.negative_control
  791. predefined_negative_control = args.predefined_negative_control
  792. log_gate = args.log_gate
  793. if log_gate != 'None':
  794. log_gate = True
  795. else:
  796. log_gate = False
  797. plt.cla()
  798. plt.clf()
  799. plt.close()
  800. plt.style.use('default')
  801. inpts = [[join(fcs_folder, samplename), samplename,
  802. join(fcs_folder, samplename.rsplit('.')[0]), leg_g1,
  803. [negative_control, [[None, None], [None, None], [None, None]]], log_gate] if
  804. (samplename != '.DS_Store' and not os.path.isdir(join(fcs_folder, samplename.rsplit('.')[0])))
  805. else None for samplename in os.listdir(fcs_folder)]
  806. inpts = list(filter(None, inpts))
  807. outstuff = []
  808. if negative_control != "None":
  809. for i_pos, el in enumerate(inpts):
  810. if el[1].rsplit('.fcs')[0] == negative_control:
  811. neg_outpt = FETCH_analysis(inpts.pop(i_pos))
  812. if len(neg_outpt) == 3: #the three fluorescent channels condition
  813. first_res = neg_outpt[0][4]
  814. second_res = neg_outpt[1][4]
  815. third_res = neg_outpt[2][4]
  816. for subel in neg_outpt:
  817. outstuff.append(subel)
  818. inpts = [[el[0], el[1], el[2], el[3], [el[4], [first_res, second_res, third_res]], log_gate] for el in inpts]
  819. else:
  820. outstuff.append(neg_outpt)
  821. inpts = [[el[0], el[1], el[2], el[3], [el[4], [neg_outpt[4], [None, None], [None, None]]], log_gate] for el in inpts]
  822. break
  823. if predefined_negative_control != "None":
  824. x_ax_thresh, y_ax_thresh = [int(val) for val in predefined_negative_control.split(',')]
  825. inpts = [[el[0], el[1], el[2], el[3], [el[4], [[x_ax_thresh, y_ax_thresh], [None, None], [None, None]]], log_gate] for el in inpts]
  827. for inp in inpts:
  828. res = FETCH_analysis(inp)
  829. if len(res) == 3:
  830. for subel in res:
  831. outstuff.append(subel)
  832. else:
  833. outstuff.append(res)
  834. #draws the aggregate plot figure and table comparing FETCH scores
  835. summarize(outstuff, fcs_folder, project_name, skip_renaming)
  836. #this is where the code actually starts; runs 'main', above
  837. if __name__ == '__main__':
  838. args = parse_arguments()
  839. main(args)

FETCH.py at commit 89e9901, no license · at the source

Overview

Authors: Elizabeth Ransey1,2, Gwenaëlle E. Thomas1,3, Elias M. Wisdom4, Agustin Almoril-Porras4, Ryan Bowman2, Elise Adamson1,5, Kathryn K. Walder-Christensen1,2, Jesse A. White1,6, Dalton N. Hughes1,3, Hannah Schwennesen2, Caly Ferguson1,2, Kay M. Tye1,6, Stephen D. Mague1,2, Longgang Niu7, Zhao-Wen Wang7, Daniel Colón-Ramos4,8, Rainbo Hultman9, Nenad Bursac5, Kafui Dzirasa1,2,3,5,10
ORCID iDs: Elias M. Wisdom, Kathryn K. Walder-Christensen, Longgang Niu, Zhao-Wen Wang, Daniel Colón-Ramos, Nenad Bursac, Kafui Dzirasa
  1. Howard Hughes Medical Institute,Chevy Chase, MD USA
  2. Deparment of Psychiatry and Behavioral Sciences, Duke University Medical Center,Durham, NC USA
  3. Department of Neurobiology, Duke University Medical Center,Durham, NC USA
  4. Department of Neuroscience and Department of Cell Biology, Program in Cellular Neuroscience, Neurodegeneration and Repair, Yale University School of Medicine,New Haven, CT USA
  5. Department of Biomedical Engineering, Duke University,Durham, NC USA
  6. Salk Institute for Biological Studies,La Jolla, CA USA
  7. Department of Neuroscience, University of Connecticut School of Medicine,Farmington, CT USA
  8. Instituto de Neurobiología, Recinto de Ciencias Médicas, Universidad de Puerto Rico,San Juan, Puerto Rico
  9. Department of Molecular Physiology and Biophysics, Department of Psychiatry, University of Iowa,Iowa City, IA USA
  10. Department of Neurosurgery, Duke University Medical Center,Durham, NC USA
Institutions: Howard Hughes Medical Institute (United States); Duke Medical Center (United States); Yale University (United States); Duke University (United States); Salk Institute for Biological Studies (United States); University of Connecticut (United States); University of Puerto Rico System (Puerto Rico); University of Iowa (United States)
Journal: Nature, volume 655, issue 8123, pages 703-715
Dates: received 9 August 2023; accepted 7 April 2026; published online 13 May 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41586-026-10501-y · PMID 42129559 · PMCID PMC13372691 · OpenAlex W7160996005
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), C. elegans (organism), cellular / molecular (subfield)
Methods: Spectral & time-frequency, Preprocessing, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Evoked potentials, fMRI & imaging, Single-unit activity, calcium imaging
Keywords: Emotion, Molecular neuroscience
MeSH: Brain*, Connexins*, Electrical Synapses*, Amino Acid Motifs, Animals, Caenorhabditis elegans, Mice, Models, Molecular, Time Factors (* major topic)
Topic: Connexins and lens biology (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: NHLBI NIH HHS (R01 HL126524, R01 HL132389, R01 HL164013, U01 HL134764); NIMH NIH HHS (R01 MH125430, R01 MH120158, DP1 MH132709, R01 MH085927); NIBIB NIH HHS (R01 EB032726); NINDS NIH HHS (R01 NS076558, DP1 NS111778); NEI NIH HHS (R21 EY029451)
Citations: cited by 3 papers (Europe PMC); 98 references in the paper

Abstract

Electrical signalling across distinct populations of brain cells underpins cognitive and emotional function. However, approaches that selectively regulate electrical signalling between two cellular components of a mammalian neural circuit remain sparse. Here we engineered an electrical synapse composed of two connexin proteins1 found in Morone americana (white perch fish)—connexin 34.7 and connexin 35—to accomplish mammalian circuit modulation. By exploiting protein mutagenesis, devising a new in vitro system for assaying connexin hemichannel docking, and performing computational modelling of hemichannel interactions, we uncovered a structural motif that contributes to electrical synapse formation. Targeting this motif, we designed connexin 34.7 and connexin 35 hemichannels that dock with each other to form an electrical synapse but not with other major connexins expressed in the mammalian central nervous system. We validated this electrical synapse in vivo using worms (Caenorhabditis elegans) and mice (Mus musculus). We demonstrate that it can strengthen communication across neural circuits composed of pairs of distinct cell types and modify behaviour accordingly. Thus, we establish ‘long-term integration of circuits using connexins’ (LinCx) for precision circuit editing in mammals.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 13 matches between paragraphs and lines of code.

carlson-lab/FETCH

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 89e99011dd920cb70f1ea87cf26c894bf9b4e3c3, 21 July 2024
Languages: Python (1)
Size: 5 files, 1 script
Software Heritage: not archived
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
2 files
At the source: github.com/carlson-lab/FETCH

carlson-lab/VMD-and-NAMD-Connexin-Protein-Simulation-Protocol

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 0c48394e07514b514bf671a383587a89017c46e8, 15 August 2025
Size: 10 files
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
2 files
At the source: github.com/carlson-lab/VMD-and-NAMD-Connexin-Protein-Simulation-Protocol

carlson-lab/OptoLinCx

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: fa62a83d2fe30cfd0e3b744174fb6cf5fcae5be4, 4 November 2024
Languages: MATLAB (65), Jupyter (2), C (1), C/C++ (1)
Size: 115 files, 69 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 2 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Signal Processing Toolbox (2 files), Matplotlib (2 files), Neo (2 files), NumPy (2 files), pandas (2 files), Pingouin (2 files), SciPy (2 files), seaborn (2 files), Statistics and Machine Learning Toolbox (1 file), scikit-learn (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
70 files
At the source: github.com/carlson-lab/OptoLinCx

Code availability

The following codes for the methods implemented in this study are available from GitHub: (1) https://github.com/carlson-lab/FETCH; (2) https://github.com/carlson-lab/VMD-and-NAMD-Connexin-Protein-Simulation-Protocol; (3) https://github.com/carlson-lab/OptoLinCx. All are available stored in an online repository (10.7924/r4r486).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 70 scripts, each with its path and the digest of its content;
  • 13 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

All data generated in support of the findings of this study are available from the corresponding author for academic purposes upon reasonable request. Such data will be made available under a material transfer agreement. Connexin gene information was procured from the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov) and the Ensembl genome browser (http://ensembl.org), and are readily accessible.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 19 authors, 2 keywords, 9 MeSH terms, 5 funders, 95 references.

Cite

This paper

Ransey, E., Thomas, G. E., Wisdom, E. M., Almoril-Porras, A., Bowman, R., Adamson, E., Walder-Christensen, K. K., White, J. A., Hughes, D. N., Schwennesen, H., Ferguson, C., Tye, K. M., Mague, S. D., Niu, L., Wang, Z.-W., Colón-Ramos, D., Hultman, R., Bursac, N., & Dzirasa, K. (2026). Long-term editing of brain circuits using an engineered electrical synapse. Nature, 655(8123), 703-715. https://doi.org/10.1038/s41586-026-10501-y

BibTeX

@article{ransey2026long,
author = {Ransey, Elizabeth and Thomas, Gwenaëlle E. and Wisdom, Elias M. and Almoril-Porras, Agustin and Bowman, Ryan and Adamson, Elise and Walder-Christensen, Kathryn K. and White, Jesse A. and Hughes, Dalton N. and Schwennesen, Hannah and Ferguson, Caly and Tye, Kay M. and Mague, Stephen D. and Niu, Longgang and Wang, Zhao-Wen and Colón-Ramos, Daniel and Hultman, Rainbo and Bursac, Nenad and Dzirasa, Kafui},
title = {{Long-term editing of brain circuits using an engineered electrical synapse}},
journal = {Nature},
year = {2026},
month = may,
volume = {655},
number = {8123},
pages = {703--715},
publisher = {Nature Portfolio},
issn = {0028-0836},
doi = {10.1038/s41586-026-10501-y},
url = {https://doi.org/10.1038/s41586-026-10501-y},
pmid = {42129559},
pmcid = {PMC13372691}
}

RIS

TY - JOUR
AU - Ransey, Elizabeth
AU - Thomas, Gwenaëlle E.
AU - Wisdom, Elias M.
AU - Almoril-Porras, Agustin
AU - Bowman, Ryan
AU - Adamson, Elise
AU - Walder-Christensen, Kathryn K.
AU - White, Jesse A.
AU - Hughes, Dalton N.
AU - Schwennesen, Hannah
AU - Ferguson, Caly
AU - Tye, Kay M.
AU - Mague, Stephen D.
AU - Niu, Longgang
AU - Wang, Zhao-Wen
AU - Colón-Ramos, Daniel
AU - Hultman, Rainbo
AU - Bursac, Nenad
AU - Dzirasa, Kafui
TI - Long-term editing of brain circuits using an engineered electrical synapse
T2 - Nature
J2 - Nature
PY - 2026
DA - 2026/05/13
VL - 655
IS - 8123
SP - 703
EP - 715
SN - 0028-0836
PB - Nature Portfolio
DO - 10.1038/s41586-026-10501-y
UR - https://doi.org/10.1038/s41586-026-10501-y
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41586-026-10501-y",
"type": "article-journal",
"title": "Long-term editing of brain circuits using an engineered electrical synapse",
"container-title": "Nature",
"author": [
{
"family": "Ransey",
"given": "Elizabeth"
},
{
"family": "Thomas",
"given": "Gwenaëlle E."
},
{
"family": "Wisdom",
"given": "Elias M."
},
{
"family": "Almoril-Porras",
"given": "Agustin"
},
{
"family": "Bowman",
"given": "Ryan"
},
{
"family": "Adamson",
"given": "Elise"
},
{
"family": "Walder-Christensen",
"given": "Kathryn K."
},
{
"family": "White",
"given": "Jesse A."
},
{
"family": "Hughes",
"given": "Dalton N."
},
{
"family": "Schwennesen",
"given": "Hannah"
},
{
"family": "Ferguson",
"given": "Caly"
},
{
"family": "Tye",
"given": "Kay M."
},
{
"family": "Mague",
"given": "Stephen D."
},
{
"family": "Niu",
"given": "Longgang"
},
{
"family": "Wang",
"given": "Zhao-Wen"
},
{
"family": "Colón-Ramos",
"given": "Daniel"
},
{
"family": "Hultman",
"given": "Rainbo"
},
{
"family": "Bursac",
"given": "Nenad"
},
{
"family": "Dzirasa",
"given": "Kafui"
}
],
"container-title-short": "Nature",
"volume": "655",
"issue": "8123",
"page": "703-715",
"DOI": "10.1038/s41586-026-10501-y",
"PMID": "42129559",
"PMCID": "PMC13372691",
"ISSN": "0028-0836",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41586-026-10501-y",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
13
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10464-w [code]
Quantifying electrostatic control of docking and binding energetics in functional Cx36 gap junctions.
Journal: Communications biology
In common: seaborn, scikit-learn, pandas, 3 other tools, cellular / molecular, 10 references
[2] doi:10.7554/elife.106496
AFD thermosensory neurons mediate tactile-dependent locomotion modulation in &lt;i&gt;C. elegans&lt;/i&gt;.
Journal: eLife
In common: C. elegans, cellular / molecular, 7 references
[3] doi:10.1016/j.celrep.2026.117590 [code]
Impaired behavioral inhibition in Fmr1 KO mice is linked to disrupted visual cortex theta oscillations.
Journal: Cell reports
In common: Neo, Pingouin, statsmodels, 7 other tools, mouse, 1 reference
[4] doi:10.1016/j.patter.2026.101590 [code]
Density-based longitudinal neuron tracking in high-density electrophysiological recordings.
Journal: Patterns (New York, N.Y.)
In common: Neo, Signal Processing Toolbox, statsmodels, 7 other tools, 1 reference
[5] doi:10.1126/sciadv.aef0343 [code]
Learning induces activation-mechanism-dependent neural plasticity in an intracortical microstimulation task.
Journal: Science advances
In common: Neo, Signal Processing Toolbox, statsmodels, 7 other tools
[6] doi:10.1016/j.isci.2026.117375 [code]
Motor priming is associated with widespread recruitment into neural ensembles and more rapid ensemble transitions.
Journal: iScience
In common: Pingouin, Signal Processing Toolbox, statsmodels, 7 other tools, 1 reference
[7] doi:10.7554/elife.108675 [code]
SynaptoTagMe, a toolkit for in vivo mapping and modulating neurotransmission at single-cell resolution.
Journal: eLife
In common: C. elegans, cellular / molecular, 3 references, author Daniel A Colón-Ramos
[8] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: Pingouin, Signal Processing Toolbox, statsmodels, 7 other tools, mouse
[9] doi:10.1038/s41467-026-73818-2 [code]
Prefrontal parvalbumin neurons mediate working memory in a task demand-dependent manner.
Journal: Nature communications
In common: statsmodels, scikit-learn, pandas, 3 other tools, mouse, 3 references
[10] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: Pingouin, Signal Processing Toolbox, statsmodels, 7 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.