OSCR

A High-Throughput Live Imaging Platform to Investigate Circuit-Dependent Regulation of Circadian Rhythms in Brain Tissue.

Code ↔ Paper

8 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 8 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › Data Analysis › Network Analysis ↔ R/processing_foo.R, lines 807–852 · score 0.80 · n_iterations, community detection, trace correlation, Leiden, isolate, igraph
  2. [2] § Methods › Data Analysis › Period Estimation ↔ R/processing_foo.R, lines 525–566 · score 0.66 · minpack.lm, nls.lm, offset, model, amplitude, phase
  3. [3] § Methods › Data Analysis › Network Analysis ↔ R/ClockCyteR.spatial-package.R, the whole file · a weak match · score 0.66 · convex hull, Leiden, deleted, igraph, segregation, weighting
  4. [4] § Methods › Statistical Analysis ↔ R/plotting_foo.R, lines 1555–1625 · score 0.60 · Shapiro Wilk normality, ANOVA, ggplot2, vectors
  5. [5] § Methods › Data Analysis › ClockCyteR ↔ R/processing_foo.R, lines 200–237 · score 0.60 · relative amplitude error, preprocess, RAE, NLLS, detrending, export
  6. [6] § Results › BMAL1 Deletion Impairs Intercluster Network Organization and Period Coherence of Axonal Ca2+ in the SCN ↔ R/processing_foo.R, lines 855–926 · score 0.53 · cluster phase variance, spatial segregation, strength, correlating, coherence, network
  7. [7] § Results › Enabling High‐Throughput Detection and Analysis of Circadian Rhythms in SCN Tissues ↔ R/processing_foo.R, lines 477–523 · score 0.53 · raw traces, FFT NLLS, smoothing, signal, amplitude
  8. [8] § Results › Spatiotemporal Characterization of Intracellular vs. Axonally‐Enriched Ca2+ Reporters in SCN Slices ↔ R/processing_foo.R, lines 1714–1791 · score 0.50 · circular variance, circadian phases, peak, Rayleigh, amplitude, traces

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

R · 1,791 lines · 56 KB · other · 6 matches

  1. # collection of functions to process data
  3. # function to align phases to a standard
  4. #' Align a vector of phases to a target Circadian Time
  5. #'
  6. #' @param data Numeric vector of phases in radians.
  7. #' @param CT Numeric; target Circadian Time in hours. Defaults to \code{6}.
  8. #' @param ... Currently unused.
  9. #'
  10. #' @return Numeric vector of aligned phases in radians.
  11. #' @keywords internal
  12. align_phase <- function(data, CT = 6, ...) {
  13. phases = data
  14. mean_phase = mean(phases, na.rm = TRUE)
  15. CT = CT
  16. radCT = (CT / 12) * pi
  17. difference = mean_phase - radCT
  18. if (radCT >= mean_phase) {
  19. phases_aligned = phases + difference
  20. } else {
  21. phases_aligned = phases - difference
  22. }
  23. return(phases_aligned)
  24. }
  27. # Function to calculate AUC using the trapz function
  28. #' Calculate the area under a fluorescence trace
  29. #'
  30. #' @param y Numeric vector of fluorescence values.
  31. #' @param normalize Logical; if \code{TRUE}, divides the AUC by the number of
  32. #' time points. Defaults to \code{TRUE}.
  33. #'
  34. #' @return Numeric; the (optionally normalised) trapezoidal AUC.
  35. #' @keywords internal
  36. calculate_auc <- function(y, normalize = TRUE) {
  37. # Trapezoidal method
  38. time <- seq(0, length.out = length(y), by = 0.5)
  39. auc_trapz <- trapz(time, y)
  40. if (normalize) {
  41. auc_trapz <- auc_trapz / length(time)
  42. }
  44. return(AUC1 = auc_trapz)
  45. }
  47. #' Compute circular statistics from a period table
  48. #'
  49. #' @param period_table Data frame of period-analysis results; must contain a
  50. #' \code{phase_rad} column.
  51. #'
  52. #' @return A named list with elements \code{mean_phase}, \code{RT} (Rayleigh
  53. #' test result), \code{vectorLength}, \code{circStats} (Rayleigh p-value),
  54. #' and \code{phase_var} (angular variance).
  55. #' @keywords internal
  56. circular_stats <- function(period_table) {
  57. period_table_clean <- period_table[complete.cases(period_table),]
  58. phases <- period_table_clean$phase_rad
  59. # Compute and plot mean direction
  60. # Calculate the mean circular statistic
  61. suppressWarnings(
  62. meanTest <- circular::mean.circular(phases,
  63. Rotation = "counter",
  64. na.rm = TRUE))
  66. # Perform a Rayleigh test on RCaMP
  67. suppressWarnings(
  68. RT <- rayleigh.test(phases))
  70. # Calculate the vector length statistic
  71. vectorLength <- RT$statistic[1]
  73. # Get the p-value from the Rayleigh test
  74. circStats <- RT$p.value[1]
  76. # Round the p-value to 3 decimal places
  77. try(circStats <- round(circStats, digits = 3),
  78. silent = TRUE,
  79. circStats <- NA
  80. )
  83. # Angular variance
  84. variance <- circular::angular.variance(phases, na.rm = TRUE)
  86. list(mean_phase = meanTest,
  87. RT = RT,
  88. vectorLength = vectorLength,
  89. circStats = circStats,
  90. phase_var = variance)
  91. }
  93. #' Composite Interaction Score (CI) between two variables
  94. #'
  95. #' @param tbl1 Data frame containing the first variable and an \code{ID}
  96. #' column.
  97. #' @param var1 Character; column name of the first variable in \code{tbl1}.
  98. #' @param tbl2 Data frame containing the second variable.
  99. #' @param var2 Character; column name of the second variable in \code{tbl2}.
  100. #' @param weight Logical; if \code{TRUE}, weights the score by
  101. #' \code{1 - |var1 - var2|}. Defaults to \code{FALSE}.
  102. #' @param merge Logical; reserved for future use. Defaults to \code{FALSE}.
  103. #'
  104. #' @return A data frame with columns \code{ID} and \code{CI_score}.
  105. #' @keywords internal
  106. CI_score <- function(tbl1, var1, tbl2, var2, weight = FALSE, merge = FALSE) {
  107. # TODO add NA vals filtering
  108. score <- sqrt(tbl1[var1] * tbl2[var2])
  109. if (weight) {
  110. weight = 1 - abs(tbl1[var1] - tbl2[var2])
  111. score <- score * data$weight
  112. }
  114. CI_results <- data.frame(ID = tbl1$ID, CI_score <- score) %>%
  115. `colnames<-`(c("ID", "CI_score"))
  117. return(CI_results)
  118. }
  120. #function to compute spatial segregation index
  121. #' Compute a spatial segregation index across clusters
  122. #'
  123. #' @param node_data Data frame of node metadata with columns \code{cluster},
  124. #' \code{X}, and \code{Y}.
  125. #' @param min_cluster_size Integer; clusters with fewer cells than this are
  126. #' excluded before computing the index. Defaults to \code{3}.
  127. #'
  128. #' @return A named list with elements \code{global_index} (numeric) and
  129. #' \code{overlap_matrix} (square matrix of pairwise cluster overlaps).
  130. #' @keywords internal
  131. compute_spatial_segregation <- function(node_data, min_cluster_size = 3) {
  133. node_data <- node_data[which(node_data$cluster != "NULL"), ]
  134. if(nrow(node_data) > 0) {
  135. node_data$cluster <- droplevels(node_data$cluster)
  136. }
  137. # Filter out small clusters if needed
  138. cluster_counts <- table(node_data$cluster)
  139. keep_clusters <- names(cluster_counts[cluster_counts >= min_cluster_size])
  140. node_data <- node_data[node_data$cluster %in% keep_clusters, ]
  142. if (length(unique(node_data$cluster)) < 2) {
  143. warning("Less than 2 clusters remaining after filtering")
  144. return(list(global_index = NA, overlap_matrix = NULL))
  145. }
  147. # Convert to sf
  148. points_sf <- st_as_sf(node_data, coords = c("X", "Y"))
  150. # Compute convex hulls
  151. hulls <- points_sf |>
  152. group_by(cluster) |>
  153. summarise(geometry = st_combine(geometry), .groups = "drop") |>
  154. mutate(geometry = st_convex_hull(geometry))
  156. n <- nrow(hulls)
  157. overlap_matrix <- matrix(
  158. NA,
  159. nrow = n,
  160. ncol = n,
  161. dimnames = list(hulls$cluster, hulls$cluster)
  162. )
  164. # Pairwise segregation
  165. for (i in 1:(n - 1)) {
  166. for (j in (i + 1):n) {
  167. inter <- st_intersection(
  168. st_geometry(hulls[i, ]),
  169. st_geometry(hulls[j, ])
  170. )
  171. if (length(inter) > 0) {
  172. area_inter <- st_area(inter)
  173. area_union <- st_area(st_union(
  174. st_geometry(hulls[i, ]),
  175. st_geometry(hulls[j, ]))
  176. )
  177. segregation <- 1 - (as.numeric(area_inter) / as.numeric(area_union))
  178. } else {
  179. segregation <- 1
  180. }
  181. overlap_matrix[i, j] <- segregation
  182. overlap_matrix[j, i] <- segregation
  183. }
  184. }
  186. # Cluster sizes for weighting
  187. cluster_sizes <- table(node_data$cluster)
  188. weights <- outer(cluster_sizes, cluster_sizes)
  189. weights <- weights[upper.tri(weights)]
  191. # Global weighted segregation
  192. seg_values <- overlap_matrix[upper.tri(overlap_matrix, diag = FALSE)]
  193. global_segregation <- weighted.mean(seg_values, weights, na.rm = TRUE)
  195. return(list(
  196. global_index = global_segregation,
  197. overlap_matrix = overlap_matrix
  198. ))
  199. }
  200. # foo to analyse TS and make circular plot
  201. #' Run FFT-NLLS period analysis on a fluorescence trace matrix
  202. #'
  203. #' @param df Numeric matrix or data frame of fluorescence traces (cells ×
  204. #' frames), or a long-format table if \code{prep_tbl = FALSE}.
  205. #' @param excludeNC Logical; if \code{TRUE}, excludes non-circadian cells
  206. #' from the returned table. Defaults to \code{FALSE}.
  207. #' @param top Numeric; upper period bound in hours. Defaults to \code{30}.
  208. #' @param bottom Numeric; lower period bound in hours. Defaults to
  209. #' \code{18}.
  210. #' @param save.trace Logical; whether to store individual fitted traces.
  211. #' Defaults to \code{FALSE}.
  212. #' @param rm.start Integer; number of initial time points to discard before
  213. #' analysis. Defaults to \code{0}.
  214. #' @param prep_tbl Logical; if \code{TRUE}, calls \code{prep_table()} to
  215. #' reshape \code{df} before analysis. Defaults to \code{TRUE}.
  216. #' @param preprocess Logical; if \code{TRUE}, smooths and detrends traces
  217. #' during table preparation. Defaults to \code{TRUE}.
  218. #' @param transp Logical; if \code{TRUE}, transposes \code{df} during table
  219. #' preparation. Defaults to \code{TRUE}.
  220. #' @param add_t Logical; if \code{TRUE}, adds a time column during table
  221. #' preparation. Defaults to \code{TRUE}.
  222. #' @param time_res Numeric; temporal resolution in hours per frame. Defaults
  223. #' to \code{0.5}.
  224. #' @param smooth Logical; if \code{TRUE}, uses the smoothed trace for
  225. #' analysis. Defaults to \code{TRUE}.
  226. #' @param invert Logical; if \code{TRUE}, inverts the fluorescence signal
  227. #' before analysis. Defaults to \code{FALSE}.
  228. #' @param filterRAE Numeric threshold for the Relative Amplitude Error;
  229. #' cells with RAE above this value are excluded. Pass \code{NULL} to
  230. #' disable filtering. Defaults to \code{NULL}.
  231. #' @param ... Currently unused.
  232. #'
  233. #' @return A named list with elements \code{period_table} (data frame of
  234. #' per-cell results), \code{period_table_unfiltered}, and optionally
  235. #' \code{traces} and \code{fit_traces}.
  236. #' @export
  237. computePeriod <- function(
  238. df,
  239. excludeNC = FALSE,
  240. top = 30,
  241. bottom = 18,
  242. save.trace = FALSE,
  243. rm.start = 0,
  244. prep_tbl = TRUE,
  245. preprocess = TRUE,
  246. transp = TRUE,
  247. add_t = TRUE,
  248. time_res = 0.5,
  249. smooth = TRUE,
  250. invert = FALSE,
  251. filterRAE = NULL,
  252. ...
  253. ) {
  255. #prepare table
  256. if (prep_tbl) {
  257. suppressWarnings(
  258. data_df <- prep_table(df, rm.start, preprocess, transp, add_t, time_res)
  259. )
  260. } else {
  261. data_df <- df
  262. }
  263. # analyse period and store result in list
  264. unique_ids = unique(data_df$ID)
  265. results_list <- list()
  266. traces_list <- list()
  267. fit_trace_list <- list()
  268. smooth_trace_list <- list()
  269. plot_list <- list()
  271. message(" |--- Computing period...")
  272. cli::cli_progress_bar("Processing", total = length(unique_ids))
  273. #pb <- txtProgressBar(min = 2, max = length(unique_ids), style = 3)
  274. for (i in seq_len(length(unique_ids))) {
  275. # Filter data for the current ID
  276. ID = unique_ids[i]
  277. toget = which(data_df$ID == ID)
  278. data_ID <- data_df[toget, ]
  280. # perform linear detrending
  281. fit <- lm(data_ID$value ~ data_ID$t)
  282. detrended_data <- residuals(fit)
  284. # if required, invert the data
  285. if (invert == TRUE) {
  286. detrended_data = detrended_data * -1
  287. }
  289. if (length(detrended_data == length(data_ID$values))) {
  290. data_ID$value <- detrended_data
  291. }
  293. # smooth data with loess
  294. smooth_param <- loess(data_ID$value ~ data_ID$t, span = 0.08)
  295. data_ID$smooth <- smooth_param$fitted
  297. # Estimate period and other parameters
  298. result <- FFT_NLLS_analyse(
  299. data_ID,
  300. minPer = bottom,
  301. maxPer = top,
  302. smooth = smooth,
  303. invert = invert
  304. )
  306. # Add the ID to the result
  307. result$ID <- as.integer(ID)
  309. # add the detrended trace to a list for export
  310. traces_list[[ID]] <- detrended_data
  311. fit_trace_list[[ID]] <- result$fitted_trace
  312. smooth_trace_list[[ID]] <- data_ID$smooth
  314. # subset results list to omit fitted trace
  315. result$fitted_trace <- NULL
  317. # Append the result to the results list
  318. results_list[[ID]] <- result
  320. # setTxtProgressBar(pb, i)
  321. # cli::cli_inform("Processing ID {i}")
  322. cli::cli_progress_update()
  323. }
  325. # close(pb)
  326. cli::cli_progress_done()
  327. message(" |--- Completed")
  328. # merge all period data and traces into tables
  330. period_table = do.call(rbind, lapply(results_list, as.data.frame))
  331. traces_table = do.call(cbind, lapply(traces_list, as.data.frame)) %>%
  332. t() %>%
  333. `rownames<-`(names(traces_list))
  334. fit_table = do.call(cbind, lapply(fit_trace_list, as.data.frame)) %>%
  335. t() %>%
  336. `rownames<-`(names(fit_trace_list))
  337. smooth_table = do.call(cbind, lapply(smooth_trace_list, as.data.frame)) %>%
  338. t() %>%
  339. `rownames<-`(names(smooth_trace_list))
  340. colnames(fit_table) = colnames(traces_table)
  341. colnames(smooth_table) = colnames(traces_table)
  344. # count total traces
  345. period_counts <- sum(!is.na(period_table$phase_h))
  346. # clean table of all values outside of limits and count valid traces
  347. period_tbl_allvals <- period_table
  348. if (excludeNC == TRUE) {
  349. topvals <- which(period_table$period >= top)
  350. bottomvals <- which(period_table$period <= bottom)
  351. exclude = c(topvals, bottomvals)
  352. period_table[exclude, ] <- NA
  353. period_tbl_allvals$keep[exclude] <- 0
  355. #print(paste("Values >", top, "h | <", bottom, "h removed", sep = ""))
  356. }
  357. if (!is.null(filterRAE)) {
  358. badvals <- which(period_table$RAE >= filterRAE)
  359. period_table[badvals, ] <- NA
  360. period_tbl_allvals$keep[badvals] <- 0
  362. #print(paste("RAE vals exceeding ", filterRAE, " removed", sep = ""))
  363. }
  364. non_na_count <- sum(!is.na(period_table$phase_h))
  366. # clean table of NA values
  367. period_table <- period_table[complete.cases(period_table),]
  369. # add a column for normalized phase values
  370. period_table$phase_norm <- normalize_phase(period_table$phase_rad)
  371. # round all digits in the table
  372. period_table <- period_table %>%
  373. dplyr::mutate(dplyr::across(.cols = -all_of("ID"), ~ round(., 2)))
  375. # list to return results and traces of period analysis
  376. return_list = list(
  377. traces = traces_table,
  378. fit_traces = fit_table,
  379. smooth_traces = smooth_table,
  380. period_table = period_table,
  381. period_table_unfiltered = period_tbl_allvals
  382. )
  384. return(return_list)
  385. }
  387. # function to generate a coherence analysis
  388. #' Compute local spatial coherence for a period-analysis variable
  389. #'
  390. #' @param period_table Data frame of period-analysis results; must contain an
  391. #' \code{ID} column and the column named by \code{variable}.
  392. #' @param grid_coord Data frame of grid coordinates with three columns: ID,
  393. #' X, Y.
  394. #' @param variable Character; name of the column in \code{period_table} to
  395. #' use as the coherence metric.
  396. #' @param radius Numeric; neighbourhood radius in micrometres. Defaults to
  397. #' \code{100}.
  398. #' @param threshold Numeric; minimum coherence value to retain a cell.
  399. #' Defaults to \code{0.1}.
  400. #' @param merge Logical; if \code{TRUE}, merges coherence results back into
  401. #' \code{period_table}. Defaults to \code{TRUE}.
  402. #' @param ... Currently unused.
  403. #'
  404. #' @return A data frame of per-cell coherence scores merged with the input
  405. #' \code{period_table}.
  406. #' @keywords internal
  407. coherence_analysis <- function(
  408. period_table,
  409. grid_coord,
  410. variable,
  411. radius = 100,
  412. threshold = 0.1,
  413. merge = TRUE,
  414. ...
  415. ) {
  416. plot_data <- grid_coord %>% `colnames<-`(c("ID", "X", "Y"))
  417. # Merge tables to include period data
  418. merged_data <- left_join(plot_data, period_table, by = "ID") %>%
  419. .[!is.na(.$keep), ]
  421. # Choose the variable to visualize (e.g., period, amplitude, error)
  422. variable <- variable
  424. # list all cells within radius value
  425. coherence_results <- lapply(seq_len(nrow(merged_data)), function(i) {
  427. cell_id <- merged_data$ID[i]
  428. x_coord <- merged_data$X[i]
  429. y_coord <- merged_data$Y[i]
  431. distances <- sqrt(
  432. (merged_data$X - x_coord)^2 +
  433. (merged_data$Y - y_coord)^2
  434. )
  436. var_value <- merged_data[[variable]][i]
  438. nearby_mask <- distances <= radius & merged_data$ID != cell_id
  439. # Get periods of nearby cells
  440. nearby_var <- merged_data[[variable]][nearby_mask]
  442. if (length(nearby_var) == 0) return(NA) # No nearby cells
  444. # Calculate the absolute differences in periods
  445. var_diffs <- abs(nearby_var - var_value)
  447. # Count how many differences are within the threshold
  448. coherence_count <- sum(var_diffs <= threshold)
  450. # calculate the ratio between the coherence count and the number of nearby cells
  451. if (length(nearby_var) > 0) {
  452. coherence_ratio <- coherence_count / length(nearby_var)
  453. } else {
  454. coherence_ratio <- NA # No nearby cells
  455. }
  456. # return a table containing the ID and coherence ratio
  457. results <- data.frame(ID = cell_id, coherence_ratio = coherence_ratio)
  459. return(results)
  460. }
  461. )
  463. coherence_results <- do.call(rbind, coherence_results)
  464. # safety check to avoid returning an incompatible output
  465. if(!identical(colnames(coherence_results), c("ID", "coherence_ratio"))) {
  466. coherence_results <- data.frame(ID = NA,
  467. coherence_ratio = NA)
  468. }
  469. if (merge) {
  470. # Merge coherence results with the original period table
  471. coherence_results <- left_join(coherence_results, period_table, by = "ID")
  472. }
  474. return(coherence_results)
  475. }
  477. # perform period analysis
  478. #' Estimate period, phase, and amplitude using FFT initialisation and NLLS fitting
  479. #'
  480. #' @param data Data frame with columns \code{t} (time in hours), \code{value}
  481. #' (raw fluorescence), and \code{smooth} (smoothed fluorescence).
  482. #' @param minPer Numeric; lower period bound in hours. Defaults to \code{16}.
  483. #' @param maxPer Numeric; upper period bound in hours. Defaults to \code{32}.
  484. #' @param smooth Logical; if \code{TRUE}, fits the smoothed trace; if
  485. #' \code{FALSE}, fits the raw trace. Defaults to \code{TRUE}.
  486. #' @param invert Logical; if \code{TRUE}, inverts the signal before fitting.
  487. #' Defaults to \code{FALSE}.
  488. #'
  489. #' @return A named list with fitted \code{period}, \code{phase},
  490. #' \code{amplitude}, \code{RAE}, and \code{fitted_trace}.
  491. #' @keywords internal
  492. FFT_NLLS_analyse <- function(
  493. data,
  494. minPer = 16,
  495. maxPer = 32,
  496. smooth = TRUE,
  497. invert = FALSE
  498. ) {
  499. if (smooth) {
  500. vals <- data$smooth
  501. } else {
  502. vals <- data$value
  503. }
  504. # Perform FFT
  505. n <- length(vals)
  506. dt <- mean(diff(data$t))
  507. fft_result <- stats::fft(vals)
  508. frequencies <- seq(0, 1 / dt, length.out = n)
  510. # Identify the dominant frequency (excluding the zero frequency)
  511. dominant_frequency <- frequencies[
  512. which.max(base::Mod(fft_result)[2:(n / 2)]) + 1
  513. ]
  514. initial_period <- 1 / dominant_frequency
  516. # Define the sinusoidal model function
  517. sinusoidal_model <- function(params, t) {
  518. amplitude <- params[1]
  519. period <- params[2]
  520. phase <- params[3]
  521. offset <- params[4]
  522. return(amplitude * sin(2 * pi * t / period + phase) + offset)
  523. }
  525. # Define the residual function for nonlinear least squares
  526. residuals <- function(params, t, y) {
  527. return(y - sinusoidal_model(params, t))
  528. }
  530. # Initial parameter estimates
  531. initial_amplitude <- (max(vals) - min(vals)) / 2
  532. initial_phase <- 0
  533. initial_offset <- mean(vals)
  534. initial_params <- c(
  535. initial_amplitude,
  536. initial_period,
  537. initial_phase,
  538. initial_offset
  539. )
  541. # Perform nonlinear least squares fitting
  542. suppressWarnings(
  543. fit <- minpack.lm::nls.lm(
  544. par = initial_params,
  545. fn = residuals,
  546. t = data$t,
  547. y = vals,
  548. lower = c(-Inf, minPer, -Inf, -Inf),
  549. upper = c(Inf, maxPer, Inf, Inf)
  550. )
  551. )
  553. # Extract fitted parameters
  554. fitted_params <- fit$par
  555. fitted_amplitude <- fitted_params[1]
  556. fitted_period <- fitted_params[2]
  557. fitted_phase <- fitted_params[3]
  558. fitted_offset <- fitted_params[4]
  560. # Ensure the amplitude is positive
  561. if (fitted_amplitude < 0) {
  562. fitted_amplitude <- abs(fitted_amplitude)
  563. fitted_params[1] <- fitted_amplitude
  564. fitted_phase <- fitted_phase + pi # Adjust the phase by 180 degrees (π radians)
  565. fitted_params[3] <- fitted_phase
  566. }
  568. # save fitted phase before manipulation
  569. fitted_phase_original <- fitted_phase
  571. if (invert) {
  572. fitted_phase <- fitted_phase #+ pi
  573. }
  575. # Wrap the phase into the 0-2pi radians range and mirror it
  576. fitted_phase <- (2 * pi - fitted_phase) %% (2 * pi)
  577. if (fitted_phase < 0) {
  578. fitted_phase <- (2 * pi - fitted_phase) + 2 * pi
  579. }
  581. # Convert phase to time units (hours)
  582. fitted_phase <- circular::as.circular(
  583. fitted_phase,
  584. type = "angles",
  585. units = "radians",
  586. rotation = "clock",
  587. template = "none",
  588. modulo = "asis",
  589. zero = 0
  590. )
  591. phase_circadian <- circular::conversion.circular(
  592. fitted_phase,
  593. units = "hours"
  594. )
  595. phase_absolute <- phase_circadian * (fitted_period / 24)
  597. # Peak-based Phase Calculation #
  599. # Find peaks using the raw data
  600. peaks <- pracma::findpeaks(vals, nups = 5, ndowns = 5)
  602. if (!is.null(peaks)) {
  603. peak_times <- data$t[peaks[, 2]] # Get peak times from peak indices
  605. # First peak
  606. first_peak_time <- min(peak_times)
  607. first_peak_phase <- (first_peak_time %% fitted_period) / fitted_period * 24 # Convert to circadian hours
  609. # Last peak
  610. last_peak_time <- max(peak_times)
  611. last_peak_phase <- (last_peak_time %% fitted_period) / fitted_period * 24
  613. # Mean phase
  614. mean_peak_phase <- mean((peak_times %% fitted_period) / fitted_period * 24)
  615. } else {
  616. first_peak_phase <- NA
  617. last_peak_phase <- NA
  618. mean_peak_phase <- NA
  619. }
  621. # Compute residual error
  622. fitted_values <- sinusoidal_model(fitted_params, data$t)
  623. residual_error <- sqrt(mean((vals - fitted_values)^2))
  625. # Calculate R-squared (GOF)
  626. ss_total <- sum((vals - mean(vals))^2)
  627. ss_res <- sum((vals - fitted_values)^2)
  628. r_squared <- 1 - (ss_res / ss_total)
  630. # Compute standard deviation of residuals
  631. residual_std <- sqrt(
  632. sum((vals - fitted_values)^2) / (length(vals) - length(fitted_params))
  633. )
  635. # Compute RAE
  636. RAE <- residual_std / fitted_amplitude
  638. # Return the results
  639. result <- list(
  640. keep = TRUE,
  641. period = fitted_period,
  642. amplitude = fitted_amplitude,
  643. FFT_phase = fitted_phase_original,
  644. phase_h = phase_absolute,
  645. phase_rad = fitted_phase,
  646. phase_circ = phase_circadian,
  647. first_peak_phase_h = first_peak_phase,
  648. last_peak_phase_h = last_peak_phase,
  649. mean_peak_phase_h = mean_peak_phase,
  650. offset = fitted_offset,
  651. error = residual_error,
  652. GOF = r_squared,
  653. RAE = RAE,
  654. fitted_trace = fitted_values
  655. )
  656. return(result)
  657. }
  659. #' Run period analysis on the mean fluorescence trace of a channel
  660. #'
  661. #' @param channel_ctx List; channel context object containing elements
  662. #' \code{grid_vals} (numeric matrix, cells × frames) and \code{invert}
  663. #' (logical).
  664. #' @param filename Character; file identifier passed to \code{computePeriod()}.
  665. #' @param time_res Numeric; temporal resolution in hours per frame.
  666. #'
  667. #' @return A one-row data frame of period-analysis results for the mean trace.
  668. #' @keywords internal
  669. mean_trace_period <- function(channel_ctx, filename, time_res) {
  671. mean_trace <- colMeans(channel_ctx$grid_vals, na.rm = TRUE)
  672. matrix_onetrace <- matrix(rep(mean_trace, each = 2), nrow = 2)
  673. colnames(matrix_onetrace) <- colnames(channel_ctx$grid_vals)
  674. rownames(matrix_onetrace) <- seq(1, nrow(matrix_onetrace))
  675. suppressMessages(
  676. period_mean_res <- computePeriod(
  677. matrix_onetrace,
  678. filename,
  679. excludeNC = TRUE,
  680. top = 36,
  681. bottom = 16,
  682. time_res = time_res,
  683. invert = channel_ctx$invert,
  684. filterRAE = 0.90)
  685. )
  686. results_mean_trace <- period_mean_res$period_table[1, ]
  688. return(results_mean_trace)
  689. }
  691. #' Compute mutual information between two normalised variables
  692. #'
  693. #' @param var1 Numeric vector of values in \code{[0, 1]}.
  694. #' @param var2 Numeric vector of values in \code{[0, 1]}, same length as
  695. #' \code{var1}.
  696. #' @param n_bins Integer; number of bins for discretisation. Defaults to
  697. #' \code{10}.
  698. #' @param normalize Logical; if \code{TRUE}, normalises MI by the joint
  699. #' entropy. Defaults to \code{FALSE}.
  700. #'
  701. #' @return Numeric; the mutual information value, or \code{NA} if any input
  702. #' values are \code{NA}.
  703. #' @keywords internal
  704. mutual_information <- function(var1, var2, n_bins = 10, normalize = FALSE) {
  706. x = var1
  707. y = var2
  708. # Check if inputs are valid
  709. if (length(x) != length(y)) {
  710. stop("x and y must be the same length")
  711. }
  712. if (any(x < 0 | x > 1, na.rm = TRUE) || any(y < 0 | y > 1, na.rm = TRUE)) {
  713. stop("Inputs must be normalized to [0,1]")
  714. }
  715. if (any(is.na(x) | is.na(y))) return(mi = NA)
  717. # Discretize
  718. x_bin <- cut(
  719. x,
  720. breaks = seq(0, 1, length.out = n_bins + 1),
  721. include.lowest = TRUE,
  722. labels = FALSE
  723. )
  724. y_bin <- cut(
  725. y,
  726. breaks = seq(0, 1, length.out = n_bins + 1),
  727. include.lowest = TRUE,
  728. labels = FALSE
  729. )
  731. # Build joint frequency table
  732. joint_table <- table(x_bin, y_bin)
  733. joint_probs <- joint_table / sum(joint_table)
  735. # Marginal probabilities
  736. px <- rowSums(joint_probs)
  737. py <- colSums(joint_probs)
  739. # Compute mutual information
  740. mi <- 0
  741. for (i in seq_along(px)) {
  742. for (j in seq_along(py)) {
  743. pxy <- joint_probs[i, j]
  744. if (pxy > 0 && px[i] > 0 && py[j] > 0) {
  745. mi <- mi + pxy * log2(pxy / (px[i] * py[j]))
  746. }
  747. }
  748. }
  750. # Optional normalization: normalized MI in [0, 1]
  751. if (normalize) {
  752. h_x <- -sum(px[px > 0] * log2(px[px > 0]))
  753. h_y <- -sum(py[py > 0] * log2(py[py > 0]))
  754. mi <- mi / min(h_x, h_y) # or (h_x + h_y)/2 for symmetric NMI
  755. }
  757. return(mi)
  758. }
  760. #' Build and cluster a cell correlation network
  761. #'
  762. #' @param filename Character; file identifier used in plot filenames.
  763. #' @param interval_name Character; interval label used in plot filenames.
  764. #' @param cell_traces Numeric matrix of fluorescence traces (cells × frames);
  765. #' row names must be cell IDs.
  766. #' @param grid_coord Data frame of grid coordinates with columns \code{ID},
  767. #' \code{X}, and \code{Y}.
  768. #' @param period_table Data frame of period-analysis results containing an
  769. #' \code{ID} column.
  770. #' @param network_dir Character; directory for saving network output files.
  771. #' @param ... Optional: \code{plot_matrix} (logical, save correlation
  772. #' heatmap); \code{filter_nonrhythmic} (logical, remove NA-weight edges).
  773. #'
  774. #' @return A named list with elements \code{igraph_obj}, \code{node_data},
  775. #' and \code{network_vars}.
  776. #'
  777. #' @export
  778. network_analysis <- function(filename,
  779. interval_name,
  780. cell_traces,
  781. grid_coord,
  782. period_table,
  783. network_dir,
  784. ...
  785. ) {
  787. args <- list(...)
  789. plot_matrix <- args$plot_matrix
  790. plot_matrix <- ifelse(is.null(plot_matrix), FALSE, plot_matrix)
  792. filter_nonrhythmic <- args$filter_nonrhythmic
  793. filter_nonrhythmic <- ifelse(is.null(filter_nonrhythmic),
  794. TRUE,
  795. filter_nonrhythmic)
  797. # remove rows containing NA values
  798. period_table <- period_table[complete.cases(period_table), ]
  800. cell_names <-data.frame("del" = 0 ,"ID" = as.integer(rownames(cell_traces)))
  801. colnames(grid_coord) <- c("ID", "X", "Y")
  802. rownames(grid_coord) <- grid_coord$ID
  804. cell_coords <- left_join(cell_names, grid_coord, by = "ID") %>% select(-del)
  805. cell_nodes <- right_join(period_table, cell_coords, by = "ID")
  807. # # Now merge with period_table
  808. cell_nodes <- cell_nodes[, c("ID", setdiff(names(cell_nodes), "ID"))]
  809. names(cell_nodes)[1] <- "name"
  810. cell_nodes$name <- as.numeric(cell_nodes$name)
  812. # Reorder cell_nodes to match cell_traces row order:
  813. # TODO add the matching as a test
  814. trace_ids <- as.integer(rownames(cell_traces))
  815. cell_nodes <- cell_nodes[match(trace_ids, cell_nodes$name), ]
  817. # correlate traces
  818. edges_corr <- traces_correlation(method = "combined",
  819. traces = cell_traces,
  820. filename = filename,
  821. interval_name = interval_name,
  822. network_dir = network_dir,
  823. grid_coord = grid_coord,
  824. cell_nodes = cell_nodes,
  825. filter_nonrhythmic = TRUE,
  826. plot_matrix = plot_matrix
  827. )
  829. # create igraph
  830. g <- igraph::graph_from_data_frame(d = edges_corr,
  831. vertices = cell_nodes,
  832. directed = FALSE)
  834. # Calculate node metrics
  835. node_metrics_tbl <- node_metrics(g)
  836. node_data <- left_join(cell_nodes, node_metrics_tbl, "name")
  838. # Local clustering coefficient for each node
  839. clust_coef <- transitivity(g, type = "localundirected", isolates = "zero")
  840. # TODO create table with cluster summary coefficients
  842. # Global clustering coefficient
  843. global_clust <- transitivity(g, type = "global")
  845. # Community detection using Louvain method (good for weighted graphs)
  846. # TODO add additional methods for community detection (i.e. Kmeans)
  847. comm <- cluster_leiden(g,
  848. weights = E(g)$weight,
  849. n_iterations = 200,
  850. resolution = 0.75,
  851. objective_function = "modularity"
  852. )
  855. # Membership vector (which cluster each node belongs to)
  856. membership_vec <- membership(comm)
  858. # Modularity score (how strong the community structure is)
  859. modularity_score <- modularity(g, membership_vec, weights = E(g)$weight)
  861. # add back to the node data
  862. node_data$cluster <- as.integer(membership_vec)
  863. node_data$clust_coef <- clust_coef
  865. # returns g, node_data and mean_phase_per_cluster
  866. ordered_clusters <- order_clusters(graph_obj = g,
  867. node_data = node_data,
  868. cluster_order = "phase",
  869. network_dir = network_dir)
  871. g <- ordered_clusters$g
  873. node_data <- ordered_clusters$node_data
  875. color_map <- assign_cluster_colors(node_data) # TODO find function
  877. # create 3D plot representing amplitude, phase coherence (stored in node_data)
  878. # and cluster size (stored in cluster_size)
  880. cluster_summary_original <- summarise_cluster_data(node_data)
  882. cluster_summary_original <- left_join(cluster_summary_original,
  883. color_map,
  884. by = c("cluster")
  885. )
  887. segregation_ind <- compute_spatial_segregation(node_data)
  889. # Filter edges based on module
  890. g_filtered <- filter_edges(g)
  892. # Prepare network layout
  893. layout_list <- network_plot_layout(g_filtered, cell_nodes, color_map)
  895. # calculate cluster statistics
  896. cluster_stats <- cluster_stats(node_data)
  898. # correlate mean cluster phase and mean cluster edges strength
  899. summary_nodes <- node_data %>%
  900. group_by(cluster) %>%
  901. summarise(
  902. size = n(),
  903. mean_phase = as.numeric(mean(phase_norm, na.rm = TRUE)),
  904. mean_strength = mean(strength, na.rm = TRUE),
  905. mean_degree = mean(degree, na.rm = TRUE),
  906. mean_clust_coef = mean(clust_coef, na.rm = TRUE)
  907. )
  909. network_vars_list = list(global_clust = global_clust,
  910. modularity_score = modularity_score,
  911. node_metrics = node_metrics_tbl,
  912. cluster_df = cluster_summary_original,
  913. summary_nodes = summary_nodes,
  914. mean_cluster_phases = cluster_stats$mean_cluster_phases,
  915. variance_cluster_phases = cluster_stats$variance_cluster_phases,
  916. cluster_phases_variance = cluster_stats$cluster_phases_variance,
  917. segregation_idx = segregation_ind
  918. )
  921. network_return_list <- list(node_data = node_data,
  922. igraph = g,
  923. network_vars = network_vars_list
  924. )
  925. return(network_return_list)
  926. }
  928. #' Normalise phases to hours centred around zero
  929. #'
  930. #' @param angles Numeric vector or circular object of phase values in radians.
  931. #'
  932. #' @return Numeric vector of phase values in hours, centred around zero in
  933. #' the range \code{[-12, 12]}.
  934. #' @keywords internal
  935. normalize_phase <- function(angles) {
  937. if(!(length(angles)>0)) return(angles_in_hours <- numeric(0))
  939. # Ensure input is a circular object
  940. angles_circ <- circular(
  941. angles,
  942. type = "angles",
  943. units = "radians",
  944. modulo = "asis"
  945. )
  947. # Compute circular mean
  948. mean_angle <- mean.circular(angles_circ, na.rm = TRUE)
  950. # Normalize angles to be centered around 0 in the [-pi, pi] range
  951. normalized_angles <- (angles - as.numeric(mean_angle) + pi) %% (2 * pi) - pi
  953. # Convert radians to hours (scale π to 12)
  954. angles_in_hours <- normalized_angles * (12 / pi)
  956. return(angles_in_hours)
  957. }
  959. # prepare data into long table format
  960. #' Reshape and optionally preprocess a trace matrix into long format
  961. #'
  962. #' @param data Numeric matrix of fluorescence traces (cells × frames), with
  963. #' an optional leading ID column.
  964. #' @param rm.start Integer; number of initial time points to remove. Defaults
  965. #' to \code{0}.
  966. #' @param preprocess Logical; if \code{TRUE}, calls \code{preprocess_data()}
  967. #' to smooth and detrend. Defaults to \code{TRUE}.
  968. #' @param transp Logical; if \code{TRUE}, transposes \code{data} before
  969. #' processing. Defaults to \code{TRUE}.
  970. #' @param add_t Logical; if \code{TRUE}, adds a time column based on
  971. #' \code{time_res}. Defaults to \code{TRUE}.
  972. #' @param time_res Numeric; temporal resolution in hours per frame.
  973. #'
  974. #' @return A long-format \code{data.table} with columns \code{ID}, \code{t},
  975. #' \code{value}, and \code{smooth}.
  976. #' @keywords internal
  977. prep_table <- function(
  978. data,
  979. rm.start = 0,
  980. preprocess = TRUE,
  981. transp = TRUE,
  982. add_t = TRUE,
  983. time_res
  984. ) {
  985. # re-transpose table
  986. if (transp) {
  987. data = t(data)
  988. }
  990. # remove initial data
  991. if (rm.start != 0) {
  992. print("removing data")
  993. data = data[-c(seq(1, rm.start)), ]
  994. }
  996. # add time column
  997. if (add_t) {
  998. t = seq(0, length.out = length(data[, 1]), by = time_res)
  999. data = cbind(t, data)
  1000. }
  1002. # perform outliers removal, smoothening and detrending
  1003. if (preprocess) {
  1004. data_clean <- preprocess_data(data, grade = 3)
  1005. } else {
  1006. data_clean <- as.data.frame(data)
  1007. }
  1009. #make it into long table format
  1010. data_long <- tidyr::pivot_longer(
  1011. data_clean,
  1012. cols = colnames(data_clean)[-1],
  1013. names_to = "ID"
  1014. ) %>%
  1015. data.table::as.data.table(key = "ID")
  1017. return(data_long)
  1018. }
  1020. # function to remove outliers, smooth and detrend
  1021. #' Remove outliers, smooth, and detrend fluorescence traces
  1022. #'
  1023. #' @param data Numeric matrix of fluorescence traces with a leading time
  1024. #' column (time × cells).
  1025. #' @param grade Integer; polynomial degree used for detrending.
  1026. #' @param mode Character; smoothing method — currently only
  1027. #' \code{"mov_avg"} (moving average) is supported. Defaults to
  1028. #' \code{"mov_avg"}.
  1029. #' @param parallel Logical; reserved for future parallelisation. Defaults to
  1030. #' \code{FALSE}.
  1031. #'
  1032. #' @return A numeric matrix of the same dimensions as \code{data} with
  1033. #' smoothed and detrended traces.
  1034. #' @export
  1035. preprocess_data <- function(data, grade, mode = "mov_avg", parallel = FALSE) {
  1036. # load table
  1037. data_to_clear <- data
  1039. # create table to store smoothened data
  1040. n <- nrow(data)
  1041. p <- ncol(data)
  1043. clear_mat <- matrix(NA_real_, n, p)
  1044. clear_mat[,1] <- data[,1]
  1045. # get time values
  1046. x_vals <- data[, 1]
  1047. t_index <- seq_len(nrow(data))
  1048. #message(" |--- Preprocessing data...")
  1049. # pb <- txtProgressBar(
  1050. # min = 2,
  1051. # max = length(colnames(data_to_clear)),
  1052. # style = 3
  1053. # )
  1054. pad <- 2
  1056. # for cycle to go through each column and smooth the data
  1057. # TODO implement parallelisation
  1058. for (i in 2:p) {
  1060. timeserie <- data_to_clear[, i]
  1062. # create vector for cleaned data
  1063. cleaned_ts <- timeserie
  1065. switch (mode,
  1066. "mov_avg" = {
  1068. trend <- stats::lm(timeserie ~ poly(t_index, grade))
  1069. detrended <- resid(trend)
  1071. # outliers detection and removal
  1072. sigma <- mad(detrended, constant = 1)
  1073. outliers <- abs(detrended) > 2 * sigma
  1075. # delete outliers from data
  1076. detrended[outliers] <- NA
  1077. # interpolate missing values
  1078. interpolated_ts <- imputeTS::na_interpolation(detrended, option = "spline")
  1080. # # perform smoothening of the data
  1081. # fit <- smooth.spline(t_index, interpolated_ts, spar = 0.6)
  1082. # smoothed_values <- predict(fit, t_index)$y
  1084. # add edges padding
  1085. x_pad <- c(
  1086. interpolated_ts[1:pad],
  1087. interpolated_ts,
  1088. interpolated_ts[(length(interpolated_ts) - pad + 1):length(interpolated_ts)]
  1089. )
  1090. # moving average method
  1091. averaged_pad <- stats::filter(
  1092. x_pad,
  1093. rep(1 / 5, 5),
  1094. sides = 2
  1095. )
  1096. # remove edges
  1097. averaged_ts <- averaged_pad[(pad + 1):(length(averaged_pad) - pad)]
  1099. # add smoothened data to matrix
  1100. clear_mat[, i] <- averaged_ts
  1102. },
  1103. "loess" = {
  1107. # perform first loess approximation
  1108. smoothed_series <- loess(timeserie ~ t_index, span = 0.08)
  1109. smoothed_values <- predict(
  1110. smoothed_series,
  1111. newdata = data.frame(x = x_vals)
  1112. )
  1114. # calculate outliers based on distance from loess curve and sd
  1115. stdev_ts <- sd(smoothed_values)
  1116. outliers <- which(abs(smoothed_values - timeserie) > stdev_ts * 0.6)
  1118. # delete outliers from data
  1119. cleaned_ts[outliers] <- NA
  1121. # perform second loess approximation on clean data
  1122. smoothed_clean_series <- loess(
  1123. cleaned_ts ~ seq_along(cleaned_ts),
  1124. span = 0.08
  1125. )
  1126. smoothed_clean_values <- predict(
  1127. smoothed_clean_series,
  1128. newdata = data.frame(x = x_vals)
  1129. )
  1131. # interpolate missing values before detrending
  1132. smoothed_clean_values <- imputeTS::na_interpolation(
  1133. smoothed_clean_values,
  1134. option = "spline"
  1135. )
  1137. # add detrending step
  1138. smoothed_clean_detr_values <- astsa::detrend(smoothed_clean_values, grade)
  1140. # add smoothened data to matrix
  1141. clear_mat[, i] <- smoothed_clean_detr_values
  1143. # update progressbar
  1144. #setTxtProgressBar(pb, i)
  1145. }
  1146. )
  1147. }
  1149. #close(pb)
  1150. #message(" |--- Completed")
  1151. clear_data <- as.data.frame(clear_mat)
  1152. # re add colnames
  1153. colnames(clear_data) <- colnames(data)
  1154. # return new table
  1155. return(clear_data)
  1156. }
  1159. # function to align traces to the mean phase of the group
  1160. #' Circularly shift traces so their phases align to the group mean
  1161. #'
  1162. #' @param traces_table Numeric matrix of fluorescence traces (cells × frames).
  1163. #' @param period_table Data frame of period-analysis results with a
  1164. #' \code{phase_rad} column.
  1165. #' @param remove_start Integer; frames to discard from the start after
  1166. #' shifting. Defaults to \code{0}.
  1167. #' @param remove_end Integer; frames to discard from the end after shifting.
  1168. #' Defaults to \code{0}.
  1169. #' @param align_to Numeric; target phase in hours; if \code{NA}, uses the
  1170. #' circular mean of the group. Defaults to \code{NA}.
  1171. #' @param debug Logical; if \code{TRUE}, enables debugging output. Defaults
  1172. #' to \code{FALSE}.
  1173. #'
  1174. #' @return A numeric matrix of phase-aligned traces, same dimensions as
  1175. #' \code{traces_table} (minus any removed frames).
  1176. #' @keywords internal
  1177. phase_align_trace <- function(
  1178. traces_table,
  1179. period_table,
  1180. remove_start = 0,
  1181. remove_end = 0,
  1182. align_to = NA,
  1183. debug = FALSE
  1184. ) {
  1185. # TODO check align phase value
  1187. # transform rad values in plus minus Pi
  1188. phases = circular::minusPiPlusPi(circular(
  1189. period_tbl$phase_rad,
  1190. units = "rad"
  1191. ))
  1192. if (is.na(align_to)) {
  1193. align_phase = circular::mean.circular(phases, na.rm = TRUE)
  1194. } else {
  1195. align_phase = circular::as.circular(
  1196. (align_to / 12) * pi,
  1197. type = "angles",
  1198. units = "radians",
  1199. rotation = "clock",
  1200. template = "none",
  1201. modulo = "asis",
  1202. zero = 0
  1203. )
  1204. }
  1205. ph_diff = align_phase - phases
  1206. adjust = circular::minusPiPlusPi(circular(ph_diff, units = "rad"))
  1207. adjust_frames = round(
  1208. circular::conversion.circular(adjust, units = "hours") * 2,
  1209. 0
  1210. )
  1211. # initialize new table to store vecs
  1212. traces_aligned <- replace(traces_table, all(), NA)
  1213. for (i in 1:dim(traces_table)[1]) {
  1214. # get vector
  1215. trace = traces_table[i, ]
  1216. adjust_fct = adjust_frames[i]
  1217. if (is.na(adjust_fct)) {
  1218. clean_trace = rep(NA, dim(traces_table)[2])
  1219. } else if (adjust_fct == 0) {
  1220. clean_trace = trace
  1221. } else if (adjust_fct < 0) {
  1222. clean_trace = c(
  1223. rep_len(NA, abs(adjust_fct)),
  1224. trace[1:(length(trace) - abs(adjust_fct))]
  1225. )
  1226. } else if (adjust_fct > 0) {
  1227. clean_trace = c(
  1228. trace[(abs(adjust_fct) + 1):length(trace)],
  1229. rep_len(NA, abs(adjust_fct))
  1230. )
  1231. }
  1232. traces_aligned[i, ] = clean_trace
  1233. if (debug) {
  1234. # browser()
  1235. plot(x = 0:(dim(traces_table)[2] - 1), y = trace, col = "red", type = "l")
  1236. lines(
  1237. x = 0:(dim(traces_table)[2] - 1),
  1238. y = clean_trace,
  1239. col = "blue",
  1240. type = "l"
  1241. )
  1242. }
  1243. }
  1244. # trim start and end of table
  1245. # evaluate if all have been shifted in one direction
  1246. # min_shift = min(adjust_frames, na.rm = TRUE)
  1247. # max_shift = max(adjust_frames, na.rm = TRUE)
  1248. # shift = min_shift*max_shift
  1249. # if(shift > 0){
  1250. # # case when all the shift happens in the same direction (need to cut only from one of the two ends)
  1251. # # TODO complete trimming in another function
  1252. # } else if(shift < 0){
  1253. # # case when shift happens in both direction (you can use the min-max rune simply)
  1254. # traces_aligned_trim = traces_aligned[ , (max(adjust_frames, na.rm = TRUE)+1):(dim(traces_aligned)[2]+min(adjust_frames, na.rm = TRUE))]
  1255. # }
  1257. #traces_aligned_trim = traces_aligned[ , (max(adjust_frames, na.rm = TRUE)+1):(dim(traces_aligned)[2]+min(adjust_frames, na.rm = TRUE))]
  1259. return(traces_aligned)
  1260. }
  1261. ranges_copy <- function(params) {
  1262. # TODO remove if and make function proper
  1263. if(!is.null(params$plotting$ranges)){
  1264. message("\nLoading plot ranges from previous experiment")
  1265. old_params_path <- file.path(previous_exp, "summary_stats", "plot_limits.rds")
  1266. old_params <- readRDS(old_params_path)
  1268. params$plotting$ranges <- old_params$plotting$ranges
  1269. '
  1270. # assign all period ranges from previous experiment into params
  1271. if(Channel2_an){
  1272. period_range_ch2 <- plot_ranges$period_range_ch2
  1273. amplitude_range_ch2 <- plot_ranges$amplitude_range_ch2
  1274. RAE_range_ch2 <- plot_ranges$RAE_range_ch2
  1275. AUC_range_ch2 <- plot_ranges$AUC_range_ch2
  1276. phase_y_lims_Ch2 <- plot_ranges$phase_y_lims_Ch2
  1277. }
  1278. if(Channel1_an & Channel2_an){
  1279. phase_y_lims_Chx <- plot_ranges$phase_y_lims_Chx
  1280. }
  1281. '
  1282. Y_range_loaded = TRUE
  1283. }
  1284. }
  1286. #' Compute per-variable display ranges across all files and intervals
  1287. #'
  1288. #' @param params List of analysis parameters produced by \code{make_params()};
  1289. #' must contain \code{channels} and \code{time$intervals}.
  1290. #' @param file_rows Data frame of file metadata; must contain columns
  1291. #' \code{file_id}, \code{file_path}, and \code{folder_path}.
  1292. #'
  1293. #' @return A named list (by channel) of per-variable min/max ranges used
  1294. #' to set consistent colour scales across plots.
  1295. #' @export
  1296. ranges_calculation <- function(params, file_rows) {
  1298. channel_ranges <- purrr::map(names(params$channels), function(ch_id) {
  1299. if (!params$channels[[ch_id]]$enabled) return(NULL)
  1300. # get all period tables
  1301. # TODO create a save function for parameters and retrieve params from path
  1303. # TODO make these two into functions
  1304. period_tbl_list <- unlist(
  1305. purrr::pmap(
  1306. file_rows,
  1307. function(file_id, file_path, folder_path,...) {
  1309. file_row <- list(file_id = file_id,
  1310. file_path = file_path,
  1311. folder_path = folder_path)
  1313. # Create a list for each interval
  1314. interval_tables <- setNames(
  1315. lapply(names(params$time$intervals), function(int_name) {
  1316. period_tbl_path <- file.path(file_row$folder_path, "rds",
  1317. paste0(file_row$file_id, "_",
  1318. int_name, "_",
  1319. ch_id, "_period_tbl_clean.rds"))
  1321. if(file.exists(period_tbl_path)) {
  1322. readRDS(period_tbl_path)
  1323. } else {
  1324. warning("File not found: ", period_tbl_path)
  1325. NULL
  1326. }
  1327. }),
  1328. paste0(file_id, "_", names(params$time$intervals))
  1329. )
  1331. return(interval_tables)
  1332. }),
  1333. recursive = FALSE
  1334. )
  1336. # Apply the same pattern to auc_tbl_list:
  1337. auc_tbl_list <- unlist(
  1338. purrr::pmap(
  1339. file_rows,
  1340. function(file_id, file_path, folder_path,...) {
  1342. file_row <- list(file_id = file_id,
  1343. file_path = file_path,
  1344. folder_path = folder_path)
  1346. interval_tables <- setNames(
  1347. lapply(names(params$time$intervals), function(int_name) {
  1348. auc_tbl_path <- file.path(file_row$folder_path, "rds",
  1349. paste0(file_row$file_id, "_",
  1350. int_name, "_",
  1351. ch_id, "_auc_results.rds"))
  1353. if(file.exists(auc_tbl_path)) {
  1354. readRDS(auc_tbl_path)
  1355. } else {
  1356. warning("File not found: ", auc_tbl_path)
  1357. NULL
  1358. }
  1359. }),
  1360. paste0(file_id, "_", names(params$time$intervals))
  1361. )
  1363. return(interval_tables)
  1364. }),
  1365. recursive = FALSE
  1366. )
  1367. # period_tbl_list <- setNames( # TODO make these two into functions
  1368. # purrr::pmap(
  1369. # file_rows,
  1370. # function(file_id, file_path, folder_path, ...) {
  1371. #
  1372. # file_row <- list(file_id = file_id,
  1373. # file_path = file_path,
  1374. # folder_path = folder_path)
  1375. #
  1376. # # generate period_tbl path
  1377. # # TODO fix case with multiple intervals, where it should gather period
  1378. # # tables from all intervals analyzed
  1379. # browser()
  1380. # period_tbl_path <- file.path(file_row$folder_path, "rds",
  1381. # paste0(file_row$file_id, "_",
  1382. # names(params$time$intervals), "_",
  1383. # ch_id, "_period_tbl_clean.rds"
  1384. # ))
  1385. # period_tbl <- readRDS(period_tbl_path)
  1386. #
  1387. # return(period_tbl)
  1388. # }),
  1389. # file_rows$file_id)
  1390. #
  1391. # auc_tbl_list <- setNames(
  1392. # purrr::pmap(
  1393. # file_rows,
  1394. # function(file_id, file_path, folder_path, ...) {
  1395. #
  1396. # file_row <- list(file_id = file_id,
  1397. # file_path = file_path,
  1398. # folder_path = folder_path)
  1399. #
  1400. # # generate auc_tbl path
  1401. # # TODO fix case with multiple intervals, where it should gather period
  1402. # # tables from all intervals analyzed
  1403. # auc_tbl_path <- file.path(file_row$folder_path, "rds",
  1404. # paste0(file_row$file_id, "_",
  1405. # names(params$time$intervals), "_",
  1406. # ch_id, "_auc_results.rds"
  1407. # ))
  1408. # auc_tbl <- readRDS(auc_tbl_path)
  1409. #
  1410. # return(auc_tbl)
  1411. # }),
  1412. # file_rows$file_id)
  1413. '
  1414. period_tbl_allpaths_Ch1 = readRDS(
  1415. file.path(wd, "summary_stats",
  1416. "period_tbl_allpaths_Ch1.rds"))
  1417. Y_range_loaded = FALSE
  1419. # load all period tables into a list
  1421. period_tbl_list_Ch1 <- setNames(lapply(seq_along(period_tbl_allpaths_Ch1),
  1422. function(i) {
  1423. readRDS(period_tbl_allpaths_Ch1[i])
  1424. }),
  1425. filenames)
  1426. '
  1428. # In each object of the list period_tbl_list, calculate the min, max
  1429. # and sd for the required variables
  1431. period_ranges <- lapply(period_tbl_list, function(tbl) {
  1432. if(!nrow(tbl) > 0) return(NULL)
  1433. data.frame(
  1434. period_min = min(tbl$period, na.rm = TRUE),
  1435. period_max = max(tbl$period, na.rm = TRUE),
  1436. period_median = median(tbl$period, na.rm = TRUE),
  1437. period_sd = safe_sd(tbl$period, na.rm = TRUE),
  1438. phase_distr_max = density(tbl$phase_norm, bw = 5) |>
  1439. (\(d) max(d$y))(),
  1440. amplitude_min = min(tbl$amplitude, na.rm = TRUE),
  1441. amplitude_max = max(tbl$amplitude, na.rm = TRUE),
  1442. amplitude_median = median(tbl$amplitude, na.rm = TRUE),
  1443. amplitude_sd = safe_sd(tbl$amplitude, na.rm = TRUE),
  1444. error_min = min(tbl$error, na.rm = TRUE),
  1445. error_max = max(tbl$error, na.rm = TRUE),
  1446. error_median = median(tbl$error, na.rm = TRUE),
  1447. error_sd = safe_sd(tbl$error, na.rm = TRUE),
  1448. RAE_min = min(tbl$RAE, na.rm = TRUE),
  1449. RAE_max = max(tbl$RAE, na.rm = TRUE),
  1450. RAE_median = median(tbl$RAE, na.rm = TRUE),
  1451. RAE_sd = safe_sd(tbl$RAE, na.rm = TRUE)
  1452. )
  1453. })
  1455. # collate all into a table
  1456. period_tbl_ranges_df <- do.call(rbind, period_ranges)
  1458. # if AUC gets computed, get all AUC results
  1459. AUC_ranges <- lapply(auc_tbl_list, function(vec) {
  1461. if(!nrow(vec) > 0) return(NULL)
  1462. data.frame(
  1463. AUC_min = min(vec[,2], na.rm = TRUE),
  1464. AUC_max = max(vec[,2], na.rm = TRUE),
  1465. AUC_median = median(vec[,2], na.rm = TRUE),
  1466. AUC_sd = sd(vec[,2], na.rm = TRUE)
  1467. )
  1468. })
  1470. # collate all into a table
  1471. AUC_tbl_ranges_df <- do.call(rbind, AUC_ranges)
  1473. # determine ranges
  1475. period_range <-
  1476. c(
  1477. round(
  1478. (mean(period_tbl_ranges_df$period_min, na.rm = TRUE)
  1479. - (mean(period_tbl_ranges_df$period_sd, na.rm = TRUE))),
  1480. 1),
  1481. round(
  1482. (mean(period_tbl_ranges_df$period_max, na.rm = TRUE)
  1483. + (mean(period_tbl_ranges_df$period_sd, na.rm = TRUE))),
  1484. 1)
  1485. )
  1487. amplitude_range <-
  1488. c(
  1489. round(
  1490. (mean(period_tbl_ranges_df$amplitude_min, na.rm = TRUE)
  1491. -(mean(period_tbl_ranges_df$amplitude_sd, na.rm = TRUE)/2)),
  1492. 0),
  1493. round(
  1494. (median(period_tbl_ranges_df$amplitude_max, na.rm = TRUE)
  1495. + (mean(period_tbl_ranges_df$amplitude_sd, na.rm = TRUE)/2)),
  1496. 0)
  1497. )
  1499. RAE_range <-
  1500. c(
  1501. round(
  1502. (min(period_tbl_ranges_df$RAE_min, na.rm = TRUE)
  1503. -(mean(period_tbl_ranges_df$RAE_sd, na.rm = TRUE)/2)),
  1504. 1),
  1505. round(
  1506. (max(period_tbl_ranges_df$RAE_max, na.rm = TRUE)
  1507. + (mean(period_tbl_ranges_df$RAE_sd, na.rm = TRUE)/2)),
  1508. 1)
  1509. )
  1511. error_range <-
  1512. c(
  1513. round(
  1514. (min(period_tbl_ranges_df$error_min, na.rm = TRUE)
  1515. -(mean(period_tbl_ranges_df$error_sd, na.rm = TRUE)/2)),
  1516. 1),
  1517. round(
  1518. (max(period_tbl_ranges_df$error_max, na.rm = TRUE)
  1519. + (mean(period_tbl_ranges_df$error_sd, na.rm = TRUE)/2)),
  1520. 1)
  1521. )
  1523. AUC_range <-
  1524. c(
  1525. round(
  1526. (min(AUC_tbl_ranges_df$AUC_min, na.rm = TRUE)
  1527. - (mean(AUC_tbl_ranges_df$AUC_sd, na.rm = TRUE))),
  1528. 0),
  1529. round(
  1530. (median(AUC_tbl_ranges_df$AUC_max, na.rm = TRUE)
  1531. + (mean(AUC_tbl_ranges_df$AUC_sd, na.rm = TRUE))),
  1532. 0)
  1533. )
  1535. phase_distr_range <-
  1536. c(0,
  1537. (round(
  1538. max(period_tbl_ranges_df$phase_distr_max, na.rm = TRUE),
  1539. 2)*1.3)
  1540. )
  1542. # assign to list to be added to params
  1543. channel_ranges <- list(
  1544. period_range = period_range,
  1545. amplitude_range = amplitude_range,
  1546. phase_distr_range = phase_distr_range,
  1547. RAE_range = RAE_range,
  1548. error_range = error_range,
  1549. AUC_range = AUC_range
  1550. )
  1551. return(channel_ranges)
  1553. }) # end of purrr::walk
  1554. channel_ranges <- setNames(channel_ranges, names(params$channels))
  1555. return(channel_ranges)
  1556. }
  1558. # create a similarity scores between two tables that share a key variable
  1559. #' Compute a pairwise similarity score between two variables
  1560. #'
  1561. #' @param tbl1 Data frame containing \code{var1} and \code{key}.
  1562. #' @param tbl2 Data frame containing \code{var2} and \code{key}.
  1563. #' @param var1 Character; column name of the first variable in \code{tbl1}.
  1564. #' @param var2 Character; column name of the second variable in \code{tbl2}.
  1565. #' @param key Character; name of the shared ID column used to join the
  1566. #' tables. Defaults to \code{"ID"}.
  1567. #'
  1568. #' @return A data frame with the \code{key} column and a \code{similarity}
  1569. #' column (\code{1 - |var1 - var2|}).
  1570. #' @keywords internal
  1571. similarity_score <- function(tbl1, tbl2, var1, var2, key = "ID") {
  1572. # browser()
  1573. # Ensure both inputs are data frames
  1574. stopifnot(is.data.frame(tbl1), is.data.frame(tbl2))
  1576. # keyvar = !!sym(key)
  1577. # merge the tables by key
  1578. merged_data <- left_join(tbl1, tbl2, by = key) %>%
  1579. `colnames<-`(c(key, var1, var2))
  1581. # Extract the vectors
  1582. x <- merged_data[[var1]]
  1583. y <- merged_data[[var2]]
  1585. # # Check for values outside 0–1
  1586. # if (any(x < 0 | x > 1, na.rm = TRUE) || any(y < 0 | y > 1, na.rm = TRUE)) {
  1587. # warning("One or both variables have values outside [0,1].")
  1588. # }
  1589. #
  1590. # # Remove NA pairs
  1591. # valid <- complete.cases(x, y)
  1592. # if (!all(valid)) {
  1593. # warning(sprintf("Removing %d incomplete rows (NA)", sum(!valid)))
  1594. # x <- x[valid]
  1595. # y <- y[valid]
  1596. # }
  1598. # Compute similarity score
  1599. similarity <- 1 - abs(merged_data[[var1]] - merged_data[[var2]])
  1600. merged_data$similarity <- similarity
  1601. #remove var1 and var2 columns
  1602. merged_data <- merged_data %>% select(-!!sym(var1), -!!sym(var2))
  1603. return(merged_data)
  1604. }
  1606. #' Summarise period analysis results into a single-row table
  1607. #'
  1608. #' @param period_res Named list as returned by \code{computePeriod()}; must
  1609. #' contain \code{period_table} and \code{period_table_unfiltered}.
  1610. #' @param circ_stats Named list as returned by \code{circular_stats()};
  1611. #' must contain \code{vectorLength} and \code{circStats}.
  1612. #' @param Ch_rep Character; channel label used as a column prefix. Defaults
  1613. #' to \code{"Ch"}.
  1614. #' @param coherence Logical; if \code{TRUE}, includes coherence metrics in
  1615. #' the summary. Defaults to \code{FALSE}.
  1616. #'
  1617. #' @return A one-row data frame of summary statistics for the analysis.
  1618. #' @keywords internal
  1619. summarizePeriod <- function(
  1620. period_res,
  1621. circ_stats,
  1622. # filename = "sample",
  1623. # interval_name = "Int",
  1624. Ch_rep = "Ch",
  1625. coherence = FALSE
  1626. ) {
  1628. period_table <- period_res$period_table
  1629. # TODO add check that
  1630. # trace stats
  1631. period_counts <- length(period_table$phase_h)
  1632. non_na_count <- sum(!is.na(period_res$period_table_unfiltered$phase_h))
  1633. if(non_na_count > 0) {
  1635. # period stats
  1636. period_var = stats::var(period_table$period, na.rm = TRUE) |> round(3)
  1637. period_mean = mean(period_table$period, na.rm = TRUE) |> round(2)
  1638. # phase stats
  1639. circ_phase_mean_rad = circular::mean.circular(
  1640. circular(
  1641. period_table$phase_rad,
  1642. units = "rad",
  1643. rotation = "clock",
  1644. template = "none",
  1645. modulo = "asis",
  1646. zero = 0
  1647. ),
  1648. na.rm = TRUE
  1649. ) |> round(3)
  1651. circ_phase_mean = circular::mean.circular(
  1652. circular(
  1653. period_table$phase_circ,
  1654. units = "hour",
  1655. rotation = "clock",
  1656. template = "none",
  1657. modulo = "asis",
  1658. zero = 0
  1659. ),
  1660. na.rm = TRUE
  1661. ) |> round(2)
  1663. phase_rad_var = circular::angular.variance(
  1664. circular(
  1665. period_table$phase_rad,
  1666. units = "radians",
  1667. rotation = "clock",
  1668. template = "none",
  1669. modulo = "asis",
  1670. zero = 0
  1671. ),
  1672. na.rm = TRUE
  1673. ) |> round(3)
  1675. phase_var = circular::angular.variance(
  1676. circular(
  1677. period_table$phase_circ,
  1678. units = "radians",
  1679. rotation = "clock",
  1680. template = "none",
  1681. modulo = "asis",
  1682. zero = 0
  1683. ),
  1684. na.rm = TRUE
  1685. ) |> round(3)
  1687. ' Deprecated
  1688. # first peak phase
  1689. first_peak_phase_mean = circular::mean.circular(
  1690. circular(
  1691. period_table$first_peak_phase_h,
  1692. units = "hour",
  1693. rotation = "clock",
  1694. template = "none",
  1695. modulo = "asis",
  1696. zero = 0
  1697. ),
  1698. na.rm = TRUE
  1699. )
  1701. first_phase_var = circular::angular.variance(
  1702. circular(
  1703. period_table$first_peak_phase_h,
  1704. units = "hour",
  1705. rotation = "clock",
  1706. template = "none",
  1707. modulo = "asis",
  1708. zero = 0
  1709. ),
  1710. na.rm = TRUE
  1711. )
  1712. '
  1714. # other stats
  1715. vector_length <- ifelse(!is.null(circ_stats$vectorLength),
  1716. round(circ_stats$vectorLength, 3),
  1717. NA)
  1718. amplitude_mean <- mean(abs(period_table$amplitude), na.rm = TRUE) |> round(2)
  1719. amplitude_var <- stats::var(period_table$amplitude, na.rm = TRUE) |> round(3)
  1720. RAE_mean <- mean(period_table$RAE, na.rm = TRUE) |> round(2)
  1721. RAE_var <- stats::var(period_table$RAE, na.rm = TRUE) |> round(2)
  1722. error_mean <- mean(period_table$error, na.rm = TRUE) |> round(2)
  1724. } else {
  1725. period_var = NA
  1726. period_mean = NA
  1727. circ_phase_mean_rad = NA
  1728. phase_rad_var = NA
  1729. vector_length = NA
  1730. circ_phase_mean = NA
  1731. phase_var = NA
  1732. amplitude_mean = NA
  1733. amplitude_var = NA
  1734. RAE_mean = NA
  1735. RAE_var = NA
  1736. error_mean = NA
  1737. }
  1738. # return values outside
  1739. outreturn_df <- data.frame(
  1740. # filename = as.factor(filename),
  1741. # interval = as.factor(interval_name),
  1743. channel_label = as.factor(Ch_rep), # channel label
  1744. trace_no = non_na_count, # number of rhythmic cells
  1745. trace_tot = period_counts, # number of cells sampled
  1746. period_var = period_var, # variance of period values
  1747. period_mean = period_mean, # mean of period values in hours
  1748. phase_rad_mean = circ_phase_mean_rad, # mean of period values in radians
  1749. phase_rad_var = phase_rad_var, # circular variance of phases in radians
  1750. vector_length = vector_length, # vector length of rayleigh plot
  1751. circ_phase_mean = circ_phase_mean, # circular mean of circadian phase values in hours
  1752. circ_phase_var = phase_var, # circular variance of circadian phases
  1753. # first_peak_phase_mean = first_peak_phase_mean,
  1754. # first_phase_var = first_phase_var,
  1755. amplitude_mean = amplitude_mean,
  1756. amplitude_var = amplitude_var,
  1757. RAE_mean = RAE_mean,
  1758. RAE_var = RAE_var,
  1759. error_mean = error_mean
  1760. )
  1762. if(coherence){
  1763. # spatial coherence stats
  1765. #period
  1766. period_crn_var = stats::var(period_table$period_crn, na.rm = TRUE) |>
  1767. round(3)
  1768. period_crn_mean = mean(period_table$period_crn, na.rm = TRUE) |>
  1769. round(3)
  1771. #phase
  1772. phase_crn_var = stats::var(period_table$phase_crn, na.rm = TRUE) |>
  1773. round(3)
  1774. phase_crn_mean = mean(period_table$phase_crn, na.rm = TRUE) |>
  1775. round(3)
  1777. #amplitude
  1778. amplitude_crn_var = stats::var(period_table$amp_crn, na.rm = TRUE) |>
  1779. round(3)
  1780. amplitude_crn_mean = mean(period_table$amp_crn, na.rm = TRUE) |>
  1781. round(3)
  1783. # add to table
  1784. outreturn_df <- cbind(outreturn_df,
  1785. period_crn_mean, period_crn_var,
  1786. phase_crn_mean, phase_crn_var,
  1787. amplitude_crn_mean, amplitude_crn_var)
  1788. }
  1790. return(outreturn_df)
  1791. }

processing_foo.R at commit 0c2d463, under other · at the source

Overview

Authors: Marco Ferrari1,2, Natalie Ness1,2, Julieta Acosta1,2, Marco Brancaccio1,2
ORCID iDs: Natalie Ness, Marco Brancaccio
  1. UK Dementia Research Institute at Imperial College London London United Kingdom
  2. Department of Brain Science Imperial College London London United Kingdom
Institutions: UK Dementia Research Institute (United Kingdom); Imperial College London (United Kingdom)
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany), volume 13, issue 40, article e75427
Dates: received 3 February 2026; accepted 11 April 2026; published online 28 April 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/advs.75427 · PMID 42048014 · PMCID PMC13335426 · OpenAlex W7158352060
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism)
Methods: Statistics, Connectivity, Preprocessing, Spectral & time-frequency, fMRI & imaging
Keywords: brain tissue, circadian rhythms, live‐imaging
MeSH: Brain*, Circadian Rhythm*, Neurons*, Suprachiasmatic Nucleus*, Animals, ARNTL Transcription Factors, Astrocytes, Calcium, Mice (* major topic)
Topic: Circadian rhythm and melatonin (Endocrine and Autonomic Systems, Neuroscience), according to OpenAlex
Funding: UK Dementia Research Institute (UKDRI‐5007); Michael Uren Foundation
Citations: not cited yet (Europe PMC); 62 references in the paper
Research resources: RRID:IMSR_JAX:006852

Abstract

Circadian function in multicellular organisms arises from coordinated interactions amongst diverse cellular tissue populations. Existing approaches for long‐term imaging of within‐tissue circadian regulation remain low‐throughput, highly specialized, and largely inaccessible. Here, we developed ClockCyte, a high‐content fluorescent live‐imaging platform that enables continuous monitoring of circadian rhythms in up to 144 brain tissue samples. Using the mouse suprachiasmatic nucleus as a model, ClockCyte captures the differential circadian tissue regulation of neurons and astrocytes. We further identified a previously uncharacterized oscillatory circadian compartment in axonal calcium, showing highly homogeneous activity, opposed to waves of intracellular neuronal calcium. By deleting Bmal1 in neurons, we reveal the network underpinnings connecting clock gene expression to network‐wide axonal regulation. The discovery of distinct circadian properties of axonal calcium and their disruption by Bmal1 ablation highlights the potential to reveal new principles of intra‐tissue network‐level circadian organization. More broadly, this approach will enable systematic explorations of how cell‐type‐specific and compartmentalized subcellular rhythms contribute to brain physiology.

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

Repository

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

cabaJr/ClockCyteR.spatial

License: other
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 0c2d463886c3398bb6e55d8211cc4e60c7f08a55, 1 May 2026
Languages: R (19)
Size: 141 files, 19 scripts
Software Heritage: not archived
Found in: the text, “Spatiotemporal Circadian Activity Analysis–Clock”
Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 4 notebooks
Not found: CITATION.cff
Tools: tidyverse (10 files), ggplot2 (3 files), data.table (2 files), igraph (2 files), ggpubr (1 file), Plotly (1 file), rstatix (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
21 files
At the source: github.com/cabaJr/ClockCyteR.spatial

The paper's code and data availability statement is in the Data section.

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 19 scripts, each with its path and the digest of its content;
  • 8 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 Statement

The ClockCyteR and ClockCyteR.spatial packages are currently deposited on GitHub repositories (https://github.com/cabaJr/ClockCyteR) (https://github.com/cabaJr/ClockCyteR.spatial). The raw data used in the generation of the figures are available upon request.

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 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 3 keywords, 9 MeSH terms, 2 funders, 53 references, 1 RRID.

Cite

This paper

Ferrari, M., Ness, N., Acosta, J., & Brancaccio, M. (2026). A High-Throughput Live Imaging Platform to Investigate Circuit-Dependent Regulation of Circadian Rhythms in Brain Tissue. Advanced science (Weinheim, Baden-Wurttemberg, Germany), 13(40), e75427. https://doi.org/10.1002/advs.75427

BibTeX

@article{ferrari2026high,
author = {Ferrari, Marco and Ness, Natalie and Acosta, Julieta and Brancaccio, Marco},
title = {{A High-Throughput Live Imaging Platform to Investigate Circuit-Dependent Regulation of Circadian Rhythms in Brain Tissue}},
journal = {Advanced science (Weinheim, Baden-Wurttemberg, Germany)},
year = {2026},
month = apr,
volume = {13},
number = {40},
pages = {e75427},
publisher = {Wiley},
issn = {2198-3844},
doi = {10.1002/advs.75427},
url = {https://doi.org/10.1002/advs.75427},
pmid = {42048014},
pmcid = {PMC13335426}
}

RIS

TY - JOUR
AU - Ferrari, Marco
AU - Ness, Natalie
AU - Acosta, Julieta
AU - Brancaccio, Marco
TI - A High-Throughput Live Imaging Platform to Investigate Circuit-Dependent Regulation of Circadian Rhythms in Brain Tissue
T2 - Advanced science (Weinheim, Baden-Wurttemberg, Germany)
J2 - Adv Sci (Weinh)
PY - 2026
DA - 2026/04/28
VL - 13
IS - 40
SP - e75427
SN - 2198-3844
PB - Wiley
DO - 10.1002/advs.75427
UR - https://doi.org/10.1002/advs.75427
LA - en
ER -

CSL-JSON

{
"id": "10.1002/advs.75427",
"type": "article-journal",
"title": "A High-Throughput Live Imaging Platform to Investigate Circuit-Dependent Regulation of Circadian Rhythms in Brain Tissue",
"container-title": "Advanced science (Weinheim, Baden-Wurttemberg, Germany)",
"author": [
{
"family": "Ferrari",
"given": "Marco"
},
{
"family": "Ness",
"given": "Natalie"
},
{
"family": "Acosta",
"given": "Julieta"
},
{
"family": "Brancaccio",
"given": "Marco"
}
],
"container-title-short": "Adv Sci (Weinh)",
"volume": "13",
"issue": "40",
"page": "e75427",
"DOI": "10.1002/advs.75427",
"PMID": "42048014",
"PMCID": "PMC13335426",
"ISSN": "2198-3844",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/advs.75427",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
28
]
]
}
}

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.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: rstatix, igraph, Plotly, 4 other tools, mouse
[2] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: rstatix, igraph, Plotly, 4 other tools, mouse
[3] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: rstatix, igraph, Plotly, 4 other tools
[4] doi:10.1038/s41467-026-74753-y [code]
A human-specific microRNA controls the timing of excitatory synaptogenesis.
Journal: Nature communications
In common: rstatix, igraph, Plotly, 4 other tools
[5] doi:10.1073/pnas.2609132123 [code]
A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: rstatix, igraph, Plotly, 4 other tools
[6] doi:10.1126/sciadv.aeg3223 [code]
The extreme diversity of retinal amacrine cells has deep evolutionary roots.
Journal: Science advances
In common: rstatix, igraph, Plotly, 3 other tools, 1 reference
[7] doi:10.1371/journal.pcbi.1014573 [code]
Cell-type-specific m1A dynamics are associated with microglial phenotypic transition and neuronal metabolic adaptation during spinal cord injury.
Journal: PLoS computational biology
In common: rstatix, igraph, ggpubr, 3 other tools, mouse
[8] doi:10.1038/s41586-026-10755-6 [code]
An intrinsic cytoskeletal oscillator establishes neuronal polarity.
Journal: Nature
In common: rstatix, Plotly, ggpubr, 3 other tools, mouse
[9] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: igraph, Plotly, ggpubr, 3 other tools, mouse
[10] doi:10.1073/pnas.2523130123 [code]
FABP7 controls radial glial scaffold stability during human cortical development.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: igraph, Plotly, ggpubr, 3 other tools, mouse

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.