OSCR

Diurnal choroid plexus function in mice depends on sex, age, and amyloid-β status.

Code ↔ Paper

9 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 9 matches
  1. [1] § Results › Diurnal regulation of Aqp1 and NKCC1 at the transcript level in the aged ChP ↔ deseq2 feb 2023.Rmd, lines 1088–1113 · score 1.00 · carbonic anhydrases, ATP synthase, Atp5b, Slc12a4, Slc12a7, Slc25a28
  2. [2] § Results › Comparison of gene expression across different datasets ↔ deseq2 feb 2023.Rmd, lines 1088–1113 · score 0.72 · Slc22a8, Nr1d2, Nr1d1, Bhlhe41, Ciart, Per3
  3. [3] § Methods › Transcriptomic analysis ↔ deseq2 feb 2023.Rmd, lines 49–163 · score 0.68 · AnnotationDbi, DESeq2, log2 fold change, pheatmap, threshold, dplyr
  4. [4] § Results › Aged mice exhibit altered diurnal transcriptional profile ↔ CP circadian gene validation.Rmd, lines 491–555 · score 0.67 · circadian related genes, nr1d2, nr1d1, Cry, amplitude, mesor
  5. [5] § Methods › Statistics and reproducibility ↔ CP circadian gene validation.Rmd, lines 84–215 · score 0.65 · Kruskall Wallis, Shapiro Wilk
  6. [6] § Results › Young mice exhibit diurnal, transcriptional regulation of the ChP, that is not significantly impacted by acute sleep disruption ↔ CP circadian gene validation.Rmd, lines 491–555 · score 0.60 · nr1d2, nr1d1, related genes, Cry, validated, Dbp
  7. [7] § Results › Comparison of gene expression across different datasets ↔ genecompare.R, lines 891–976 · score 0.54 · biological processes, Molecular function, downregulated, pseudotime, ZT, compartment
  8. [8] § Results › Diurnal regulation of Aqp1 and NKCC1 at the transcript level in the aged ChP ↔ deseq2 feb 2023.Rmd, lines 1190–1317 · score 0.54 · slc12a2, log2 fold change, related genes, junctional, amyloid, proteins
  9. [9] § Results › Comparison of gene expression across different datasets ↔ genecompare.R, lines 891–976 · score 0.53 · cellular compartment, biological processes, molecular function, Pathway, clusterProfiler, young

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 Markdown · 1,396 lines · 58 KB · CC-BY-4.0 · 4 matches

  1. ---
  2. title: "DESEQ CP circadian data"
  3. output: html_document
  4. date: "2023-02-22"
  5. editor_options:
  6. chunk_output_type: console
  7. ---
  9. #Analysis of RNAseq data from sleep/circadian choroid plexus
  10. ###There are several different ways to analyze this data as there are 3 conditions, awake, asleep and sleep disrupted. Previous analysis showed no significant difference between sleep and sleep disrupted although there was a difference in the sig genes that were changed compared to awake state. So will explore all these options.
  11. ###Note: Awake (AW) is Night, sleep (SL) and sleep disrupted (SD) are Day. comparison will be done Awake vs SL or SD, OR Night vs Day
  13. Install required packages
  15. if (!require("BiocManager", quietly = TRUE))
  16. install.packages("BiocManager")
  18. BiocManager::install("DESeq2")
  19. BiocManager::install("org.Mm.eg.db")
  20. install.packages()
  21. install.packages("Hmisc")
  22. BiocManager::install("ComplexHeatmap")
  24. if (!require("BiocManager", quietly = TRUE))
  25. install.packages("BiocManager")
  27. BiocManager::install("EnhancedVolcano")
  28. install.packages("vsn")
  29. setwd("~/Documents/Iliff lab /Experiments/Expt 13 sleep study mouse choroid plexus/Circadian R folder/deseq2 feb2023 r markdown file output")
  32. ```{r setup, include=FALSE}
  33. knitr::opts_chunk$set(echo = TRUE)
  35. library(ComplexHeatmap)
  36. library(pheatmap)
  37. library(tidyverse)
  38. library(dplyr)
  39. library(DESeq2)
  40. library(RColorBrewer)
  41. library(AnnotationDbi)
  42. library(org.Mm.eg.db)
  43. library(reshape)
  44. library(viridis)
  45. library(Hmisc)
  46. library(EnhancedVolcano)
  47. ```
  49. #QC of the samples
  50. ```{r }
  51. dds_raw<-read.csv("3mocounts.csv", row.names = 1)
  52. metadata<-read.csv("3mometadata.csv", row.names = 1)
  53. metadata$Condition<-factor(metadata$Condition, levels=c("AW", "SL", "SD"))
  54. metadata$TimeOfDay<-factor(metadata$TimeOfDay, levels=c("Night", "Day"))
  55. View(metadata)
  57. #explore
  58. rownames(metadata)
  59. colnames(dds_raw)
  60. #make sure the rows and columns match
  61. all(rownames(metadata) == colnames(dds_raw))
  63. #creating the DESeq2 object, using sleep/wake state
  64. dds_raw<-DESeqDataSetFromMatrix(countData = dds_raw, colData = metadata, design = ~ Condition )
  67. #pre-filtering
  68. #smallestGroupSize<-4
  69. #keep<-rowSums(counts(dds_raw) >= smallestGroupSize)
  70. #dds_raw<-dds_raw[keep,]
  71. #cuts from 29179 genes to 17701- this is too strict as circadian differences are subtle. if using this cutoff no sig DEGs
  73. #run DESeq2
  74. dds_raw<-DESeq(dds_raw)
  75. res_raw<-results(dds_raw)
  76. res_raw
  78. #determine size factors to use for normalization
  79. dds_raw <- estimateSizeFactors(dds_raw)
  80. dds_raw <- estimateDispersions(dds_raw)
  81. dds_raw <- nbinomWaldTest(dds_raw)
  83. sizeFactors(dds_raw)
  85. #extract the normalized counts
  86. normalized_counts <-counts(dds_raw, normalized = TRUE)
  88. ###Getting the ENTREZID as rownames for the damn dataframe.
  89. normalized_counts <- data.frame(ENTREZID=row.names(normalized_counts), normalized_counts)
  90. View(normalized_counts)
  92. write.csv(normalized_counts, file="normalized_counts_3mo_all.csv")
  94. #perform unsupervised clustering analysis:log transformation
  95. nrow(dds_raw)
  96. vsd_dds_raw <- vst(dds_raw, blind = TRUE)
  97. vsd_mat_dds_raw <- assay(vsd_dds_raw)
  98. vsd_cor_dds_raw <- cor(vsd_mat_dds_raw)
  99. View(vsd_cor_dds_raw)
  101. pheatmap(vsd_cor_dds_raw, cluster_rows = TRUE, cluster_cols = TRUE)
  102. pheatmap(vsd_cor_dds_raw, annotation = dplyr::select(metadata, TimeOfDay))
  103. pheatmap(vsd_cor_dds_raw, annotation = dplyr::select(metadata, Condition))
  105. plotPCA(vsd_dds_raw, intgroup = "Condition")
  106. plotPCA(vsd_dds_raw, intgroup=c("Condition", "TimeOfDay"))
  109. ###############################################################
  110. padj.cutoff <- 0.05
  111. log2FoldChange.cutoff <- 0.58
  112. dds_raw_AW_SL<-dds_raw
  113. design(dds_raw_AW_SL)<-formula(~Condition)
  114. dds_raw_AW_SL$Condition<-factor(dds_raw_AW_SL$Condition, levels=c("AW", "SL", "SD"))
  115. dds_raw_AW_SL<-DESeq(dds_raw_AW_SL)
  116. res_rawAWvsSL<-results(dds_raw_AW_SL, contrast = c("Condition", "SL", "AW"))
  118. res_rawAWvsSL_df <- as.data.frame(res_rawAWvsSL)
  119. res_rawAWvsSL_df<-tibble::rownames_to_column(res_rawAWvsSL_df, "ENTREZID")
  120. res_rawAWvsSL_df$ENTREZID<- as.character(res_rawAWvsSL_df$ENTREZID)
  121. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_rawAWvsSL_df$ENTREZID,
  122. columns=c("SYMBOL","GENENAME"),
  123. keytype="ENTREZID")
  124. res_rawAWvsSL_df <- left_join(res_rawAWvsSL_df,anno,by="ENTREZID")
  125. #for some reason the annotation put gene ID 11865 to Bmal instead of Arntl so need to replace this gene and make consistent with the old/AD datasets
  126. res_rawAWvsSL_df$SYMBOL[res_rawAWvsSL_df$ENTREZID == "11865"]<-"Arntl"
  127. res_rawAWvsSL_df$GENENAME[res_rawAWvsSL_df$ENTREZID == "11865"]<-"aryl hydrocarbon receptor nuclear translocator-like"
  129. save(res_rawAWvsSL_df, file="raw_AWvsSL_3mo.Rdata")
  131. threshold <- res_rawAWvsSL_df$padj < padj.cutoff & abs(res_rawAWvsSL_df$log2FoldChange) > log2FoldChange.cutoff
  132. length(which(threshold))
  133. res_rawAWvsSL_df$threshold <- threshold
  134. sig_res_rawAWvsSL <- data.frame(subset(res_rawAWvsSL_df, threshold==TRUE), row.names = 1)
  135. dim(sig_res_rawAWvsSL)
  136. sig_res_rawAWvsSL<-na.omit(sig_res_rawAWvsSL)
  137. dim(sig_res_rawAWvsSL)
  138. write.csv(sig_res_rawAWvsSL, file="sig_res_rawAWvsSL.csv")
  139. #save(sig_res_3moAWvsSL, file="sigrawAWvsSL_3mo.Rdata")
  141. par(mfrow=c(2,3))
  142. plotCounts(dds_raw_AW_SL, "217166", intgroup = "TimeOfDay", main = "Nr1d1")
  143. plotCounts(dds_raw_AW_SL, "11865", intgroup = "TimeOfDay", main = "Arntl")
  144. plotCounts(dds_raw_AW_SL, "12952", intgroup = "TimeOfDay", main = "Cry1")
  145. plotCounts(dds_raw_AW_SL, "233908", intgroup = "TimeOfDay", main = "Fus")
  146. plotCounts(dds_raw_AW_SL, "74840", intgroup = "TimeOfDay", main = "Manf")
  147. plotCounts(dds_raw_AW_SL, "18030", intgroup = "TimeOfDay", main = "Nfil3")
  149. plotCounts(dds_raw_AW_SL, "217166", intgroup = "Condition", main = "Nr1d1")
  150. plotCounts(dds_raw_AW_SL, "11865", intgroup = "Condition", main = "Arntl")
  151. plotCounts(dds_raw_AW_SL, "12952", intgroup = "Condition", main = "Cry1")
  152. plotCounts(dds_raw_AW_SL, "233908", intgroup = "Condition", main = "Fus")
  153. plotCounts(dds_raw_AW_SL, "74840", intgroup = "Condition", main = "Manf")
  154. plotCounts(dds_raw_AW_SL, "18030", intgroup = "Condition", main = "Nfil3")
  156. assays(dds_raw)[["cooks"]]
  157. summary(res_raw)
  159. par(mar=c(8,5,2,2))
  160. boxplot(log10(assays(dds_raw)[["cooks"]]), range=0, las=2)
  161. #removing sample #12 going forward as it is an outlier by PCA, and failed RNAseq QC
  163. ```
  166. ```{r running DESEQ2 on count data obtained from galaxy alignment}
  168. #make a dds file from galaxy generated normalized counts (rlog)
  169. #sample 12 is an outlier so remove from anlaysis (from previous analysis)
  170. dds_3mo<-read.csv("3mocounts.csv", row.names = 1)
  171. dds_3mo<-dds_3mo[ , c(-7)]
  172. View(dds_3mo)
  174. metadata<-read.csv("3mometadata.csv", row.names = 1)
  175. metadata<-metadata[c(-7), ]
  176. metadata$Condition<-factor(metadata$Condition, levels=c("AW", "SL", "SD"))
  177. metadata$TimeOfDay<-factor(metadata$TimeOfDay, levels=c("Night", "Day"))
  178. View(metadata)
  180. #explore
  181. rownames(metadata)
  182. colnames(dds_3mo)
  183. #make sure the rows and columns match
  184. all(rownames(metadata) == colnames(dds_3mo))
  187. ```
  189. #Comparing the AW vs SL or SD, including all 3 conditions for now
  190. ```{r}
  191. #creating the DESeq2 object, using *sleep condition as condition
  192. dds_3mo_cond<-DESeqDataSetFromMatrix(countData = dds_3mo, colData = metadata, design = ~ Condition)
  193. #determine size factors to use for normalization
  194. dds_3mo_cond <- estimateSizeFactors(dds_3mo_cond)
  196. sizeFactors(dds_3mo_cond)
  197. #extract the normalized counts
  198. normalized_counts_3mo_cond <-counts(dds_3mo_cond, normalized = TRUE)
  200. View(normalized_counts_3mo)
  202. write.csv(normalized_counts_3mo_cond, file="normalized_counts_3mo_cond_no12.csv")
  204. ###Getting the ENTREZID as rownames for the damn dataframe.
  205. normalized_counts_3mo_cond_df <- data.frame(ENTREZID=row.names(normalized_counts_3mo_cond), normalized_counts_3mo_cond)
  207. #perform unsupervised clustering analysis:log transformation
  208. nrow(dds_3mo_cond)
  209. vsd_3mo_cond <- vst(dds_3mo_cond, blind = TRUE)
  210. vsd_mat_3mo_cond <- assay(vsd_3mo_cond)
  211. vsd_cor_3mo_cond <- cor(vsd_mat_3mo_cond)
  212. View(vsd_cor_3mo_cond)
  213. ```
  215. ```{r}
  216. #DEG
  217. dds_3mo_cond<-DESeq(dds_3mo_cond)
  218. res<-results(dds_3mo_cond, contrast = c("Condition", "SD", "SL")) #this should give AW vs SL AW vs SD and SD vs SL
  219. ```
  221. Time of Day analysis comparing the AW animals to the grouped SL and SD animals together
  222. ```{r}
  223. #creating the DESeq2 object, using *TimeOfDay as condition
  224. dds_3mo<-DESeqDataSetFromMatrix(countData = dds_3mo, colData = metadata, design = ~ TimeOfDay)
  225. #determine size factors to use for normalization
  226. dds_3mo <- estimateSizeFactors(dds_3mo)
  228. sizeFactors(dds_3mo)
  229. #extract the normalized counts
  230. normalized_counts_3mo <-counts(dds_3mo, normalized = TRUE)
  232. View(normalized_counts_3mo)
  234. write.csv(normalized_counts_3mo, file="normalized_counts_3mo_ToD_no12.csv")
  236. ###Getting the ENTREZID as rownames for the damn dataframe.
  237. normalized_counts_3mo_df <- data.frame(ENTREZID=row.names(normalized_counts_3mo), normalized_counts_3mo)
  239. #perform unsupervised clustering analysis:log transformation
  240. nrow(dds_3mo)
  241. vsd_3mo <- vst(dds_3mo, blind = TRUE)
  242. vsd_mat_3mo <- assay(vsd_3mo)
  243. vsd_cor_3mo <- cor(vsd_mat_3mo)
  244. View(vsd_cor_3mo)
  245. ```
  247. #Running the DESEQ2 with multifactorial effect of timeofday with condition as covariate
  248. ```{r}
  249. #create a copy of the DESeqDataSet to run the multi-factor design
  250. ddsMFtime<-dds_3mo
  251. levels(ddsMFtime$Condition)
  252. #make SD to AW (as this is their sleep state) to be able to have both condition and time of day as factors (no samples can be have unique factors for each variable)
  253. levels(ddsMFtime$Condition)<-sub("SD", "AW", levels(ddsMFtime$Condition))
  254. design(ddsMFtime)<-formula(~ Condition + TimeOfDay)
  255. ddsMFtime<-DESeq(ddsMFtime)
  257. resMFtime<-results(ddsMFtime)
  258. head(ddsMFtime)
  260. par(mfrow=c(2,3))
  261. plotCounts(ddsMFtime, "217166", intgroup = "TimeOfDay", main = "Nr1d1")
  262. plotCounts(ddsMFtime, "11865", intgroup = "TimeOfDay", main = "Arntl")
  263. plotCounts(ddsMFtime, "12952", intgroup = "TimeOfDay", main = "Cry1")
  264. plotCounts(ddsMFtime, "233908", intgroup = "TimeOfDay", main = "Fus")
  265. plotCounts(ddsMFtime, "74840", intgroup = "TimeOfDay", main = "Manf")
  266. plotCounts(ddsMFtime, "18030", intgroup = "TimeOfDay", main = "Nfil3")
  267. ```
  268. #Running the DESEQ2 with multifactorial effect of condition with time of day as covariate
  269. ```{r}
  270. #create a copy of the DESeqDataSet to run the multi-factor design
  271. ddsMFcond<-dds_3mo
  272. levels(ddsMFcond$Condition)
  273. levels(ddsMFcond$Condition)<-sub("SD", "AW", levels(ddsMFcond$Condition))
  274. #this one is doing the comparison of timeofday taking into account the condition.
  275. #later we'll do the other way around
  276. design(ddsMFcond)<-formula(~TimeOfDay + Condition)
  277. ddsMFcond<-DESeq(ddsMFcond)
  279. resMFcond<-results(ddsMFcond)
  280. head(resMFcond)
  282. par(mfrow=c(2,3))
  283. plotCounts(ddsMFcond, "217166", intgroup = "Condition", main = "Nr1d1")
  284. plotCounts(ddsMFcond, "11865", intgroup = "Condition", main = "Arntl")
  285. plotCounts(ddsMFcond, "12952", intgroup = "Condition", main = "Cry1")
  286. plotCounts(ddsMFcond, "233908", intgroup = "Condition", main = "Fus")
  287. plotCounts(ddsMFcond, "74840", intgroup = "Condition", main = "Manf")
  288. plotCounts(ddsMFcond, "18030", intgroup = "Condition", main = "Nfil3")
  289. ```
  291. ## Visualizing sample distribution
  292. ```{r pressure, echo=FALSE}
  293. pheatmap(vsd_cor_3mo, cluster_rows = TRUE, cluster_cols = TRUE)
  294. pheatmap(vsd_cor_3mo, annotation = dplyr::select(metadata, TimeOfDay))
  295. plotPCA(vsd_3mo, intgroup = "Condition")
  296. plotPCA(vsd_3mo, intgroup=c("Condition", "TimeOfDay"))
  297. ```
  299. ###Visualizing the data with norm counts and heatmaps from top genes
  300. ```{r}
  301. ### Set thresholds
  302. padj.cutoff <- 0.05
  303. log2FoldChange.cutoff <- 0.58
  304. #need to make into data frame to get gene symbols
  306. library(tibble)
  307. resMF_df <- as.data.frame(resMFtime)
  308. resMF_df<-tibble::rownames_to_column(resMF_df, "ENTREZID")
  309. resMF_df$ENTREZID<- as.character(resMF_df$ENTREZID)
  310. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=resMF_df$ENTREZID,
  311. columns=c("SYMBOL","GENENAME"),
  312. keytype="ENTREZID")
  313. resMF_df <- left_join(resMF_df,anno,by="ENTREZID")
  315. threshold <- resMF_df$padj < padj.cutoff & abs(resMF_df$log2FoldChange) > log2FoldChange.cutoff
  316. length(which(threshold))
  317. resMF_df$threshold <- threshold
  318. sig_resMF <- data.frame(subset(resMF_df, threshold==TRUE), row.names = 1)
  319. dim(sig_resMF)
  320. write.csv(sig_resMF, file="sig_resMF.csv")
  321. save(sig_resMF, file="sigToDMF_3mo.Rdata")
  323. ###For condition accounting for ToD
  324. resMFcond_df <- as.data.frame(resMFcond)
  325. resMFcond_df<-tibble::rownames_to_column(resMFcond_df, "ENTREZID")
  326. resMFcond_df$ENTREZID<- as.character(resMFcond_df$ENTREZID)
  327. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=resMFcond_df$ENTREZID,
  328. columns=c("SYMBOL","GENENAME"),
  329. keytype="ENTREZID")
  330. resMFcond_df <- left_join(resMFcond_df,anno,by="ENTREZID")
  332. threshold <- resMFcond_df$padj < padj.cutoff & abs(resMFcond_df$log2FoldChange) > log2FoldChange.cutoff
  333. length(which(threshold))
  334. resMFcond_df$threshold <- threshold
  335. sig_resMFcond <- data.frame(subset(resMFcond_df, threshold==TRUE), row.names = 1)
  336. dim(sig_resMFcond)
  337. save(sig_resMFcond, file="sigcondMF_3mo.Rdata")
  338. #actually nothing is significant so it gives 0 results
  340. ```
  342. ```{r volcano plot}
  343. #looking at timeofday differences accounting for condition
  344. sig_resMF %>% drop_na("SYMBOL")
  345. EnhancedVolcano(resMF_df,
  346. lab = resMF_df$SYMBOL,
  347. x = 'log2FoldChange',
  348. y = 'padj',
  349. title="Night vs Day (accounting for sleep condition)",
  350. subtitle = NULL,
  351. xlim = c(-2.5,2.5),
  352. col=c('grey','black','purple' , 'orange1'),
  353. pCutoff = 10e-6,
  354. FCcutoff = 0.58,
  355. pointSize = 2.0,
  356. labSize = 3,
  357. axisLabSize=10)
  358. #saved as 7x7 pdf "Enhanced_volcano_young_timeofday_accounting_for_condition
  359. ```
  360. #This would be to run the samples ignoring the Condition, since when comparing condition accounting for time of day there are no significant differences
  361. ```{r}
  362. #DEG
  363. dds_3mo<-DESeq(dds_3mo)
  364. res<-results(dds_3mo, contrast = c("TimeOfDay", "Night", "Day"))
  365. ```
  367. ```{r}
  368. res_df <- as.data.frame(res)
  369. res_df<-tibble::rownames_to_column(res_df, "ENTREZID")
  370. res_df$ENTREZID<- as.character(res_df$ENTREZID)
  371. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_df$ENTREZID,
  372. columns=c("SYMBOL","GENENAME"),
  373. keytype="ENTREZID")
  374. res_df <- left_join(res_df,anno,by="ENTREZID")
  376. threshold <- res_df$padj < padj.cutoff & abs(res_df$log2FoldChange) > log2FoldChange.cutoff
  377. length(which(threshold))
  378. res_df$threshold <- threshold
  379. sig_res <- data.frame(subset(res_df, threshold==TRUE), row.names = 1)
  380. dim(sig_res)
  381. save(sig_res, file="sigToD_3mo.Rdata")
  382. ```
  384. ```{r volcano plot}
  385. #plot ignoring the condition variable and just comparing time of day effects
  387. EnhancedVolcano(res_df,
  388. lab = res_df$SYMBOL,
  389. x = 'log2FoldChange',
  390. y = 'padj',
  391. title="Night vs Day (ignoring sleep condition)",
  392. subtitle = NULL,
  393. xlim = c(-2.5,2.5),
  394. col=c('grey','black','purple' , 'orange1'),
  395. pCutoff = 10e-6,
  396. FCcutoff = 0.58,
  397. pointSize = 2.0,
  398. labSize = 3,
  399. axisLabSize=10)
  400. #saved as pdf 7x7 "Enhanced_volcano_young_timeofday_ignoring_for_condition"
  401. ```
  402. #DEG analysis of SL vs AW and SD vs AW
  403. ```{r}
  404. dds_3mo_AW_SL<-dds_3mo
  405. design(dds_3mo_AW_SL)<-formula(~Condition)
  406. dds_3mo_AW_SL$Condition<-factor(dds_3mo_AW_SL$Condition, levels=c("AW", "SL", "SD"))
  407. dds_3mo_AW_SL<-DESeq(dds_3mo_AW_SL)
  408. res_3moAWvsSL<-results(dds_3mo_AW_SL, contrast = c("Condition", "SL", "AW"))
  410. res_3moAWvsSL_df <- as.data.frame(res_3moAWvsSL)
  411. res_3moAWvsSL_df<-tibble::rownames_to_column(res_3moAWvsSL_df, "ENTREZID")
  412. res_3moAWvsSL_df$ENTREZID<- as.character(res_3moAWvsSL_df$ENTREZID)
  413. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_3moAWvsSL_df$ENTREZID,
  414. columns=c("SYMBOL","GENENAME"),
  415. keytype="ENTREZID")
  416. res_3moAWvsSL_df <- left_join(res_3moAWvsSL_df,anno,by="ENTREZID")
  417. #for some reason the annotation put gene ID 11865 to Bmal instead of Arntl so need to replace this gene and make consistent with the old/AD datasets
  418. res_3moAWvsSL_df$SYMBOL[res_3moAWvsSL_df$ENTREZID == "11865"]<-"Arntl"
  419. res_3moAWvsSL_df$GENENAME[res_3moAWvsSL_df$ENTREZID == "11865"]<-"aryl hydrocarbon receptor nuclear translocator-like"
  421. save(res_3moAWvsSL_df, file="allAWvsSL_3mo.Rdata")
  423. threshold <- res_3moAWvsSL_df$padj < padj.cutoff & abs(res_3moAWvsSL_df$log2FoldChange) > log2FoldChange.cutoff
  424. length(which(threshold))
  425. res_3moAWvsSL_df$threshold <- threshold
  426. sig_res_3moAWvsSL <- data.frame(subset(res_3moAWvsSL_df, threshold==TRUE), row.names = 1)
  427. dim(sig_res_3moAWvsSL)
  428. sig_res_3moAWvsSL<-na.omit(sig_res_3moAWvsSL)
  429. dim(sig_res_3moAWvsSL)
  430. write.csv(sig_res_3moAWvsSL, file="sig_res_3moAWvsSL.csv")
  431. save(sig_res_3moAWvsSL, file="sigAWvsSL_3mo.Rdata")
  432. ```
  433. #SD vs AW
  434. ```{r}
  435. dds_3mo_AW_SD<-dds_3mo
  436. design(dds_3mo_AW_SD)<-formula(~Condition)
  437. dds_3mo_AW_SD$Condition<-factor(dds_3mo_AW_SD$Condition, levels=c("AW", "SL", "SD"))
  438. dds_3mo_AW_SD<-DESeq(dds_3mo_AW_SD)
  439. res_3moAWvsSD<-results(dds_3mo_AW_SD, contrast = c("Condition", "SD", "AW"))
  441. res_3moAWvsSD_df <- as.data.frame(res_3moAWvsSD)
  442. res_3moAWvsSD_df<-tibble::rownames_to_column(res_3moAWvsSD_df, "ENTREZID")
  443. res_3moAWvsSD_df$ENTREZID<- as.character(res_3moAWvsSD_df$ENTREZID)
  444. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_3moAWvsSD_df$ENTREZID,
  445. columns=c("SYMBOL","GENENAME"),
  446. keytype="ENTREZID")
  447. res_3moAWvsSD_df <- left_join(res_3moAWvsSD_df,anno,by="ENTREZID")
  448. #for some reason the annotation put gene ID 11865 to Bmal instead of Arntl so need to replace this gene and make consistent with the old/AD datasets
  449. res_3moAWvsSD_df$SYMBOL[res_3moAWvsSD_df$ENTREZID == "11865"]<-"Arntl"
  450. res_3moAWvsSD_df$GENENAME[res_3moAWvsSD_df$ENTREZID == "11865"]<-"aryl hydrocarbon receptor nuclear translocator-like"
  452. save(res_3moAWvsSD_df, file="allAWvsSD_3mo.Rdata")
  454. threshold <- res_3moAWvsSD_df$padj < padj.cutoff & abs(res_3moAWvsSD_df$log2FoldChange) > log2FoldChange.cutoff
  455. length(which(threshold))
  456. res_3moAWvsSD_df$threshold <- threshold
  457. sig_res_3moAWvsSD <- data.frame(subset(res_3moAWvsSD_df, threshold==TRUE), row.names = 1)
  458. dim(sig_res_3moAWvsSD)
  459. sig_res_3moAWvsSD<-na.omit(sig_res_3moAWvsSD)
  460. dim(sig_res_3moAWvsSD)
  461. write.csv(sig_res_3moAWvsSD, file="sig_res_3moAWvsSD.csv")
  462. save(sig_res_3moAWvsSD, file="sigAWvsSD_3mo.Rdata")
  464. ```
  466. #SD vs SL
  467. ```{r}
  468. dds_3mo_SL_SD<-dds_3mo
  469. design(dds_3mo_SL_SD)<-formula(~Condition)
  470. dds_3mo_SL_SD$Condition<-factor(dds_3mo_SL_SD$Condition, levels=c("AW", "SL", "SD"))
  471. dds_3mo_SL_SD<-DESeq(dds_3mo_SL_SD)
  472. res_3moSLvsSD<-results(dds_3mo_SL_SD, contrast = c("Condition", "SD", "SL"))
  474. res_3moSLvsSD_df <- as.data.frame(res_3moSLvsSD)
  475. res_3moSLvsSD_df<-tibble::rownames_to_column(res_3moSLvsSD_df, "ENTREZID")
  476. res_3moSLvsSD_df$ENTREZID<- as.character(res_3moSLvsSD_df$ENTREZID)
  477. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_3moSLvsSD_df$ENTREZID,
  478. columns=c("SYMBOL","GENENAME"),
  479. keytype="ENTREZID")
  480. res_3moSLvsSD_df <- left_join(res_3moSLvsSD_df,anno,by="ENTREZID")
  481. #for some reason the annotation put gene ID 11865 to Bmal instead of Arntl so need to replace this gene and make consistent with the old/AD datasets
  482. res_3moSLvsSD_df$SYMBOL[res_3moSLvsSD_df$ENTREZID == "11865"]<-"Arntl"
  483. res_3moSLvsSD_df$GENENAME[res_3moSLvsSD_df$ENTREZID == "11865"]<-"aryl hydrocarbon receptor nuclear translocator-like"
  484. save(res_3moSLvsSD_df, file="allSLvsSD_3mo.Rdata")
  486. threshold <- res_3moSLvsSD_df$padj < padj.cutoff & abs(res_3moSLvsSD_df$log2FoldChange) > log2FoldChange.cutoff
  487. length(which(threshold))
  488. res_3moSLvsSD_df$threshold <- threshold
  489. sig_res_3moSLvsSD_df <- data.frame(subset(res_3moSLvsSD_df, threshold==TRUE), row.names = 1)
  490. dim(sig_res_3moSLvsSD_df)
  491. sig_res_3moSLvsSD_df<-na.omit(sig_res_3moSLvsSD_df)
  492. dim(sig_res_3moSLvsSD_df)
  493. save(sig_res_3moSLvsSD_df, file="sigSLvsSD_3mo.Rdata")
  495. ```
  497. #Volcano plots of Aw vs SL or AW vs SD
  498. ```{r volcano plot}
  499. #plot ignoring the condition variable and just comparing time of day effects
  501. EnhancedVolcano(res_3moAWvsSL_df,
  502. lab = res_3moAWvsSL_df$SYMBOL,
  503. x = 'log2FoldChange',
  504. y = 'padj',
  505. title="Aw vs SL",
  506. subtitle = NULL,
  507. xlim = c(-2.5,2.5),
  508. col=c('grey','black', 'orange1' , 'red2'),
  509. pCutoff = 0.05,
  510. FCcutoff = 0.58,
  511. pointSize = 2.0,
  512. labSize = 3,
  513. axisLabSize=10)
  514. #saved as pdf 7x7 "Enhanced_volcano_young_AWvsSL"
  516. EnhancedVolcano(res_3moAWvsSD_df,
  517. lab = res_3moAWvsSD_df$SYMBOL,
  518. x = 'log2FoldChange',
  519. y = 'padj',
  520. title="Aw vs SD",
  521. subtitle = NULL,
  522. xlim = c(-2.5,2.5),
  523. col=c('grey','black', 'orange1' , 'red2'),
  524. pCutoff = 0.05,
  525. FCcutoff = 0.58,
  526. pointSize = 2.0,
  527. labSize = 3,
  528. axisLabSize=10)
  529. #saved as pdf 7x7 "Enhanced_volcano_young_AWvsSD"
  530. ```
  532. ```{r}
  533. # identifying genes unique to either AW vs SL or AW vs SD or overlapping
  534. only_AWvsSL<- subset(sig_res_3moAWvsSL, !SYMBOL %in% sig_res_3moAWvsSD$SYMBOL)
  535. save(only_AWvsSL, file="unique_AWvsSL.Rdata")
  536. only_AWvsSD<- subset(sig_res_3moAWvsSD, !SYMBOL %in% sig_res_3moAWvsSL$SYMBOL)
  537. save(only_AWvsSD, file="unique_AWvsSD.Rdata")
  538. both_AWvsSL_AWvsSD<-subset(sig_res_3moAWvsSL, SYMBOL %in% sig_res_3moAWvsSD$SYMBOL)
  539. save(both_AWvsSL_AWvsSD, file="common_AWvsSLandAWvsSD.Rdata")
  541. ```
  542. #Bar plots of RNAseq data young animals AW vs SD, AW vs SL
  543. ```{r bar plots of SL and SD data}
  544. #plotting the DS and SL changes by circadian genes
  545. #join the fold changes from each results dataframe
  546. #subset the data to get just Fc and padj
  547. sd<-res_3moAWvsSD_df[, c("SYMBOL", "padj", "log2FoldChange")] %>% na.omit()
  548. sd<-sd %>% dplyr::rename(
  549. padj_sd = padj,
  550. log2FC_sd = log2FoldChange)
  552. sl<-res_3moAWvsSL_df[, c("SYMBOL", "padj", "log2FoldChange")] %>% na.omit()
  553. sl<-sl %>% dplyr::rename(
  554. padj_sl = padj,
  555. log2FC_sl = log2FoldChange)
  556. #put the data together
  557. SDandSL_DEGs<-inner_join(sd, sl, by= "SYMBOL")
  559. #Add row that has *s or ns for significance
  560. SDandSL_DEGs_symbol<-SDandSL_DEGs %>%
  561. mutate(p_symbol_sd=case_when(padj_sd<0.001 ~ "***"
  562. , padj_sd<0.01 ~ "**"
  563. , padj_sd<0.05 ~ "*"
  564. , TRUE ~ "ns")) %>%
  565. mutate(p_symbol_sl=case_when(padj_sl<0.001 ~ "***"
  566. , padj_sl<0.01 ~ "**"
  567. , padj_sl<0.05 ~ "*"
  568. , TRUE ~ "ns"))
  569. #just getting log2FC data
  570. SDandSL_DEGs_FC<-SDandSL_DEGs[, c("SYMBOL", "log2FC_sd", "log2FC_sl")] %>% pivot_longer(col=!SYMBOL, names_to = "sleep_group", names_prefix="log2FC_", values_to = "log2FC")
  572. #getting log2FC and padj values
  573. SDandSL_DEGs_FC_padj<-SDandSL_DEGs %>%
  574. pivot_longer(col=!SYMBOL,
  575. names_to = c("measure", "sleep_group"),
  576. names_sep= "_",
  577. values_to = "value")
  579. SDandSL_DEGs_FC_padj %>% pivot_wider(names_from = "measure", values_from = "value") %>% ggplot(aes(x=SYMBOL, y=log2FC, fill=sleep_group)) + geom_bar(stat="identity")
  581. circadian_genes<-c("Nr1d1", "Nr1d2", "Manf", "Fus", "Ciart", "Dbp", "Arntl", "Clock", "Bhlhe41", "Nfil3", "Cry1", "Rorg", "Xbp1", "Per3")
  583. SDandSL_DEGs_FC_padj[SDandSL_DEGs_FC_padj$SYMBOL %in% circadian_genes, ] %>%
  584. pivot_wider(names_from = "measure", values_from = "value") %>%
  585. ggplot(aes(x=reorder(SYMBOL,-log2FC), y=log2FC, fill=sleep_group)) +
  586. geom_bar(stat="identity", position="dodge", color="black") +
  587. #geom_text(aes(label=sprintf(fmt="%0.2e", round(padj, digits = 2)))) +
  588. scale_fill_manual(values = alpha(c("white", "royalblue4"), 0.8)) +
  589. theme_classic()
  590. #saved as "barplot_SLandSD_circgenes.pdf" 6.5x4
  593. ```
  596. #Wild-type aged mouse samples, day vs night
  597. ```{r}
  598. #make a dds file from galaxy generated normalized counts (rlog)
  600. dds_12moWT<-read.csv("12mo_counts.csv", row.names = 1)
  601. View(dds_12moWT)
  602. #but this is all the samples
  603. #need to get only the WT
  604. metadata<-read.csv("12mo_metadata.csv", row.names = 1)
  606. metadata$Genotype<-factor(metadata$Genotype, levels=c("WT", "AD"))
  607. metadata_WT <- dplyr::filter(metadata, Genotype == "WT")
  608. metadata_WT$TimeOfDay<-factor(metadata_WT$TimeOfDay, levels=c("Night", "Day"))
  609. View(metadata_WT)
  611. #explore
  612. rownames(metadata_WT)
  613. colnames(dds_12moWT)
  614. #filter dds file by metadata rownames
  615. dds_12moWT<-dds_12moWT[, rownames(metadata_WT)]
  616. #make sure the rows and columns match
  617. all(rownames(metadata_WT) == colnames(dds_12moWT))
  619. ```
  621. #Normalizing Old Wt samples day vs night
  622. ```{r}
  623. #creating the DESeq2 object, using *TimeOfDay as condition
  624. dds_12moWT<-DESeqDataSetFromMatrix(countData = dds_12moWT, colData = metadata_WT, design = ~ TimeOfDay)
  625. #determine size factors to use for normalization
  626. dds_12moWT <- estimateSizeFactors(dds_12moWT)
  628. sizeFactors(dds_12moWT)
  629. #extract the normalized counts
  630. normalized_counts_12moWT <-counts(dds_12moWT, normalized = TRUE)
  632. View(normalized_counts_12moWT)
  634. write.csv(normalized_counts_12moWT, file="normalized_counts_12moWT.csv")
  636. ###Getting the ENTREZID as rownames for the damn dataframe.
  637. normalized_counts_12moWT_df <- data.frame(ENTREZID=row.names(normalized_counts_12moWT), normalized_counts_12moWT)
  639. #perform unsupervised clustering analysis:log transformation
  640. nrow(dds_12moWT)
  641. vsd_12moWT <- vst(dds_12moWT, blind = TRUE)
  642. vsd_mat_12moWT <- assay(vsd_12moWT)
  643. vsd_cor_12moWT <- cor(vsd_mat_12moWT)
  644. View(vsd_cor_12moWT)
  645. ```
  646. #DEG of wild-type day vs night
  647. ```{r}
  649. dds_12moWT<-DESeq(dds_12moWT)
  650. res_12moWT<-results(dds_12moWT, contrast = c("TimeOfDay","Night", "Day"))
  652. res_12moWT_df <- as.data.frame(res_12moWT)
  653. res_12moWT_df<-tibble::rownames_to_column(res_12moWT_df, "ENTREZID")
  654. res_12moWT_df$ENTREZID<- as.character(res_12moWT_df$ENTREZID)
  655. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_12moWT_df$ENTREZID,
  656. columns=c("SYMBOL","GENENAME"),
  657. keytype="ENTREZID")
  658. res_12moWT_df <- left_join(res_12moWT_df,anno,by="ENTREZID")
  660. threshold <- res_12moWT_df$padj < padj.cutoff & abs(res_12moWT_df$log2FoldChange) > log2FoldChange.cutoff
  661. length(which(threshold))
  662. res_12moWT_df$threshold <- threshold
  663. sig_res_12moWT_df <- data.frame(subset(res_12moWT_df, threshold==TRUE), row.names = 1)
  664. dim(sig_res_12moWT_df)
  665. sig_res_12moWT_df<-na.omit(sig_res_12moWT_df)
  666. dim(sig_res_12moWT_df)
  667. #save files
  668. save(res_12moWT_df, file="res_12moWT_df.Rdata")
  669. write.csv(res_12moWT_df, file= "res_12moWT_df.csv")
  670. save(sig_res_12moWT_df, file="sig_12moWT_df.Rdata")
  671. write.csv(sig_res_12moWT_df, file= "sig_res_12moWT_df.csv")
  672. ```
  674. #AD samples
  675. ```{r}
  676. #make a dds file from galaxy generated normalized counts (rlog)
  678. dds_12moAD<-read.csv("12mo_counts.csv", row.names = 1)
  679. View(dds_12moAD)
  680. #but this is all the samples
  681. #need to get only the AD
  682. metadata<-read.csv("12mo_metadata.csv", row.names = 1)
  684. metadata$Genotype<-factor(metadata$Genotype, levels=c("WT", "AD"))
  685. metadata_AD <- dplyr::filter(metadata, Genotype == "AD")
  686. metadata_AD$TimeOfDay<-factor(metadata_AD$TimeOfDay, levels=c("Night", "Day"))
  687. View(metadata_AD)
  689. #explore
  690. rownames(metadata_AD)
  691. colnames(dds_12moAD)
  692. #filter dds file by metadata ronames
  693. dds_12moAD<-dds_12moAD[, rownames(metadata_AD)]
  694. #make sure the rows and columns match
  695. all(rownames(metadata_AD) == colnames(dds_12moAD))
  697. ```
  699. #Normalizing Old AD samples day vs night
  700. ```{r}
  701. #creating the DESeq2 object, using *TimeOfDay as condition
  702. dds_12moAD<-DESeqDataSetFromMatrix(countData = dds_12moAD, colData = metadata_AD, design = ~ TimeOfDay)
  703. #determine size factors to use for normalization
  704. dds_12moAD <- estimateSizeFactors(dds_12moAD)
  706. sizeFactors(dds_12moAD)
  707. #extract the normalized counts
  708. normalized_counts_12moAD <-counts(dds_12moAD, normalized = TRUE)
  710. View(normalized_counts_12moAD)
  712. write.csv(normalized_counts_12moAD, file="normalized_counts_12moAD.csv")
  714. ###Getting the ENTREZID as rownames for the damn dataframe.
  715. normalized_counts_12moAD_df <- data.frame(ENTREZID=row.names(normalized_counts_12moAD), normalized_counts_12moAD)
  717. #perform unsupervised clustering analysis:log transformation
  718. nrow(dds_12moAD)
  719. vsd_12moAD <- vst(dds_12moAD, blind = TRUE)
  720. vsd_mat_12moAD <- assay(vsd_12moAD)
  721. vsd_cor_12moAD <- cor(vsd_mat_12moAD)
  722. View(vsd_cor_12moAD)
  723. ```
  725. #DEG of AD day vs night
  726. ```{r}
  728. dds_12moAD<-DESeq(dds_12moAD)
  729. res_12moAD<-results(dds_12moAD, contrast = c("TimeOfDay", "Night", "Day"))
  731. res_12moAD_df <- as.data.frame(res_12moAD)
  732. res_12moAD_df<-tibble::rownames_to_column(res_12moAD_df, "ENTREZID")
  733. res_12moAD_df$ENTREZID<- as.character(res_12moAD_df$ENTREZID)
  734. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_12moAD_df$ENTREZID,
  735. columns=c("SYMBOL","GENENAME"),
  736. keytype="ENTREZID")
  737. res_12moAD_df <- left_join(res_12moAD_df,anno,by="ENTREZID")
  739. threshold <- res_12moAD_df$padj < padj.cutoff & abs(res_12moAD_df$log2FoldChange) > log2FoldChange.cutoff
  740. length(which(threshold))
  741. res_12moAD_df$threshold <- threshold
  742. sig_res_12moAD_df <- data.frame(subset(res_12moAD_df, threshold==TRUE), row.names = 1)
  743. dim(sig_res_12moAD_df)
  744. sig_res_12moAD_df<-na.omit(sig_res_12moAD_df)
  745. dim(sig_res_12moAD_df)
  746. save(res_12moAD_df, file="res_12moAD_df.Rdata")
  747. write.csv(res_12moAD_df, file="res_12moAD_df.csv")
  748. save(sig_res_12moAD_df, file="sig_12moAD_df.Rdata")
  749. write.csv(sig_res_12moAD_df, file= "sig_res_12moAD_df.csv")
  750. ```
  752. #Including the comparison of the WT vs AD group as this is an obvious question even though it has been done already (daMesquita et al)
  753. ```{r }
  754. dds_12mo<-read.csv("12mo_counts.csv", row.names = 1)
  755. View(dds_12mo)
  756. metadata_WTAD<-read.csv("12mo_metadata.csv", row.names = 1)
  758. metadata_WTAD$Genotype<-factor(metadata_WTAD$Genotype, levels=c("WT", "AD"))
  760. metadata_WTAD$TimeOfDay<-factor(metadata_WTAD$TimeOfDay, levels=c("Night", "Day"))
  761. View(metadata_WTAD)
  763. #explore
  764. rownames(metadata_WTAD)
  765. colnames(dds_12mo)
  766. #filter dds file by metadata ronames
  767. dds_12mo<-dds_12mo[, rownames(metadata_WTAD)]
  768. #make sure the rows and columns match
  769. all(rownames(metadata_WTAD) == colnames(dds_12mo))
  771. ```
  773. #Normalizing all samples together
  774. ```{r }
  775. #creating the DESeq2 object, using *TimeOfDay and Genotype as condition
  776. dds_12mo<-DESeqDataSetFromMatrix(countData = dds_12mo, colData = metadata_WTAD, design = ~ TimeOfDay + Genotype)
  777. #determine size factors to use for normalization
  778. dds_12mo <- estimateSizeFactors(dds_12mo)
  780. sizeFactors(dds_12mo)
  781. #extract the normalized counts
  782. normalized_counts_12mo <-counts(dds_12mo, normalized = TRUE)
  784. View(normalized_counts_12mo)
  786. write.csv(normalized_counts_12mo, file="normalized_counts_12mo_all.csv")
  788. ###Getting the ENTREZID as rownames for the damn dataframe.
  789. normalized_counts_12mo_df <- data.frame(ENTREZID=row.names(normalized_counts_12mo), normalized_counts_12mo)
  791. #perform unsupervised clustering analysis:log transformation
  792. nrow(dds_12mo)
  793. vsd_12mo <- vst(dds_12mo, blind = TRUE)
  794. vsd_mat_12mo <- assay(vsd_12mo)
  795. vsd_cor_12mo <- cor(vsd_mat_12mo)
  796. View(vsd_cor_12mo)
  798. ```
  799. #Plots of variation
  800. ```{r pressure, echo=FALSE}
  801. pheatmap(vsd_cor_12mo, cluster_rows = TRUE, cluster_cols = TRUE)
  802. pheatmap(vsd_cor_12mo, annotation = dplyr::select(metadata_WTAD, c(TimeOfDay, Genotype)))
  803. plotPCA(vsd_3mo, intgroup = "Genotype")
  804. plotPCA(vsd_3mo, intgroup=c("Genotype", "TimeOfDay"))
  805. ```
  807. #Degs for all old samples night vs day correcting for genotype, then genotype correcting for night vs day
  808. ```{r}
  809. #time of day accounting for genotype
  810. design(dds_12mo)<-formula(~ Genotype + TimeOfDay)
  811. dds_12mo<-DESeq(dds_12mo)
  812. res_12moToD<-results(dds_12mo, contrast = c("TimeOfDay", "Night", "Day"))
  814. res_12moToD_df <- as.data.frame(res_12moToD)
  815. res_12moToD_df<-tibble::rownames_to_column(res_12moToD_df, "ENTREZID")
  816. res_12moToD_df$ENTREZID<- as.character(res_12moToD_df$ENTREZID)
  817. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_12moToD_df$ENTREZID,
  818. columns=c("SYMBOL","GENENAME"),
  819. keytype="ENTREZID")
  820. res_12moToD_df <- left_join(res_12moToD_df,anno,by="ENTREZID")
  822. threshold <- res_12moToD_df$padj < padj.cutoff & abs(res_12moToD_df$log2FoldChange) > log2FoldChange.cutoff
  823. length(which(threshold))
  824. res_12moToD_df$threshold <- threshold
  825. sig_res_12moToD_df <- data.frame(subset(res_12moToD_df, threshold==TRUE), row.names = 1)
  826. dim(sig_res_12moToD_df)
  827. sig_res_12moToD_df<-na.omit(sig_res_12moToD_df)
  828. dim(sig_res_12moToD_df)
  829. save(sig_res_12moToD_df, file="sig_12moToD_df.Rdata")
  830. write.csv(sig_res_12moToD_df, file="sig_res_12moToD_df.csv")
  832. #Genotype accounting for time of day
  833. design(dds_12mo)<-formula(~ TimeOfDay + Genotype)
  834. dds_12mo<-DESeq(dds_12mo)
  835. res_12moGT<-results(dds_12mo, contrast = c("Genotype", "AD", "WT"))
  837. res_12moGT_df <- as.data.frame(res_12moGT)
  838. res_12moGT_df<-tibble::rownames_to_column(res_12moGT_df, "ENTREZID")
  839. res_12moGT_df$ENTREZID<- as.character(res_12moGT_df$ENTREZID)
  840. anno <- AnnotationDbi::select(org.Mm.eg.db,keys=res_12moGT_df$ENTREZID,
  841. columns=c("SYMBOL","GENENAME"),
  842. keytype="ENTREZID")
  843. res_12moGT_df <- left_join(res_12moGT_df,anno,by="ENTREZID")
  845. threshold <- res_12moGT_df$padj < padj.cutoff & abs(res_12moGT_df$log2FoldChange) > log2FoldChange.cutoff
  846. length(which(threshold))
  847. res_12moGT_df$threshold <- threshold
  848. sig_res_12moGT_df <- data.frame(subset(res_12moGT_df, threshold==TRUE), row.names = 1)
  849. dim(sig_res_12moGT_df)
  850. sig_res_12moGT_df<-na.omit(sig_res_12moGT_df)
  851. dim(sig_res_12moGT_df)
  852. save(sig_res_12moGT_df, file="sig_12moGT_df.Rdata")
  853. write.csv(sig_res_12moGT_df, file="sig_res_12moGT_df.csv")
  855. ```
  857. #Volcano plots of old WT or AD samples
  858. ```{r volcano plot}
  859. #plot comparing time of day effects
  861. EnhancedVolcano(res_12moWT_df,
  862. lab = res_12moWT_df$SYMBOL,
  863. x = 'log2FoldChange',
  864. y = 'padj',
  865. title="Night vs Day",
  866. subtitle = NULL,
  867. xlim = c(-2.5,2.5),
  868. col=c('grey','black', 'orange1' , 'red2'),
  869. pCutoff = 0.05,
  870. FCcutoff = 0.58,
  871. pointSize = 2.0,
  872. labSize = 5,
  873. axisLabSize=10)
  874. #saved as pdf 7x7 "Enhanced_volcano_Old_WT_NightvsDay.pdf"
  876. EnhancedVolcano(res_3moAWvsSD_df,
  877. lab = res_3moAWvsSD_df$SYMBOL,
  878. x = 'log2FoldChange',
  879. y = 'padj',
  880. title="Aw vs SD",
  881. subtitle = NULL,
  882. xlim = c(-2.5,2.5),
  883. col=c('grey','black', 'orange1' , 'red2'),
  884. pCutoff = 0.05,
  885. FCcutoff = 0.58,
  886. pointSize = 2.0,
  887. labSize = 3,
  888. axisLabSize=10)
  889. #saved as pdf 7x7 "Enhanced_volcano_young_AWvsSD"
  890. ```
  892. #Plotting the normalized counts from the DEG lists of young, old and AD groups
  893. #using 3mo data that has 3 condition metadata
  894. ```{r heatmaps of normalized counts}
  895. # normalized_counts_3mo_cond_df
  896. #have to annotate each dataset to get gene SYMBOLS
  897. class(normalized_counts_3mo_cond_df$ENTREZID) #character
  898. normalized_counts_3mo_cond_df$ENTREZID<-as.character(normalized_counts_3mo_cond_df$ENTREZID)
  899. keys(org.Mm.eg.db, keytype="ENTREZID")[1:10]
  900. anno <- AnnotationDbi::select(org.Mm.eg.db,
  901. keys=normalized_counts_3mo_cond_df$ENTREZID,
  902. columns=c("SYMBOL","GENENAME"),
  903. keytype="ENTREZID")
  905. head(anno)
  906. dim(anno)
  907. dim(normalized_counts_3mo_cond_df)
  908. #they match!
  910. normalized_counts_3mo_cond_df_annotated <- left_join(normalized_counts_3mo_cond_df,anno,by="ENTREZID")
  911. head(normalized_counts_3mo_cond_df_annotated)
  913. #Extract the normalized count data for the DEGs of the AWvsSL comparison
  914. ##For AWvsSL
  915. col_order<-c("SYMBOL", "M1", "M2", "M3", "M4", "M9", "M10", "M11")
  916. norm_SYMBOL_AWvsSL<-normalized_counts_3mo_cond_df_annotated[ ,col_order]
  917. norm_SYMBOL_AWvsSL<-norm_SYMBOL_AWvsSL %>% drop_na(SYMBOL)
  919. norm_SYMBOL_AWvsSL<-data.frame(norm_SYMBOL_AWvsSL, row.names = 1)
  921. norm_sigAWvsSL <- norm_SYMBOL_AWvsSL[sig_res_3moAWvsSL$SYMBOL, ]
  922. norm_sigAWvsSL<-na.omit(norm_sigAWvsSL)
  923. dim(norm_sigAWvsSL)
  925. #need this for the complex heatmap
  926. norm_sigAWvsSL_mat<-as.matrix(norm_sigAWvsSL)
  928. #For AWvsSD
  929. col_order2<-c("SYMBOL", "M5", "M6", "M7", "M8", "M9", "M10", "M11")
  930. norm_SYMBOL_AWvsSD<-normalized_counts_3mo_cond_df_annotated[ ,col_order2]
  931. norm_SYMBOL_AWvsSD<-norm_SYMBOL_AWvsSD %>% drop_na(SYMBOL)
  932. norm_SYMBOL_AWvsSD<-data.frame(norm_SYMBOL_AWvsSD, row.names = 1)
  933. norm_sigAWvsSD <- norm_SYMBOL_AWvsSD[sig_res_3moAWvsSD$SYMBOL, ]
  934. norm_sigAWvsSD<-na.omit(norm_sigAWvsSD)
  935. dim(norm_sigAWvsSD)
  936. norm_sigAWvsSD_mat<-as.matrix(norm_sigAWvsSD)
  938. #Annotate our heatmap (optional)
  939. #annotation_SL <- data.frame(sampletype=metadata_SL[, "Condition"],
  940. #row.names=rownames(metadata_SL))
  942. # normalized_counts_12moWT_df
  943. class(normalized_counts_12moWT_df$ENTREZID) #character
  944. normalized_counts_12moWT_df$ENTREZID<-as.character(normalized_counts_12moWT_df$ENTREZID)
  946. anno <- AnnotationDbi::select(org.Mm.eg.db,
  947. keys=normalized_counts_12moWT_df$ENTREZID,
  948. columns=c("SYMBOL","GENENAME"),
  949. keytype="ENTREZID")
  951. head(anno)
  952. dim(anno)
  953. dim(normalized_counts_12moWT_df)
  954. #they match!
  956. # normalized_counts_12moWT_df
  957. normalized_counts_12moWT_df_annotated <- left_join(normalized_counts_12moWT_df,anno,by="ENTREZID")
  958. head(normalized_counts_12moWT_df_annotated)
  960. col_order_old<-c("SYMBOL", "M13", "M14", "M15", "M16", "M21", "M22", "M23", "M24")
  961. norm_SYMBOL_12moWT<-normalized_counts_12moWT_df_annotated[ ,col_order_old]
  962. norm_SYMBOL_12moWT<-norm_SYMBOL_12moWT %>% drop_na(SYMBOL)
  963. norm_SYMBOL_12moWT<-data.frame(norm_SYMBOL_12moWT, row.names = 1)
  965. #to get only the sig genes from the 3mo AWvsSL list
  966. norm_sig12moWT <- norm_SYMBOL_12moWT[sig_res_3moAWvsSL$SYMBOL, ]
  967. norm_sig12moWT<-na.omit(norm_sig12moWT)
  968. dim(norm_sig12moWT)
  969. norm_sig12moWT_mat<-as.matrix(norm_sig12moWT)
  971. #############################
  972. # normalized_counts_12moAD_df
  973. class(normalized_counts_12moAD_df$ENTREZID) #character
  974. normalized_counts_12moAD_df$ENTREZID<-as.character(normalized_counts_12moAD_df$ENTREZID)
  976. anno <- AnnotationDbi::select(org.Mm.eg.db,
  977. keys=normalized_counts_12moAD_df$ENTREZID,
  978. columns=c("SYMBOL","GENENAME"),
  979. keytype="ENTREZID")
  981. head(anno)
  982. dim(anno)
  983. dim(normalized_counts_12moAD_df)
  984. #they match!
  986. # normalized_counts_12moAD_df
  987. normalized_counts_12moAD_df_annotated <- left_join(normalized_counts_12moAD_df,anno,by="ENTREZID")
  988. head(normalized_counts_12moAD_df_annotated)
  990. col_order_AD<-c("SYMBOL", "M17", "M18", "M19", "M20", "M25", "M26", "M27", "M28")
  991. norm_SYMBOL_12moAD<-normalized_counts_12moAD_df_annotated[ ,col_order_AD]
  992. norm_SYMBOL_12moAD<-norm_SYMBOL_12moAD %>% drop_na(SYMBOL)
  993. norm_SYMBOL_12moAD<-data.frame(norm_SYMBOL_12moAD, row.names = 1)
  995. norm_sig12moAD <- norm_SYMBOL_12moAD[sig_res_3moAWvsSL$SYMBOL, ]
  996. norm_sig12moAD<-na.omit(norm_sig12moAD)
  997. dim(norm_sig12moAD)
  998. norm_sig12moAD_mat<-as.matrix(norm_sig12moAD)
  999. ```
  1000. #Heatmaps of normalized counts
  1001. ```{r heatmaps of normalized counts for sig genes}
  1002. #Set a color palette
  1003. heat.colors <- brewer.pal(11, "RdBu")
  1004. norm_sigAWvsSL_mat_scaled<-t(scale((t(norm_sigAWvsSL_mat))))
  1005. #sig genes for AWvsSL
  1006. hm<-pheatmap(norm_sigAWvsSL_mat, color = heat.colors, cluster_rows = T, show_rownames=T,
  1007. border_color=NA, fontsize = 10, scale="row",
  1008. fontsize_row = 3, height=20, main="AWvsSL sig genes")
  1010. hm2<-Heatmap(norm_sigAWvsSL_mat_scaled, col=heat.colors, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1013. draw(hm2, heatmap_legend_side="left")
  1014. #save image as "pheatmap_AWvsSLnormcounts_AWvsSLgenelist.pdf"
  1016. #to get row order for subsequent heatmaps
  1017. ro1<-row_order(hm2)
  1019. # to manually set which rows are labelled
  1020. ro_label<-rownames(norm_sigAWvsSL_mat_scaled)[row_order(hm2)]
  1021. label_file<-read.csv("3mo_heatmap_label.csv")
  1022. rownumbers<-label_file$row_number
  1023. rowgenes<-label_file$SYMBOL
  1024. ha <-rowAnnotation(foo=anno_mark(at = rownumbers, labels = rowgenes, labels_gp=gpar(fontsize=6)))
  1026. hm2<-Heatmap(norm_sigAWvsSL_mat_scaled, col=heat.colors, show_column_dend = FALSE, show_row_dend = FALSE, show_column_names = FALSE, right_annotation=ha, show_row_names = FALSE, row_names_gp=gpar(fontsize=5), clustering_distance_rows = "euclidean")
  1027. draw(hm2, heatmap_legend_side="left")
  1028. #saved as pdf "pheatmap_AWvsSLnormcounts_AWvsSLgenelist_select_labels_2.pdf" 14x4
  1030. #sig genes for AwvsSD
  1031. norm_sigAWvsSD_mat_scaled<-t(scale((t(norm_sigAWvsSD_mat))))
  1032. hm3<-Heatmap(norm_sigAWvsSD_mat_scaled, col=heat.colors, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1033. #manually set which genes are labelled as above
  1034. ro_label<-rownames(norm_sigAWvsSD_mat_scaled)[row_order(hm3)]
  1035. label_file<-read.csv("3moAWvsSD_heatmap_label.csv")
  1036. rownumbers<-label_file$row_number
  1037. rowgenes<-label_file$SYMBOL
  1038. ha <-rowAnnotation(foo=anno_mark(at = rownumbers, labels = rowgenes, labels_gp=gpar(fontsize=10)))
  1039. hm3<-Heatmap(norm_sigAWvsSD_mat_scaled, col=heat.colors, show_column_dend = FALSE, show_row_dend = FALSE, show_column_names = FALSE, right_annotation=ha, show_row_names = FALSE, row_names_gp=gpar(fontsize=5), clustering_distance_rows = "euclidean")
  1040. draw(hm3, heatmap_legend_side="left")
  1041. #save image as "pheatmap_AWvsSDnormcounts_AWvsSDgenelist.pdf" 4x4 portrait
  1043. #sig genes from AWvsSL in old WT norm counts
  1044. norm_sig12moWT_mat_scaled<-t(scale((t(norm_sig12moWT_mat))))
  1045. hm4<-Heatmap(norm_sig12moWT_mat_scaled, col=heat.colors, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1046. draw(hm4, heatmap_legend_side="left")
  1047. hm4<-Heatmap(norm_sig12moWT_mat_scaled, col=heat.colors, row_order=ro1, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1048. draw(hm4, heatmap_legend_side="left")
  1049. #save image as "pheatmap_oldWTnormcounts_AWvsSLgenelist_order.pdf" 14x4 portrait or pheatmap_oldWTnormcounts_AWvsSLgenelist.pdf not in order
  1051. #heatmap of top 500 genes of the WT old data
  1052. #need to first get the sig genes in order
  1053. sig_res_12moWT_df_ordered <- sig_res_12moWT_df[order(sig_res_12moWT_df$padj), ]
  1054. top500_sig_res_12moWT <- sig_res_12moWT_df_ordered[1:500, ]
  1055. top500_sig_res_12moWT_genes<-norm_SYMBOL_12moWT[top500_sig_res_12moWT$SYMBOL, ]
  1056. top500_sig_res_12moWT_genes<-na.omit(top500_sig_res_12moWT_genes)
  1057. dim(top500_sig_res_12moWT_genes)
  1058. #make into matrix and scale
  1059. top500_sig_res_12moWT_genes_mat<-as.matrix(top500_sig_res_12moWT_genes)
  1060. top500_sig_res_12moWT_genes_mat_scaled<-t(scale((t(top500_sig_res_12moWT_genes_mat))))
  1062. #now make the heatmap!
  1063. hm4_1<-Heatmap(top500_sig_res_12moWT_genes_mat_scaled, col=heat.colors, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1064. draw(hm4_1, heatmap_legend_side="left")
  1065. #but to select only some genes to be labelled need to save to csv then do in excel and get row #s too
  1066. write.csv(top500_sig_res_12moWT_genes, file = "top500_sig_res_12moWT_genes.csv")
  1067. #save file with log2FC as well
  1068. write.csv(top500_sig_res_12moWT, file="top500_sig_res_12moWT.csv")
  1069. label_file<-read.csv("12moWT_heatmap_label.csv")
  1070. rownumbers<-label_file$row_number
  1071. rowgenes<-label_file$SYMBOL
  1072. ha <-rowAnnotation(foo=anno_mark(at = rownumbers, labels = rowgenes, labels_gp=gpar(fontsize=10)))
  1073. hm4_1<-Heatmap(top500_sig_res_12moWT_genes_mat_scaled, col=heat.colors, show_column_dend = FALSE, show_row_dend = FALSE, right_annotation=ha, show_row_names = FALSE, row_names_gp=gpar(fontsize=5), clustering_distance_rows = "euclidean")
  1074. draw(hm4_1, heatmap_legend_side="left")
  1076. #sig genes from AWvsSL in old AD norm counts
  1077. norm_sig12moAD_mat_scaled<-t(scale((t(norm_sig12moAD_mat))))
  1078. hm5<-Heatmap(norm_sig12moAD_mat_scaled, col=heat.colors, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1079. draw(hm5, heatmap_legend_side="left")
  1080. hm5<-Heatmap(norm_sig12moAD_mat_scaled, col=heat.colors, row_order=ro1, row_names_gp = gpar(fontsize = 4), clustering_distance_rows = "euclidean")
  1081. draw(hm5, heatmap_legend_side="left")
  1082. #save image as "pheatmap_oldADnormcounts_AWvsSLgenelist_order.pdf" 14x4 portrait or pheatmap_oldWTnormcounts_AWvsSLgenelist.pdf not in order
  1085. ```
  1088. #Genes of interest across young, aged, AD groups
  1089. ```{r plotting genes of interest by fuctional groups}
  1090. SLC_family_genes<-c("Slc26a2", "Slc4a10", "Slc4a5", "Slc12a2", "Slc25a28", "Slc12a7", "Slc12a4", "Slc4a2", "Slc38a3", "Slc22a8")
  1091. #"Slc12a1", "Slc12a6", "Slc25a37", "Slc5a5", "Slco4a1" not sig in any
  1092. all_SLC_genes<-wide_young_old_ad %>% dplyr::filter(grepl("Slc", SYMBOL))
  1093. all_SLC_genes_list<-all_SLC_genes$SYMBOL
  1094. Na_K_ATPase_family_genes<-c("Atp1a2", "Atp1a1", "Atp1b1")
  1095. # "Atp1b1", "Atp1b2", "Atp1a4" not sig in any
  1097. Carbonic_anhydrase_family_genes<-c("Ca4", "Ca2", "Ca3","Ca13", "Ca8", "Ca14")
  1098. #none of these are significant in the groups
  1100. ATP_synthase_family<-c("Atp5l", "Atp23", "Atp5g3", "Atpaf1", "Atp5b" ) #none are significant
  1102. tight_junctions<-c("Cldn5", "Cldn2", "Cldn3", "Cldn19", "Ocln")
  1103. #"Cldn18" not expressed, "Cldn11", "Cldn1" not sig in any
  1105. amyloid_proteins<-c("Abpa3", "Apbb1ip", "App") #only APP is different
  1107. circadian_genes<-c("Nr1d1", "Nr1d2", "Manf", "Fus", "Ciart", "Dbp", "Arntl", "Clock", "Bhlhe41", "Nfil3", "Cry1", "Rorg", "Xbp1", "Per3")
  1108. # "Sirt1" not sig in any
  1110. #DEGs from Dani et al, age-dependent
  1111. age_genes<-read.csv("age_genes_dani_et_al.csv")
  1112. age_gene_names<-age_genes$Age_dependent_genes
  1113. ```
  1115. ```{r plots from above}
  1116. res_3moAWvsSL_df[res_3moAWvsSL_df$SYMBOL %in% SLC_family_genes, ] %>% dplyr::filter(padj<0.05) %>%
  1117. ggplot(aes(x=SYMBOL, y=log2FoldChange)) +
  1118. geom_bar(stat="identity") +
  1119. geom_text(aes(label= sprintf("%.2f",round(padj, digits=2))))
  1120. ```
  1123. ```{r}
  1124. #combine young AWvsSL, oldWT and old AD datasets with FC and pvalues to have on same graph
  1125. #need to take out the NA gene symbols from each dataset to be able to pivot wide later, otherwise numbers are converted to list and then nothing works
  1127. young_fc_padj<-res_3moAWvsSL_df[, c("SYMBOL", "log2FoldChange", "padj")]
  1128. #add column to designate the dataset group
  1129. young_fc_padj<-young_fc_padj %>% mutate(Group= "Young")
  1130. young_fc_padj<-young_fc_padj %>% drop_na(SYMBOL)
  1132. oldWT_fc_padj<-res_12moWT_df[, c("SYMBOL", "log2FoldChange", "padj") ]
  1133. oldWT_fc_padj<-oldWT_fc_padj %>% mutate(Group= "OldWT")
  1134. oldWT_fc_padj<-oldWT_fc_padj%>% drop_na(SYMBOL)
  1136. oldAD_fc_padj<-res_12moAD_df[, c("SYMBOL", "log2FoldChange", "padj") ]
  1137. oldAD_fc_padj<-oldAD_fc_padj %>% mutate(Group= "OldAD")
  1138. oldAD_fc_padj<-oldAD_fc_padj%>% drop_na(SYMBOL)
  1140. young_old_ad<-dplyr::bind_rows(young_fc_padj,oldWT_fc_padj)
  1141. young_old_ad<-dplyr::bind_rows(young_old_ad, oldAD_fc_padj)
  1143. #Add row that has *s or ns for significance
  1144. young_old_ad_symbol<-young_old_ad %>%
  1145. mutate(padj_symbol=case_when(padj<0.001 ~ "***"
  1146. , padj<0.01 ~ "**"
  1147. , padj<0.05 ~ "*"
  1148. , TRUE ~ NA))
  1149. write.csv(young_old_ad_symbol, file="young_old_ad_long_FC_padj.csv")
  1150. #to place p-values in the middle of the bars, need to calculate half the fold change
  1151. young_old_ad_symbol<-young_old_ad_symbol %>%
  1152. mutate(half_fc=log2FoldChange*0.5)
  1154. #need to filter the data so that we have only genes left that are significant in at least one group
  1155. #need to do this by gene
  1157. wide_young_old_ad<-pivot_wider(young_old_ad, names_from = Group, values_from = c(log2FoldChange, padj))
  1161. wide_young_old_ad<-wide_young_old_ad %>%
  1162. mutate(Keep = case_when(padj_Young <0.05 ~ "SIG",
  1163. padj_OldWT <0.05 ~ "SIG",
  1164. padj_OldAD <0.05 ~ "SIG",
  1165. .default = "Not_SIG") )
  1167. #now filter the long data by values to keep from wide data
  1168. wide_mod<-wide_young_old_ad[wide_young_old_ad$Keep =="SIG",]
  1169. #make list with only sig genes
  1170. sig_list<-wide_mod$SYMBOL
  1171. #filter the whole dataset by the sig list
  1172. young_old_ad_filtered_bysig<-young_old_ad_symbol[young_old_ad_symbol$SYMBOL %in% sig_list, ]
  1173. young_old_ad_filtered_bysig$Group<-factor(young_old_ad_filtered_bysig$Group, levels=c("Young", "OldWT", "OldAD"))
  1174. #want a list of ones that are sig in the young group
  1175. young_sig_genes<-young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$Group=="Young", ] %>% dplyr::filter(padj<0.05) %>% dplyr::filter(log2FoldChange>0.58 | log2FoldChange < -0.58)
  1176. young_sig_genes_list<-young_sig_genes$SYMBOL
  1178. young_old_conserved<-young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% young_sig_genes_list, ] %>% subset(Group=="OldWT") %>% dplyr::filter(padj<0.05) %>% dplyr::filter(log2FoldChange>0.58 | log2FoldChange < -0.58)
  1179. young_old_conserved_list<-young_old_conserved$SYMBOL
  1181. write.csv(young_old_conserved, file= "young_old_conserved.csv")
  1182. #make matrix with young/old data to plot as heatmap
  1183. young_old_wide<-wide_young_old_ad[wide_young_old_ad$SYMBOL %in% young_old_conserved_list, c(1:4)]
  1184. young_old_mat<-as.matrix(young_old_wide[,-1])
  1185. rownames(young_old_mat)<-young_old_wide$SYMBOL
  1186. write.csv(wide_young_old_ad, file = "wide_young_aged_ad_FC_padj.csv")
  1188. ```
  1190. ```{r plotting gene expression across groups}
  1192. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% circadian_genes, ] %>%
  1193. ggplot(aes(fill=Group, x=SYMBOL, y=log2FoldChange)) +
  1194. geom_bar(position="dodge", stat="identity") +
  1195. geom_text(aes(label= sprintf("%.2f",round(padj, digits=2))), size=3,vjust=1.5, position=position_dodge(0.9)) +
  1196. theme_classic() +
  1197. scale_fill_viridis_d(option="inferno", begin=0.2, end=0.8)
  1199. #label the p values as symbols instead
  1200. #circadian related genes
  1201. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% circadian_genes, ] %>%
  1202. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1203. geom_bar(position="dodge", stat="identity", color="black") +
  1204. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1205. theme_classic() +
  1206. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1207. coord_flip()
  1209. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% circadian_genes, ] %>%
  1210. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1211. geom_bar(position="dodge", stat="identity", color="black") +
  1212. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1213. theme_classic() +
  1214. theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  1215. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7))
  1216. #saved as "barplot_young_old_ad_circadian_family_genes_rotated.pdf"
  1218. #plotting the conserved genes in aged animals
  1219. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% young_old_conserved_list, ] %>%
  1220. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1221. geom_bar(position="dodge", stat="identity", color="black") +
  1222. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1223. theme_classic() +
  1224. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1225. coord_flip()
  1226. # cool but can't see anything. try another heatmap
  1228. pheatmap(young_old_mat)
  1230. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% SLC_family_genes, ] %>%
  1231. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1232. geom_bar(position="dodge", stat="identity", color="black") +
  1233. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1234. theme_classic() +
  1235. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1236. coord_flip()
  1237. #saved as "barplot_young_old_ad_slc_family_genes.pdf" 5x4
  1239. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% all_SLC_genes_list, ] %>%
  1240. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1241. geom_bar(position="dodge", stat="identity", color="black") +
  1242. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1243. theme_classic() +
  1244. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1245. coord_flip()
  1247. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% Na_K_ATPase_family_genes, ] %>%
  1248. ggplot(aes(fill=Group, x=SYMBOL, y=log2FoldChange)) +
  1249. geom_bar(position="dodge", stat="identity", color="black") +
  1250. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1251. theme_classic() +
  1252. scale_fill_viridis_d(option="inferno", alpha=0.6, begin=0.2, end=0.8)
  1254. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% amyloid_proteins, ] %>%
  1255. ggplot(aes(fill=Group, x=SYMBOL, y=log2FoldChange)) +
  1256. geom_bar(position="dodge", stat="identity", color="black") +
  1257. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1258. theme_classic() +
  1259. scale_fill_viridis_d(option="inferno", alpha=0.6, begin=0.2, end=0.8)
  1261. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% tight_junctions, ] %>%
  1262. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1263. geom_bar(position="dodge", stat="identity", color="black") +
  1264. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1265. theme_classic() +
  1266. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1267. coord_flip()
  1269. #combine these together
  1270. slc_atp_app_tj<-c(SLC_family_genes, ATP_synthase_family, amyloid_proteins, tight_junctions, "Aqp1")
  1271. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% slc_atp_app_tj, ] %>%
  1272. ggplot(aes(fill=Group, x= SYMBOL, y=log2FoldChange)) +
  1273. geom_bar(position="dodge", stat="identity", color="black") +
  1274. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1275. theme_classic() +
  1276. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1277. coord_flip()
  1278. #saved as "barplot_young_old_ad_slc_family_genes.pdf" 5x4
  1280. #rotated
  1281. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% slc_atp_app_tj, ] %>%
  1282. ggplot(aes(fill=Group, x= SYMBOL, y=log2FoldChange)) +
  1283. geom_bar(position="dodge", stat="identity", color="black") +
  1284. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1285. theme_classic() +
  1286. theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  1287. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7))
  1288. #saved as "barplot_young_old_ad_slc_atp_tj_aqp1_app_family_genes_rotate.pdf" 5x4
  1291. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% age_gene_names, ] %>%
  1292. ggplot(aes(fill=Group, x=SYMBOL, y=log2FoldChange)) +
  1293. geom_bar(position="dodge", stat="identity", color="black") +
  1294. geom_text(aes(label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1295. theme_classic() +
  1296. theme(axis.text.x=element_text(size= 2, angle = 90)) +
  1297. scale_fill_viridis_d(option="inferno", alpha=0.6, begin=0.2, end=0.8) +
  1298. coord_flip()
  1299. #this is too many genes. Need to filter by gene and have only those that are significant at least in one group be plotted
  1300. #slc12a2
  1301. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% "Slc12a2", ] %>%
  1302. ggplot(aes(fill=Group, x=SYMBOL, y=log2FoldChange)) +
  1303. geom_bar(position="dodge", stat="identity", color="black") +
  1304. geom_text(aes(label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1305. theme_classic() +
  1306. theme(axis.text.x=element_text(size= 2, angle = 90)) +
  1307. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7))
  1310. young_old_ad_filtered_bysig[young_old_ad_filtered_bysig$SYMBOL %in% "slc12a2", ] %>%
  1311. ggplot(aes(fill=Group, x=reorder(SYMBOL, -log2FoldChange), y=log2FoldChange)) +
  1312. geom_bar(position="dodge", stat="identity", color="black") +
  1313. geom_text(aes(y=half_fc, label= padj_symbol), size=3, vjust=1, position=position_dodge(0.9)) +
  1314. theme_classic() +
  1315. scale_fill_manual(values = alpha(c("darkcyan", 'midnightblue', 'orangered1'), 0.7)) +
  1316. coord_flip()
  1317. ```
  1319. #Venn diagrams to visualize overlaps
  1320. ```{r venn diagrams of gene list overlaps}
  1322. library(VennDiagram)
  1323. sig_res_3moAWvsSL_genelist<-sig_res_3moAWvsSL$SYMBOL
  1324. sig_res_12moWT_df_genelist<-sig_res_12moWT_df$SYMBOL
  1325. sig_res_12moAD_df_genelist<-sig_res_12moAD_df$SYMBOL
  1327. venn.diagram(x=list(sig_res_3moAWvsSL_genelist, sig_res_12moWT_df_genelist, sig_res_12moAD_df_genelist),
  1328. category.names=c("Young", "Aged", "AD"),
  1329. filename = "Circadian Gene Overlaps Young Aged AD.png",
  1330. output=TRUE,
  1331. imagetype="png" ,
  1332. height = 1000 ,
  1333. width = 1000 ,
  1334. resolution = 1000,
  1335. compression = "lzw",
  1336. lwd = 1,
  1337. col=c("darkcyan", 'midnightblue', 'orangered1'),
  1338. fill = c(alpha("darkcyan",0.3), alpha('midnightblue',0.3), alpha('orangered1',0.3)),
  1339. cex = 0.5,
  1340. fontfamily = "sans",
  1341. cat.cex = 0.3,
  1342. cat.default.pos = "outer",
  1343. cat.pos = c(-27, 27, 135),
  1344. cat.dist = c(0.055, 0.055, 0.085),
  1345. cat.fontfamily = "sans",
  1346. cat.col = c("darkcyan", 'midnightblue', 'orangered1'),
  1347. rotation = 1
  1348. )
  1349. #for young, AWvsSD, AWvsSL
  1350. sig_res_3moAWvsSL_genelist<-sig_res_3moAWvsSL$SYMBOL
  1351. sig_res_3moAWvsSD_genelist<-sig_res_3moAWvsSD$SYMBOL
  1353. venn.diagram(
  1354. x=list(sig_res_3moAWvsSL_genelist,sig_res_3moAWvsSD_genelist),
  1355. category.names=c("SL", "SD"),
  1356. filename = "Circadian Gene Overlaps SL vs SD.png",
  1357. output=TRUE,
  1358. imagetype="png" ,
  1359. height = 1000 ,
  1360. width = 1000 ,
  1361. resolution = 1000,
  1362. compression = "lzw",
  1363. lwd = 1,
  1364. col=c("midnightblue", 'orangered1'),
  1365. fill = c(alpha("midnightblue",0.3), alpha('orangered1',0.3)),
  1366. cex = 0.4,
  1367. fontfamily = "sans",
  1368. cat.cex = 0.4,
  1369. cat.default.pos = "outer",
  1370. cat.fontfamily = "sans",
  1371. cat.dist = c(0.055, 0.055),
  1372. cat.pos = c(-50, 30),
  1373. cat.col = c("midnightblue", 'orangered1')
  1374. )
  1375. ```
  1378. ```{r heatmaps}
  1379. #heatmap?
  1380. young_old_ad_mat<-as.matrix(wide_mod[ ,c("log2FoldChange_Young", "log2FoldChange_OldWT", "log2FoldChange_OldAD") ])
  1381. rownames(young_old_ad_mat)<-wide_mod$SYMBOL
  1383. #make matrix of p values by taking the padj symbols
  1384. wide_sig_symbol<-pivot_wider(young_old_ad_symbol[, c(1,4:5)], names_from = Group, values_from = padj_symbol)
  1385. #then have to filter by ones that are significant in at least one group
  1386. wide_filtered_sig<-wide_sig_symbol[wide_sig_symbol$SYMBOL %in% sig_list, ]
  1387. #make into a matrix
  1388. wide_filtered_sig_mat<-as.matrix(wide_filtered_sig[, c(2:4)])
  1389. rownames(wide_filtered_sig_mat)<-wide_filtered_sig$SYMBOL
  1390. #then get rid of the NA values so they don't show up on the heatmap
  1391. wide_filtered_sig_mat[is.na(wide_filtered_sig_mat)]<-""
  1393. pheatmap(young_old_ad_mat[rownames(young_old_ad_mat) %in% age_gene_names, ], cluster_cols = FALSE, display_numbers = wide_filtered_sig_mat[rownames(wide_filtered_sig_mat) %in% age_gene_names, ])
  1396. ```

deseq2 feb 2023.Rmd, under CC-BY-4.0 · at the source

Overview

Authors: Deidre Jansson1,2,3, Ryan O’Boyle2,3, Taylor J Pedersen2,3, Elizabeth Gino2,3, Mathew Sevao2,3, Ron Vered2,3, Katherine L Suchland3, Blake Zhou4, Ryann Fame4, Samantha A Keil2,3,5, Molly Braun2,3,6, Jeffrey Iliff2,3,7
ORCID iDs: Deidre Jansson, Ryan O’Boyle, Ron Vered, Blake Zhou, Ryann Fame
  1. School of Biological Sciences, University of Auckland, Auckland, New Zealand
  2. Department of Psychiatry and Behavioral Sciences, University of Washington School of Medicine, Seattle, WA USA
  3. VA Northwest Mental Illness Research, Education, and Clinical Center (MIRECC), VA Puget Sound Health Care System, Seattle, WA USA
  4. Department of Neurosurgery, Stanford University, Stanford, CA USA
  5. Department of Radiology, Brain Health Imaging Institute, Weill Cornell Medicine, New York, NY USA
  6. Department of Neurosurgery, Medical College of Georgia, Augusta University, Augusta, GA USA
  7. Department of Neurology, University of Washington School of Medicine, Seattle, WA USA
Institutions: University of Auckland (New Zealand); University of Washington (United States); VA Puget Sound Health Care System (United States); Mental Illness Research, Education and Clinical Centers (United States); Stanford Medicine (United States); Stanford University (United States); Weill Cornell Medicine (United States); Augusta University (United States)
Journal: Communications biology, volume 9, issue 1, article 1073
Dates: received 19 November 2025; accepted 13 May 2026; published online 20 May 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s42003-026-10335-4 · PMID 42162350 · PMCID PMC13457615 · OpenAlex W7161808160
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), mouse (organism), Alzheimer's / dementia (population), cellular / molecular (subfield)
Methods: Statistics
Keywords: Circadian rhythms and sleep, Ion channels in the nervous system
MeSH: Aging*, Amyloid beta-Peptides*, Choroid Plexus*, Circadian Rhythm*, Age Factors, Alzheimer Disease, Animals, Female, Male, Mice, Sex Factors, Solute Carrier Family 12, Member 2, Transcriptome (* major topic)
Topic: Cerebrospinal fluid and hydrocephalus (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: U.S. Department of Health &amp; Human Services | NIH | National Institute on Aging (P30 AG066509); NIA NIH HHS (P30 AG066509)
Citations: cited by 1 paper (Europe PMC); 91 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

figshare 32129572

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (5 files), ggplot2 (4 files), clusterProfiler (2 files), ggpubr (2 files), pheatmap (2 files), rstatix (2 files), ComplexHeatmap (1 file), DESeq2 (1 file), patchwork (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
7 files
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s42003-026-10335-4.

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;
  • 6 scripts, each with its path and the digest of its content;
  • 9 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

Datasets cited

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s42003-026-10335-4.

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

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 2 keywords, 13 MeSH terms, 2 funders, 89 references.

Cite

This paper

Jansson, D., O’Boyle, R., Pedersen, T. J., Gino, E., Sevao, M., Vered, R., Suchland, K. L., Zhou, B., Fame, R., Keil, S. A., Braun, M., & Iliff, J. (2026). Diurnal choroid plexus function in mice depends on sex, age, and amyloid-β status. Communications biology, 9(1), 1073. https://doi.org/10.1038/s42003-026-10335-4

BibTeX

@article{jansson2026diurnal,
author = {Jansson, Deidre and O’Boyle, Ryan and Pedersen, Taylor J and Gino, Elizabeth and Sevao, Mathew and Vered, Ron and Suchland, Katherine L and Zhou, Blake and Fame, Ryann and Keil, Samantha A and Braun, Molly and Iliff, Jeffrey},
title = {{Diurnal choroid plexus function in mice depends on sex, age, and amyloid-β status}},
journal = {Communications biology},
year = {2026},
month = may,
volume = {9},
number = {1},
pages = {1073},
publisher = {Nature Publishing Group},
issn = {2399-3642},
doi = {10.1038/s42003-026-10335-4},
url = {https://doi.org/10.1038/s42003-026-10335-4},
pmid = {42162350},
pmcid = {PMC13457615}
}

RIS

TY - JOUR
AU - Jansson, Deidre
AU - O’Boyle, Ryan
AU - Pedersen, Taylor J
AU - Gino, Elizabeth
AU - Sevao, Mathew
AU - Vered, Ron
AU - Suchland, Katherine L
AU - Zhou, Blake
AU - Fame, Ryann
AU - Keil, Samantha A
AU - Braun, Molly
AU - Iliff, Jeffrey
TI - Diurnal choroid plexus function in mice depends on sex, age, and amyloid-β status
T2 - Communications biology
J2 - Commun Biol
PY - 2026
DA - 2026/05/20
VL - 9
IS - 1
SP - 1073
SN - 2399-3642
PB - Nature Publishing Group
DO - 10.1038/s42003-026-10335-4
UR - https://doi.org/10.1038/s42003-026-10335-4
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s42003-026-10335-4",
"type": "article-journal",
"title": "Diurnal choroid plexus function in mice depends on sex, age, and amyloid-β status",
"container-title": "Communications biology",
"author": [
{
"family": "Jansson",
"given": "Deidre"
},
{
"family": "O’Boyle",
"given": "Ryan"
},
{
"family": "Pedersen",
"given": "Taylor J"
},
{
"family": "Gino",
"given": "Elizabeth"
},
{
"family": "Sevao",
"given": "Mathew"
},
{
"family": "Vered",
"given": "Ron"
},
{
"family": "Suchland",
"given": "Katherine L"
},
{
"family": "Zhou",
"given": "Blake"
},
{
"family": "Fame",
"given": "Ryann"
},
{
"family": "Keil",
"given": "Samantha A"
},
{
"family": "Braun",
"given": "Molly"
},
{
"family": "Iliff",
"given": "Jeffrey"
}
],
"container-title-short": "Commun Biol",
"volume": "9",
"issue": "1",
"page": "1073",
"DOI": "10.1038/s42003-026-10335-4",
"PMID": "42162350",
"PMCID": "PMC13457615",
"ISSN": "2399-3642",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s42003-026-10335-4",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
20
]
]
}
}

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, DESeq2, clusterProfiler, 7 other tools, genetics / omics, mouse, cellular / molecular
[2] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: DESeq2, clusterProfiler, ComplexHeatmap, 6 other tools, mouse, cellular / molecular, 1 reference
[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, DESeq2, clusterProfiler, 7 other tools, genetics / omics
[4] doi:10.3390/ijms27093997 [code]
Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC.
Journal: International journal of molecular sciences
In common: rstatix, DESeq2, clusterProfiler, 7 other tools, genetics / omics
[5] doi:10.1126/sciadv.aeb4265 [code]
Single-nucleus profiling reveals a core disease signature and cell type-specific vulnerabilities in early Rett syndrome.
Journal: Science advances
In common: rstatix, DESeq2, clusterProfiler, 6 other tools, genetics / omics, mouse, cellular / molecular
[6] 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, DESeq2, clusterProfiler, 6 other tools, genetics / omics, mouse, cellular / molecular
[7] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: rstatix, DESeq2, clusterProfiler, 6 other tools, genetics / omics, mouse, cellular / molecular
[8] doi:10.1172/jci.insight.207270 [code]
Progressive hypothalamic neuroinflammation in ovariectomized mice parallels aging-related transcriptomic changes in the female human hypothalamus.
Journal: JCI insight
In common: DESeq2, clusterProfiler, ComplexHeatmap, 6 other tools, genetics / omics, mouse, cellular / molecular
[9] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: DESeq2, clusterProfiler, ComplexHeatmap, 6 other tools, genetics / omics, mouse, cellular / molecular
[10] 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, clusterProfiler, ComplexHeatmap, 6 other tools, genetics / omics, cellular / molecular

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.