OSCR

Simultaneous profiling of native-state proteomes and transcriptomes of neural cell types using proximity labeling.

Code ↔ Paper

5 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 5 matches
  1. [1] § Methods › Data analysis and visualization ↔ Astrocyte-RiboTag_mRNA seq analysis.R, lines 1407 to 1451 · score 0.79 · Bader Lab website, GOparallel, DESeq2, GMT, ANOVA, GSEA · Report a mismatch
  2. [2] § Methods › Protein identification and quantification ↔ HEK293 DIA LFQ-MS analysis.R, lines 1 to 41 · score 0.75 · MaxQuant, log2 transformation, tryptic, DDA, UniProt, DIA · Report a mismatch
  3. [3] § Methods › Data analysis and visualization ↔ HEK293 DIA LFQ-MS analysis.R, lines 1542 to 1602 · score 0.74 · Bader Lab website, GOparallel, GMT, GSEA, Fisher, symbol · Report a mismatch
  4. [4] § Methods › RNA-sequencing and processing ↔ Astrocyte-RiboTag_mRNA seq analysis.R, lines 1 to 42 · score 0.70 · smRNA, DESeq2, gene expression, RiboTag, Spliced, library · Report a mismatch
  5. [5] § Methods › RNA-sequencing and processing ↔ Astrocyte and neuron-TurboID mRNA-seq analysis.R, lines 1 to 43 · score 0.63 · smRNA, DESeq2, gene expression, Spliced, library, raw · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 2,134 lines · 96 KB · MIT · 2 matches

  1. ########################################################################################################
  2. # Code for performing Differential Gene Expression Analysis on Exp. 8-1 (Aldh1l1-RiboTag) and Exp. 9-2 (Aldh1l1-TurboID)
  3. # row means >=10 filter before dds and splicing
  4. # Christina Ramelow, MS
  5. # Date: 06/18/2024
  6. ##################################################################################################################
  7. ##################################################################################################################
  8. ## Part 1: Data processing and clean up
  9. # steps:
  10. # 1. Filter out low counts
  11. # a. Venns: keep a gene if row mean count >= 10/sample group
  12. # b. PCAs, correlations, diffex: keep a gene if row mean count >= 10 across all samples
  13. # 2. DESeq2 normalization
  14. # a. median of ratios approach (one file with all samples and 1 file w/o neg.pulldowns)
  15. ##################################################################################################################
  16. #load Libraries
  17. library(DESeq2)
  18. library(tidyverse)
  19. library(biomaRt)
  20. library(pheatmap)
  21. library(ggplot2)
  22. library(ggrepel)
  24. stringsAsFactors = FALSE
  26. rootdir <- setwd("/Users/christina/Desktop/Exp. 8-1 mRNA-seq analysis")
  28. # load Exp. 9-2 feature_counts.txt file
  29. raw_counts_Aldh_Turbo <- read.table("Exp. 1-4, 9-2 and 10-1_feature_counts.txt", sep = "\t", header = TRUE, row.names = 1)
  30. dim(raw_counts_Aldh_Turbo)
  31. # [1] 35976 55
  33. # remove the first 5 columns that are primarily categorical
  34. raw_counts_TurboID_2 <- raw_counts_Aldh_Turbo[,c(6:55)]
  35. dim(raw_counts_TurboID_2)
  36. # [1] 35976 50
  38. # remove the smRNA-seq samples that were accidentally included
  39. raw_counts_TurboID_3 <- raw_counts_TurboID_2[,c(1:20, 27:50)]
  40. dim(raw_counts_TurboID_3)
  41. # [1] 35976 44
  42. # 21047FL-134-01-22 - 21047FL-134-01-27 should be removed. CR: checked
  44. # check if there are NAs in the counts df
  45. sum(is.na(raw_counts_TurboID_3))
  46. #[1] 0
  48. # create a geneID column and move it to the front of the df
  49. raw_counts_TurboID_3$geneID <- rownames(raw_counts_TurboID_3)
  50. raw_counts_TurboID_3 <- raw_counts_TurboID_3[, c("geneID", setdiff(names(raw_counts_TurboID_3), "geneID"))]
  52. # load Exp. 8-1 feature_counts.txt file
  53. raw_counts_Aldh_RiboTag <- read.table("Exp. 8-1 feature_counts.txt", sep = "\t", header = TRUE, row.names = 1)
  54. dim(raw_counts_Aldh_RiboTag)
  55. # [1] 35976 15
  57. # remove the first 5 columns that are primarily categorical
  58. raw_counts_RiboTag_2 <- raw_counts_Aldh_RiboTag[,c(6:15)]
  59. dim(raw_counts_RiboTag_2)
  60. # [1] 35976 10
  62. # check if there are NAs in the counts df
  63. sum(is.na(raw_counts_RiboTag_2))
  64. #[1] 0
  66. # create a geneID column and move it to the front of the df
  67. raw_counts_RiboTag_2$geneID <- rownames(raw_counts_RiboTag_2)
  68. raw_counts_RiboTag_2 <- raw_counts_RiboTag_2[, c("geneID", setdiff(names(raw_counts_RiboTag_2), "geneID"))]
  70. # combine the two featurecounts.txt files by shared geneIDs
  71. TurboID_RiboTag_merged_raw_counts <- merge(raw_counts_TurboID_3, raw_counts_RiboTag_2, by = "geneID")
  72. dim(TurboID_RiboTag_merged_raw_counts)
  73. # [1] 35976 55
  75. # set the geneID column back as row names
  76. rownames(TurboID_RiboTag_merged_raw_counts) <- TurboID_RiboTag_merged_raw_counts$geneID
  78. # remove the geneID column
  79. TurboID_RiboTag_merged_raw_counts <- TurboID_RiboTag_merged_raw_counts[, -which(names(TurboID_RiboTag_merged_raw_counts) == "geneID")]
  81. dim(TurboID_RiboTag_merged_raw_counts)
  82. # [1] 35976 54
  84. # load Exp. 1-4, 9-2, 10-1 and 8-1 RNA traits file
  85. traits <- read.csv("Exp. 1-4, 9-2, 10-1 and 8-1 combined traits_RNA.csv", header = TRUE, row.names = 1)
  86. dim(traits)
  87. # [1] 54 2
  89. # rename the column name of the merged_raw_counts df to match traits df
  90. # Extract row names from traits
  91. new_col_names <- rownames(traits)
  93. # Assign row names of traits as column names
  94. names(TurboID_RiboTag_merged_raw_counts) <- new_col_names
  96. # check that the row names of the traits df and the column names of the counts df are the same
  97. raw_counts_clean <- TurboID_RiboTag_merged_raw_counts[,match(rownames(traits), colnames(TurboID_RiboTag_merged_raw_counts))]
  99. #sanity check 1 colnames in counts file = rownames in traits file
  100. all(colnames(raw_counts_clean) %in% rownames(traits))
  101. #TRUE
  103. #sanity check 2 colnames in counts file = rownames in traits file
  104. all(colnames(raw_counts_clean) == rownames(traits))
  105. #TRUE
  107. # Exp. 9-2 and 8-1 splced
  108. raw_counts_clean_9_2_8_1 <- raw_counts_clean[,c(21:41, 45:54)]
  109. dim(raw_counts_clean_9_2_8_1)
  110. # [1] 35976 31
  112. colnames(raw_counts_clean_9_2_8_1)
  113. # should be "21047FL-134-01-28" - "21047FL-134-01-48" and "23238R-04-01" - "23238R-04-10" CR: checked
  115. # write .csv file for raw_counts_clean for Exp. 9-2 and 8-1
  116. write.csv(raw_counts_clean_9_2_8_1, file=paste0("0.Exp. 9-2 and 8-1_RNA-seq_raw_counts_clean_unfiltered-",dim(raw_counts_clean_9_2_8_1)[1],"x",dim(raw_counts_clean_9_2_8_1)[2],"_CR_06182024.csv"))
  118. ################################################################################
  119. ################################################################################
  120. ## Step 1a: Filter out low counts using a row mean >/= 10 by group for Venn diagrams
  121. # Exp. 9-2 and 8-1
  122. # 1. subset the traits and counts dfs by sample group to create list of genes for Venn diagrams
  123. # 2. calculate the row mean
  124. # 3. filter for genes with a row mean >= 10
  125. # 4. write a .csv file for each sample group
  126. ################################################################################
  127. # subset the traits file to only include Exp. 9-2 and 8-1
  128. traits_9_2_8_1 <- traits[c(21:41, 45:54),]
  129. dim(traits_9_2_8_1)
  130. # [1] 31 2
  131. write.csv(traits_9_2_8_1, file = "Exp. 9-2 and 8-1 RNA traits.csv")
  132. ################################################################################
  133. # row mean >= 10 Aldh1l1.bulk.cortex
  134. counts_Aldh1l1.bulk.cortex <- raw_counts_clean_9_2_8_1[,c(4:6)]
  135. dim(counts_Aldh1l1.bulk.cortex)
  136. # [1] 35976 3
  138. traits_Aldh1l1.bulk.cortex <- traits_9_2_8_1[c(4:6),]
  139. traits_Aldh1l1.bulk.cortex
  140. # Sample.ID.. Group
  141. # 21047FL-134-01-31 4-RNA 30L Aldh1l1.bulk.cortex
  142. # 21047FL-134-01-32 5-RNA 31B Aldh1l1.bulk.cortex
  143. # 21047FL-134-01-33 6-RNA 24L Aldh1l1.bulk.cortex
  145. rownames(traits_Aldh1l1.bulk.cortex)
  146. # [1] "21047FL-134-01-31" "21047FL-134-01-32" "21047FL-134-01-33"
  148. colnames(counts_Aldh1l1.bulk.cortex)
  149. # [1] "21047FL-134-01-31" "21047FL-134-01-32" "21047FL-134-01-33"
  151. # calculate row mean
  152. Aldh1l1.bulk.cortex_row_mean <- rowMeans(counts_Aldh1l1.bulk.cortex)
  154. # subset rows where row means >= 10
  155. Aldh1l1.bulk.cortex_filtered_rows <- counts_Aldh1l1.bulk.cortex[Aldh1l1.bulk.cortex_row_mean >= 10, ]
  157. dim(Aldh1l1.bulk.cortex_filtered_rows)
  158. # [1] 16737 3
  160. # write .csv file
  161. write.csv(Aldh1l1.bulk.cortex_filtered_rows, file=paste0("1a.Exp. 9-2 and 8-1_RNA-seq_Aldh1l1.bulk.cortex_Venn filter-",dim(Aldh1l1.bulk.cortex_filtered_rows)[1],"x",dim(Aldh1l1.bulk.cortex_filtered_rows)[2],"_CR_06182024.csv"))
  163. ################################################################################
  164. # row mean >= 10 RiboTag bulk.cortex
  165. counts_RiboTag.bulk.cortex <- raw_counts_clean_9_2_8_1[,c(22:23)]
  166. dim(counts_RiboTag.bulk.cortex)
  167. # [1] 35976 2
  169. traits_RiboTag.bulk.cortex <- traits_9_2_8_1[c(22:23),]
  170. traits_RiboTag.bulk.cortex
  171. # Sample.ID.. Group
  172. # 23238R-04-01 CR1 RiboTag.bulk.cortex
  173. # 23238R-04-02 CR2 RiboTag.bulk.cortex
  175. rownames(traits_RiboTag.bulk.cortex)
  176. # [1] "23238R-04-01" "23238R-04-02"
  178. colnames(counts_RiboTag.bulk.cortex)
  179. # [1] "23238R-04-01" "23238R-04-02"
  181. # calculate row mean
  182. RiboTag.bulk.cortex_row_mean <- rowMeans(counts_RiboTag.bulk.cortex)
  184. # subset rows where row means >= 10
  185. RiboTag.bulk.cortex_filtered_rows <- counts_RiboTag.bulk.cortex[RiboTag.bulk.cortex_row_mean >= 10, ]
  187. dim(RiboTag.bulk.cortex_filtered_rows)
  188. # [1] 17584 2
  190. # write .csv file
  191. write.csv(RiboTag.bulk.cortex_filtered_rows, file=paste0("1a.Exp. 9-2 and 8-1_RNA-seq_RiboTag.bulk.cortex_Venn filter-",dim(RiboTag.bulk.cortex_filtered_rows)[1],"x",dim(RiboTag.bulk.cortex_filtered_rows)[2],"_CR_06182024.csv"))
  193. ###############################################################################
  194. # row mean >= 10 Aldh1l1.RiboTag.bulk.cortex
  195. counts_Aldh1l1.RiboTag.bulk.cortex <- raw_counts_clean_9_2_8_1[,c(24:26)]
  196. dim(counts_Aldh1l1.RiboTag.bulk.cortex)
  197. # [1] 35976 3
  198. traits_Aldh1l1.RiboTag.bulk.cortex <- traits_9_2_8_1[c(24:26),]
  199. traits_Aldh1l1.RiboTag.bulk.cortex
  200. # Sample.ID.. Group
  201. # 23238R-04-03 CR3 Aldh1l1.RiboTag.bulk.cortex
  202. # 23238R-04-04 CR4 Aldh1l1.RiboTag.bulk.cortex
  203. # 23238R-04-05 CR5 Aldh1l1.RiboTag.bulk.cortex
  205. rownames(traits_Aldh1l1.RiboTag.bulk.cortex)
  206. # [1] "23238R-04-03" "23238R-04-04" "23238R-04-05"
  208. colnames(counts_Aldh1l1.RiboTag.bulk.cortex)
  209. # [1] "23238R-04-03" "23238R-04-04" "23238R-04-05"
  211. # calculate row mean
  212. Aldh1l1.RiboTag.bulk.cortex_row_mean <- rowMeans(counts_Aldh1l1.RiboTag.bulk.cortex)
  214. # subset rows where row means >= 10
  215. Aldh1l1.RiboTag.bulk.cortex_filtered_rows <- counts_Aldh1l1.RiboTag.bulk.cortex[Aldh1l1.RiboTag.bulk.cortex_row_mean >= 10, ]
  217. dim(Aldh1l1.RiboTag.bulk.cortex_filtered_rows)
  218. # [1] 17583 3
  220. # write .csv file
  221. write.csv(Aldh1l1.RiboTag.bulk.cortex_filtered_rows, file=paste0("1a.Exp. 9-2 and 8-1_RNA-seq_Aldh1l1.RiboTag.bulk.cortex_Venn filter-",dim(Aldh1l1.RiboTag.bulk.cortex_filtered_rows)[1],"x",dim(Aldh1l1.RiboTag.bulk.cortex_filtered_rows)[2],"_CR_06182024.csv"))
  223. ################################################################################
  224. # row mean >= 10 Aldh1l1.pos.pulldown
  225. counts_Aldh1l1.pos.pulldown <- raw_counts_clean_9_2_8_1[,c(14:16)]
  226. dim(counts_Aldh1l1.pos.pulldown)
  227. # [1] 35976 3
  228. traits_Aldh1l1.pos.pulldown <- traits_9_2_8_1[c(14:16),]
  229. traits_Aldh1l1.pos.pulldown
  230. # Sample.ID.. Group
  231. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  232. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  233. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  235. rownames(traits_Aldh1l1.pos.pulldown)
  236. # [1] "21047FL-134-01-41" "21047FL-134-01-42" "21047FL-134-01-43"
  238. colnames(counts_Aldh1l1.pos.pulldown)
  239. # [1] "21047FL-134-01-41" "21047FL-134-01-42" "21047FL-134-01-43"
  241. # calculate row mean
  242. Aldh1l1.pos.pulldown_row_mean <- rowMeans(counts_Aldh1l1.pos.pulldown)
  244. # subset rows where row means >= 10
  245. Aldh1l1.pos.pulldown_filtered_rows <- counts_Aldh1l1.pos.pulldown[Aldh1l1.pos.pulldown_row_mean >= 10, ]
  247. dim(Aldh1l1.pos.pulldown_filtered_rows)
  248. # [1] 16632 3
  250. # write .csv file
  251. write.csv(Aldh1l1.pos.pulldown_filtered_rows, file=paste0("1a.Exp. 9-2 and 8-1_RNA-seq_Aldh1l1.pos.pulldown_Venn filter-",dim(Aldh1l1.pos.pulldown_filtered_rows)[1],"x",dim(Aldh1l1.pos.pulldown_filtered_rows)[2],"_CR_06182024.csv"))
  253. ################################################################################
  254. # row mean >= 10 Aldh1l1.RiboTag.pos.pulldown
  255. counts_Aldh1l1.RiboTag.pos.pulldown<- raw_counts_clean_9_2_8_1[,c(29:31)]
  256. dim(counts_Aldh1l1.RiboTag.pos.pulldown)
  257. # [1] 35976 3
  259. traits_Aldh1l1.RiboTag.pos.pulldown<- traits_9_2_8_1[c(29:31),]
  260. traits_Aldh1l1.RiboTag.pos.pulldown
  261. # Sample.ID.. Group
  262. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  263. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  264. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  266. rownames(traits_Aldh1l1.RiboTag.pos.pulldown)
  267. # [1] "23238R-04-08" "23238R-04-09" "23238R-04-10"
  269. colnames(counts_Aldh1l1.RiboTag.pos.pulldown)
  270. # [1] "23238R-04-08" "23238R-04-09" "23238R-04-10"
  272. # calculate row mean
  273. Aldh1l1.RiboTag.pos.pulldown_row_mean <- rowMeans(counts_Aldh1l1.RiboTag.pos.pulldown)
  275. # subset rows where row means >= 10
  276. Aldh1l1.RiboTag.pos.pulldown_filtered_rows <- counts_Aldh1l1.RiboTag.pos.pulldown[Aldh1l1.RiboTag.pos.pulldown_row_mean >= 10, ]
  278. dim(Aldh1l1.RiboTag.pos.pulldown_filtered_rows)
  279. # [1] 16737 3
  281. # write .csv file
  282. write.csv(Aldh1l1.RiboTag.pos.pulldown_filtered_rows, file=paste0("1a.Exp. 9-2 and 8-1_RNA-seq_Aldh1l1.RiboTag.pos.pulldown_Venn filter-",dim(Aldh1l1.RiboTag.pos.pulldown_filtered_rows)[1],"x",dim(Aldh1l1.RiboTag.pos.pulldown_filtered_rows)[2],"_CR_06182024.csv"))
  285. ################################################################################
  286. ################################################################################
  287. # clean up filtered row dfs to only include gene lists
  289. # Aldh1l1.bulk.cortex_filtered_rows
  290. Aldh1l1.bulk.cortex_genes <- as.matrix(rownames(Aldh1l1.bulk.cortex_filtered_rows))
  291. dim(Aldh1l1.bulk.cortex_genes)
  292. # [1] 16737 1
  294. # Aldh1l1.RiboTag.bulk.cortex_filtered_rows
  295. Aldh1l1.RiboTag.bulk.cortex_genes <- as.matrix(rownames(Aldh1l1.RiboTag.bulk.cortex_filtered_rows))
  296. dim(Aldh1l1.RiboTag.bulk.cortex_genes)
  297. # [1] 17583 1
  299. # RiboTag.bulk.cortex_filtered_rows
  300. RiboTag.bulk.cortex_genes <- as.matrix(rownames(RiboTag.bulk.cortex_filtered_rows))
  301. dim(RiboTag.bulk.cortex_genes)
  302. # [1] 17584 1
  304. # Aldh1l1.pos.pulldown_filtered_rows
  305. Aldh1l1.pos.pulldown_genes <- as.matrix(rownames(Aldh1l1.pos.pulldown_filtered_rows))
  306. dim(Aldh1l1.pos.pulldown_genes)
  307. # [1] 16632 1
  309. # Aldh1l1.RiboTag.pos.pulldown_filtered_rows
  310. Aldh1l1.RiboTag.pos.pulldown_genes <- as.matrix(rownames(Aldh1l1.RiboTag.pos.pulldown_filtered_rows))
  311. dim(Aldh1l1.RiboTag.pos.pulldown_genes)
  312. # [1] 16737 1
  314. ################################################################################
  315. # Create Venn diagrams for each comparison
  316. list_z <- NULL
  318. # #import fonts to use Arial in output pdf
  319. # https://r-coder.com/custom-fonts-r/
  320. library(extrafont)
  321. font_import()
  322. fonts()
  323. loadfonts(device = "pdf", quiet = TRUE)
  325. # load library to create Venn's
  326. # GUI: https://www.biovenn.nl/
  327. # R package: https://cran.r-project.org/web/packages/BioVenn/index.html
  328. library(BioVenn)
  330. ## Create Venn diagrams and output as a pdf ##
  331. pdf(file = "1a. Exp. 9-2 and 8-1_RNA-seq_Venns_CR_06182024.pdf", height = 11, width = 8.5, family = "Arial")
  332. par(mfrow= c(2,1)) # centers plot on single page (I think)
  334. # Aldh1l1.bulk.cortex vs Aldh1l1.RiboTag.bulk.cortex
  335. biovenn1 <- draw.venn(Aldh1l1.bulk.cortex_genes, Aldh1l1.RiboTag.bulk.cortex_genes, list_z,
  336. t_fb = 2,
  337. t_s = 1.5,
  338. t_c = "black",
  339. t_f = "Arial",
  340. title="Aldh1l1.TurboID vs Aldh1l1.RiboTag bulk.cortex RNA",
  341. xtitle = "Aldh1l1.bulk.cortex",
  342. xt_f = "Arial",
  343. xt_fb = 2,
  344. xt_s = 1,
  345. xt_c = "black",
  346. ytitle = "Aldh1l1.RiboTag.bulk.cortex",
  347. yt_f = "Arial",
  348. yt_fb = 2,
  349. yt_s = 1,
  350. yt_c = "black",
  351. nrtype = "abs",
  352. nr_f = "Arial",
  353. nr_fb = 2,
  354. nr_s = 1,
  355. nr_c = "black",
  356. x_c = "purple",
  357. y_c = "red",
  358. # z_c = "blue",
  359. bg_c = "white",
  360. width = 1000,
  361. height = 1000,
  362. #output = "screen",
  363. filename = NULL,
  364. map2ens = FALSE)
  365. # [1] "x bulk.cortex: 16737"
  366. # [1] "y bulk.cortex: 17583"
  367. # [1] "z bulk.cortex: 0"
  368. # [1] "x only: 245"
  369. # [1] "y only: 1091"
  370. # [1] "z only: 0"
  371. # [1] "x-y bulk.cortex overlap: 16492"
  373. # Aldh1l1.RiboTag.bulk.cortex vs Aldh1l1.RiboTag.pos.pulldown
  374. biovenn2 <- draw.venn(Aldh1l1.RiboTag.bulk.cortex_genes, Aldh1l1.RiboTag.pos.pulldown_genes, list_z,
  375. t_fb = 2,
  376. t_s = 1.5,
  377. t_c = "black",
  378. t_f = "Arial",
  379. title="Aldh1l1.RiboTag.bulk.cortex vs Aldh1l1.RiboTag.pos.pulldown RNA",
  380. xtitle = "Aldh1l1.RiboTag.bulk.cortex",
  381. xt_f = "Arial",
  382. xt_fb = 2,
  383. xt_s = 1,
  384. xt_c = "black",
  385. ytitle = "Aldh1l1.RiboTag.pos.pulldown",
  386. yt_f = "Arial",
  387. yt_fb = 2,
  388. yt_s = 1,
  389. yt_c = "black",
  390. nrtype = "abs",
  391. nr_f = "Arial",
  392. nr_fb = 2,
  393. nr_s = 1,
  394. nr_c = "black",
  395. x_c = "grey",
  396. y_c = "red",
  397. # z_c = "blue",
  398. bg_c = "white",
  399. width = 1000,
  400. height = 1000,
  401. #output = "screen",
  402. filename = NULL,
  403. map2ens = FALSE)
  404. # [1] "x bulk.cortex: 17583"
  405. # [1] "y bulk.cortex: 16737"
  406. # [1] "z bulk.cortex: 0"
  407. # [1] "x only: 1319"
  408. # [1] "y only: 473"
  409. # [1] "z only: 0"
  410. # [1] "x-y bulk.cortex overlap: 16264"
  413. # Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown
  414. biovenn3 <- draw.venn(Aldh1l1.pos.pulldown_genes, Aldh1l1.RiboTag.pos.pulldown_genes, list_z,
  415. t_fb = 2,
  416. t_s = 1.5,
  417. t_c = "black",
  418. t_f = "Arial",
  419. title="Aldh1l1.TurboID vs Aldh1l1.RiboTag pulldown RNA",
  420. xtitle = "Aldh1l1.TurboID",
  421. xt_f = "Arial",
  422. xt_fb = 2,
  423. xt_s = 1,
  424. xt_c = "black",
  425. ytitle = "Aldh1l1.RiboTag",
  426. yt_f = "Arial",
  427. yt_fb = 2,
  428. yt_s = 1,
  429. yt_c = "black",
  430. nrtype = "abs",
  431. nr_f = "Arial",
  432. nr_fb = 2,
  433. nr_s = 1,
  434. nr_c = "black",
  435. x_c = "purple4",
  436. y_c = "red4",
  437. # z_c = "blue",
  438. bg_c = "white",
  439. width = 1000,
  440. height = 1000,
  441. #output = "screen",
  442. filename = NULL,
  443. map2ens = FALSE)
  445. # [1] "x bulk.cortex: 16632"
  446. # [1] "y bulk.cortex: 16737"
  447. # [1] "z bulk.cortex: 0"
  448. # [1] "x only: 666"
  449. # [1] "y only: 771"
  450. # [1] "z only: 0"
  451. # [1] "x-y bulk.cortex overlap: 15966"
  452. # [1] "x-z bulk.cortex overlap: 0"
  453. # [1] "y-z bulk.cortex overlap: 0"
  454. # [1] "x-y only overlap: 15966"
  457. dev.off()
  459. ################################################################################
  460. ################################################################################
  461. ## Step 1b: Create row mean filter >=10 across all samples for downstream diffex analysis
  463. # calculate row mean
  464. counts_clean_row_mean <- rowMeans(raw_counts_clean_9_2_8_1)
  466. # subset rows where row means >= 10
  467. counts_clean_filtered_rows_9_2_8_1 <- raw_counts_clean_9_2_8_1[counts_clean_row_mean >= 10, ]
  469. dim(counts_clean_filtered_rows_9_2_8_1)
  470. # [1] 17054 31
  472. # write .csv file
  473. write.csv(counts_clean_filtered_rows_9_2_8_1, file=paste0("1b.Exp. 9-2 and 8-1_RNA-seq_clean_counts_diffex filter-",dim(counts_clean_filtered_rows_9_2_8_1)[1],"x",dim(counts_clean_filtered_rows_9_2_8_1)[2],"_CR_06182024.csv"))
  475. ################################################################################
  476. ################################################################################
  477. ## Step 2: DESeq2 normalization of counts based on median of ratios
  479. # load library
  480. library(DESeq2)
  482. # create DESeq data set matrix (dds) for all 31 samples
  483. dds_all9_2_8_1 <- DESeqDataSetFromMatrix(countData = counts_clean_filtered_rows_9_2_8_1,
  484. colData = traits_9_2_8_1,
  485. design = ~Group)
  486. dim(dds_all9_2_8_1)
  487. # [1] 17054 31
  489. # set the factor level
  490. dds_all9_2_8_1$Group <- relevel(dds_all9_2_8_1$Group, ref = "Aldh1l1.bulk.cortex")
  492. # run DESeq
  493. dds_all9_2_8_1<- DESeq(dds_all9_2_8_1)
  495. View(counts(dds_all9_2_8_1))
  497. dds <- dds_all9_2_8_1
  499. dds <- estimateSizeFactors(dds)
  501. sizeFactors(dds)
  502. # 21047FL-134-01-28 21047FL-134-01-29 21047FL-134-01-30 21047FL-134-01-31 21047FL-134-01-32 21047FL-134-01-33
  503. # 0.8893720 0.8049612 0.7749720 0.8572228 1.0256973 1.1482553
  504. # 21047FL-134-01-34 21047FL-134-01-35 21047FL-134-01-36 21047FL-134-01-37 21047FL-134-01-38 21047FL-134-01-39
  505. # 0.9989146 0.6364181 0.7855086 0.9811747 0.9203750 1.2266314
  506. # 21047FL-134-01-40 21047FL-134-01-41 21047FL-134-01-42 21047FL-134-01-43 21047FL-134-01-44 21047FL-134-01-45
  507. # 1.5524386 1.0284815 1.0962545 0.8934085 1.0540818 0.7051570
  508. # 21047FL-134-01-46 21047FL-134-01-47 21047FL-134-01-48 23238R-04-01 23238R-04-02 23238R-04-03
  509. # 0.8453669 0.8415983 0.8265650 1.1347803 1.2736475 1.3087446
  510. # 23238R-04-04 23238R-04-05 23238R-04-06 23238R-04-07 23238R-04-08 23238R-04-09
  511. # 1.2948424 1.0700477 1.2484744 1.2103432 0.8805125 1.1037337
  512. # 23238R-04-10
  513. # 1.2092379
  515. normalized_counts_1 <- counts(dds, normalized=TRUE)
  517. write.csv(normalized_counts_1, file=paste0("2a.Exp. 9-2 and 8-1_RNA-seq_clean_filtered_DESeq2 norm counts_diffex filter-",dim(normalized_counts_1)[1],"x",dim(normalized_counts_1)[2],"_CR_06182024.csv"))
  519. ################################################################################
  520. # updated clean_counts with neg.pulldowns removed
  521. all9_2_8_1_no_negs_clean_counts <- counts_clean_filtered_rows_9_2_8_1[,c(1:11, 14:26, 29:31)]
  523. # udated traits with neg.pulldowns removed
  524. all9_2_8_1_no_negs_traits <- traits_9_2_8_1[c(1:11, 14:26, 29:31),]
  525. all9_2_8_1_no_negs_traits
  527. # create DESeq data set matrix dds for all 22 samples
  528. dds_all9_2_8_1_no_negs <- DESeqDataSetFromMatrix(countData = all9_2_8_1_no_negs_clean_counts,
  529. colData = all9_2_8_1_no_negs_traits,
  530. design = ~Group)
  531. dim(dds_all9_2_8_1_no_negs)
  532. # [1] 17054 27
  534. # set the factor level
  535. dds_all9_2_8_1_no_negs$Group <- relevel(dds_all9_2_8_1_no_negs$Group, ref = "Aldh1l1.bulk.cortex")
  537. # run DESeq
  538. dds_all9_2_8_1_no_negs<- DESeq(dds_all9_2_8_1_no_negs)
  540. View(counts(dds_all9_2_8_1_no_negs))
  542. dds <- dds_all9_2_8_1_no_negs
  544. dds <- estimateSizeFactors(dds)
  546. sizeFactors(dds)
  547. # 21047FL-134-01-28 21047FL-134-01-29 21047FL-134-01-30 21047FL-134-01-31 21047FL-134-01-32 21047FL-134-01-33
  548. # 0.9213783 0.8349613 0.8023031 0.8899217 1.0636470 1.1919804
  549. # 21047FL-134-01-34 21047FL-134-01-35 21047FL-134-01-36 21047FL-134-01-37 21047FL-134-01-38 21047FL-134-01-41
  550. # 1.0358285 0.6600755 0.8153691 1.0210895 0.9576344 1.0654161
  551. # 21047FL-134-01-42 21047FL-134-01-43 21047FL-134-01-44 21047FL-134-01-45 21047FL-134-01-46 21047FL-134-01-47
  552. # 1.1313599 0.9255927 1.0888734 0.7342807 0.8767168 0.8755041
  553. # 21047FL-134-01-48 23238R-04-01 23238R-04-02 23238R-04-03 23238R-04-04 23238R-04-05
  554. # 0.8588437 1.1870658 1.3348204 1.3673661 1.3577919 1.1199568
  555. # 23238R-04-08 23238R-04-09 23238R-04-10
  556. # 0.9121930 1.1422150 1.2556182
  558. normalized_counts_2 <- counts(dds, normalized=TRUE)
  560. write.csv(normalized_counts_2, file=paste0("2a.Exp. 9-2 and 8-1_RNA-seq_clean_filtered_DESeq2 norm counts_no negs_diffex filter-",dim(normalized_counts_2)[1],"x",dim(normalized_counts_2)[2],"_CR_06182024.csv"))
  562. ## end of data processing and clean up ##
  564. ################################################################################
  565. ################################################################################
  566. ## Part II: Data analysis and visualization
  567. # Exp. 9-2 and 8-1 mRNA-seq analysis
  569. # steps:
  570. # 1. PCAs
  571. # 1b. all groups w/o neg.pulldowns
  572. # 1c. Aldh1l1.pos.pulldown (n=3/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  573. # 2. Correlations
  574. # a. Aldh1l1.bulk.cortex vs Aldh1l1.RiboTag.bulk.cortex (n=2/group)
  575. # b. Aldh1l1.RiboTag.bulk.cortex (n=2/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  576. # c. Aldh1l1.pos.pulldown (n=3/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  577. # 3. Diffex volcanoes & GSEA
  578. # a. Aldh1l1.bulk.cortex vs Aldh1l1.RiboTag.bulk.cortex (n=2/group)
  579. # b. Aldh1l1.RiboTag.bulk.cortex (n=2/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  580. # c. Aldh1l1.pos.pulldown (n=3/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  583. ################################################################################
  584. ## Step 1: Generate PCA plots
  585. library(DESeq2)
  586. library(ggplot2)
  588. # 1a: all groups
  589. # create DESeq data set matrix dds for all 20 samples
  590. dds_all9_2_8_1 <- DESeqDataSetFromMatrix(countData = counts_clean_filtered_rows_9_2_8_1,
  591. colData = traits_9_2_8_1,
  592. design = ~Group)
  593. dim(dds_all9_2_8_1)
  594. # [1] 17054 31
  596. # set the factor level
  597. dds_all9_2_8_1$Group <- relevel(dds_all9_2_8_1$Group, ref = "Aldh1l1.bulk.cortex")
  599. # run DESeq
  600. dds_all9_2_8_1<- DESeq(dds_all9_2_8_1)
  602. # create pdf output file
  603. pdf(file = "1a-1. Exp. 9-2 and 8-1_RNA-seq_all groups_DESeq2_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  604. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  605. heights = c(1.3,0.7,1.3), # Heights of the rows
  606. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  608. # perform DESeq2-basedPCA on rlog transformed data of all 20 samples
  609. rldr <- rlog(dds_all9_2_8_1, blind = TRUE)
  610. plotPCA(rldr, intgroup = "Group") # plots based on normalized counts
  612. pcaData <- plotPCA(rldr, intgroup = "Group", returnData=TRUE)
  613. percentVar <- round(100 * attr(pcaData, "percentVar"))
  614. pca1 <- ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  615. geom_point(size=3) +
  616. xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  617. ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  618. ggtitle("DESeq2 PCA plot for log2(mean mRNA counts)") +
  619. theme(panel.background = element_rect(fill = "transparent", color = "black"),
  620. legend.background = element_rect(fill = "transparent"),
  621. legend.key = element_rect(fill = "transparent"),
  622. plot.title = element_text(hjust = 0.5)) + # Center the title
  623. labs(color = "Transcriptome groups") +
  624. coord_fixed()
  626. print(pca1)
  628. dev.off()
  630. ################################################################################
  631. # load library
  632. library(limma)
  634. # 1b: all groups w/o neg.pulldowns
  635. pdf(file = "1b-1. Exp. 9-2 and 8-1_RNA-seq_all groups_nonegs_DESeq2_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  636. # perform DESeq2-basedPCA on rlog transformed data of 18 samples w/o neg pulldowns
  638. rldr <- rlog(dds_all9_2_8_1_no_negs, blind = TRUE)
  639. plotPCA(rldr, intgroup = "Group") # plots based on normalized counts
  641. pcaData <- plotPCA(rldr, intgroup = "Group", returnData=TRUE)
  642. percentVar <- round(100 * attr(pcaData, "percentVar"))
  643. pca2 <- ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  644. geom_point(size=3) +
  645. xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  646. ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  647. ggtitle("DESeq2 PCA plot for log2(mean mRNA counts)") +
  648. theme(panel.background = element_rect(fill = "transparent", color = "black"),
  649. legend.background = element_rect(fill = "transparent"),
  650. legend.key = element_rect(fill = "transparent"),
  651. plot.title = element_text(hjust = 0.5)) + # Center the title
  652. labs(color = "Transcriptome groups") +
  653. coord_fixed()
  655. print(pca2)
  657. dev.off()
  659. ###############################################################################
  660. # CR didn't edit 06/20/2024
  661. # perform limma:plotMDS-based PCA w/o neg pulldowns
  662. library(limma)
  664. pdf(file = "1b-2. Exp. 9-2 and 8-1_RNA-seq_all groups_nonegs_Limma_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  665. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  666. heights = c(1.3,0.7,1.3), # Heights of the rows
  667. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  669. Grouping_vector <- traits_9_2_8_1$Group[c(1:11, 14:24)]
  671. pch_values <- ifelse(Grouping_vector == "control.bulk.cortex", 16,
  672. ifelse(Grouping_vector == "control.LPS.bulk.cortex", 16,
  673. ifelse(Grouping_vector == "Aldh1l1.bulk.cortex", 16,
  674. ifelse(Grouping_vector == "Aldh1l1.RiboTag.bulk.cortex", 16,
  675. ifelse(Grouping_vector == "Aldh1l1.pos.pulldown", 16,
  676. ifelse(Grouping_vector == "Aldh1l1.LPS.pos.pulldown", 16,
  677. ifelse(Grouping_vector == "RiboTag.bulk.cortex", 16, NA)))))))
  680. pt_colors <- ifelse(Grouping_vector == "control.bulk.cortex", "magenta",
  681. ifelse(Grouping_vector == "control.LPS.bulk.cortex", "pink4",
  682. ifelse(Grouping_vector == "Aldh1l1.bulk.cortex", "coral",
  683. ifelse(Grouping_vector == "Aldh1l1.RiboTag.bulk.cortex", "green3",
  684. ifelse(Grouping_vector == "Aldh1l1.pos.pulldown", "purple2",
  685. ifelse(Grouping_vector == "Aldh1l1.LPS.pos.pulldown", "darkturquoise",
  686. ifelse(Grouping_vector == "RiboTag.bulk.cortex", "dodgerblue", NA)))))))
  688. legend_groups <- c("control.bulk.cortex", "control.LPS.bulk.cortex", "Aldh1l1.bulk.cortex", "Aldh1l1.RiboTag.bulk.cortex", "Aldh1l1.pos.pulldown","Aldh1l1.LPS.pos.pulldown", "RiboTag.bulk.cortex")
  690. plotMDS_1_4_allgroups_nonegs <- plotMDS(log2(normalized_counts_2), #top = 500,
  691. labels = NULL, pch = pch_values, col = pt_colors,
  692. cex = 3, dim.plot = c(1,2), gene.selection = "pairwise",
  693. xlab = NULL, ylab = NULL, plot = TRUE, var.explained = TRUE)
  695. mtext(side=3, text="MDS Plot for log2(mean mRNA counts)", cex=1, line=2)
  697. plot.new()
  698. legend("topright", legend = legend_groups, col = c("magenta", "pink4", "coral", "green3", "purple", "darkturquoise", "dodgerblue"),
  699. pch = 16, title = "Transcriptome groups",cex=1.4)
  701. dev.off()
  703. ###############################################################################
  704. # DESeq2-based PCA for Aldh1l1.bulk.cortex, Aldh1l1.RiboTag.bulk.cortex, Aldh1l1.pos.pulldown and Aldh1l1.RiboTag.pos.pulldown
  706. dim(counts_clean_filtered_rows_9_2_8_1)
  707. # [1] 17054 31
  709. dim(traits_9_2_8_1)
  710. # [1] 31 2
  712. # updated clean_counts
  713. TurboID_Vs_RiboTag.bulk.cortexs_pulls_clean_counts <- counts_clean_filtered_rows_9_2_8_1[,c(4:6, 14:16, 24:26, 29:31)]
  715. # updated traits
  716. TurboID_Vs_RiboTag.bulk.cortexs_pulls_traits <- traits_9_2_8_1[c(4:6, 14:16, 24:26, 29:31),]
  718. # Sample.ID.. Group
  719. # 21047FL-134-01-31 4-RNA 30L Aldh1l1.bulk.cortex
  720. # 21047FL-134-01-32 5-RNA 31B Aldh1l1.bulk.cortex
  721. # 21047FL-134-01-33 6-RNA 24L Aldh1l1.bulk.cortex
  722. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  723. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  724. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  725. # 23238R-04-03 CR3 Aldh1l1.RiboTag.bulk.cortex
  726. # 23238R-04-04 CR4 Aldh1l1.RiboTag.bulk.cortex
  727. # 23238R-04-05 CR5 Aldh1l1.RiboTag.bulk.cortex
  728. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  729. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  730. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  732. # create DESeq data set matrix dds for all 20 samples
  733. dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls <- DESeqDataSetFromMatrix(countData = TurboID_Vs_RiboTag.bulk.cortexs_pulls_clean_counts,
  734. colData = TurboID_Vs_RiboTag.bulk.cortexs_pulls_traits,
  735. design = ~Group)
  736. dim(dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls)
  737. # [1] 17054 12
  739. # set the factor level
  740. dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls$Group <- relevel(dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls$Group, ref = "Aldh1l1.bulk.cortex")
  742. # run DESeq
  743. dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls <- DESeq(dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls)
  745. # create pdf output file
  746. pdf(file = "1c-1. Exp. 9-2 and 8-1_RNA-seq_TurboID_Vs_RiboTag.bulk.cortexs_pulls_DESeq2_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  747. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  748. heights = c(1.3,0.7,1.3), # Heights of the rows
  749. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  751. # perform DESeq2-basedPCA on rlog transformed data of all 20 samples
  752. rldr <- rlog(dds_TurboID_Vs_RiboTag.bulk.cortexs_pulls, blind = TRUE)
  753. plotPCA(rldr, intgroup = "Group") # plots based on normalized counts
  755. pcaData <- plotPCA(rldr, intgroup = "Group", returnData=TRUE)
  756. percentVar <- round(100 * attr(pcaData, "percentVar"))
  757. pca3 <- ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  758. geom_point(size=3) +
  759. xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  760. ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  761. ggtitle("DESeq2 PCA plot for log2(mean mRNA counts)") +
  762. theme(panel.background = element_rect(fill = "transparent", color = "black"),
  763. legend.background = element_rect(fill = "transparent"),
  764. legend.key = element_rect(fill = "transparent"),
  765. plot.title = element_text(hjust = 0.5)) + # Center the title
  766. labs(color = "Transcriptome groups") +
  767. coord_fixed()
  769. print(pca3)
  771. dev.off()
  773. ################################################################################
  774. # 1d Aldh1l1.pos.pulldown vs Aldh1l1.Ribotag.pos.pulldown
  775. # subset DESeq2 norm counts to only include Aldh1l1 samples
  777. #Updated clean_counts
  778. Aldh1l1_clean_counts <- counts_clean_filtered_rows_9_2_8_1[,c(14:16, 29:31)]
  780. #Updated traits
  781. Aldh1l1_traits <- traits_9_2_8_1[c(14:16, 29:31),]
  783. rownames(Aldh1l1_traits)
  784. # [1] "21047FL-134-01-41" "21047FL-134-01-42" "21047FL-134-01-43" "23238R-04-08" "23238R-04-09"
  785. # [6] "23238R-04-10"
  787. colnames(Aldh1l1_clean_counts)
  788. # [1] "21047FL-134-01-41" "21047FL-134-01-42" "21047FL-134-01-43" "23238R-04-08" "23238R-04-09"
  789. # [6] "23238R-04-10"
  792. #create DESeq data set (dds) matrix for the 3 Aldh1l1 and 3 Camk2a pos.pulldowns
  793. Aldh1l1_dds <- DESeqDataSetFromMatrix(countData = Aldh1l1_clean_counts, #filtered, unnorm counts is the input here
  794. colData = Aldh1l1_traits,
  795. design = ~Group)
  796. dim(Aldh1l1_dds)
  797. # [1] 17054 6
  799. # set the factor level
  800. Aldh1l1_dds$Group <- relevel(Aldh1l1_dds$Group, ref = "Aldh1l1.RiboTag.pos.pulldown")
  802. #Run DESeq
  803. Aldh1l1_dds<- DESeq(Aldh1l1_dds)
  806. # Create PCA plots
  807. pdf(file = "1d-1. Exp. 9-2 and 8-1_RNA-seq_Aldh1l1 pulldown groups_DESeq2_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  808. # perform DESeq2-basedPCA on rlog transformed data of 16 samples w/o neg pulldowns
  810. rldr <- rlog(Aldh1l1_dds, blind = TRUE)
  811. plotPCA(rldr, intgroup = "Group") # plots based on normalized counts
  813. pcaData <- plotPCA(rldr, intgroup = "Group", returnData=TRUE)
  814. percentVar <- round(100 * attr(pcaData, "percentVar"))
  815. pca4 <- ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  816. geom_point(size=3) +
  817. xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  818. ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  819. ggtitle("DESeq2 PCA plot for log2(mean mRNA counts)") +
  820. theme(panel.background = element_rect(fill = "transparent", color = "black"),
  821. legend.background = element_rect(fill = "transparent"),
  822. legend.key = element_rect(fill = "transparent"),
  823. plot.title = element_text(hjust = 0.5)) + # Center the title
  824. labs(color = "Transcriptome groups") +
  825. coord_fixed()
  827. print(pca4)
  829. dev.off()
  831. ################################################################################
  832. # perform limma:plotMDS-based PCA of 18 samples w/o neg pulldowns
  833. pdf(file = "1c-2. Exp. 9-2 and 8-1_RNA-seq_Aldh1l1 groups_Limma_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  834. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  835. heights = c(1.3,0.7,1.3), # Heights of the rows
  836. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  838. Grouping_vector_Aldh1l1 <- Aldh1l1_traits$Group
  840. Aldh1l1_norm_counts <- normalized_counts_2[,c(4:11, 14:21)]
  842. pch_values <- ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.bulk.cortex", 16,
  843. ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.RiboTag.bulk.cortex", 16,
  844. ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.pos.pulldown", 16,
  845. ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.LPS.pos.pulldown", 16, NA))))
  848. pt_colors <- ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.bulk.cortex", "coral",
  849. ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.RiboTag.bulk.cortex", "green3",
  850. ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.pos.pulldown", "purple",
  851. ifelse(Grouping_vector_Aldh1l1== "Aldh1l1.LPS.pos.pulldown", "darkturquoise", NA))))
  853. legend_groups <- c("Aldh1l1.bulk.cortex", "Aldh1l1.RiboTag.bulk.cortex", "Aldh1l1.pos.pulldown","Aldh1l1.LPS.pos.pulldown")
  855. plotMDS_Aldh1l1<- plotMDS(log2(Aldh1l1_norm_counts), #top = 500,
  856. labels = NULL, pch = pch_values, col = pt_colors,
  857. cex = 3, dim.plot = c(1,2), gene.selection = "pairwise",
  858. xlab = NULL, ylab = NULL, plot = TRUE, var.explained = TRUE)
  859. mtext(side=3, text="MDS Plot for log2(mean mRNA counts)", cex=1, line=2)
  861. plot.new()
  862. legend("topright", legend = legend_groups, col = c("coral", "green3", "purple", "darkturquoise"),
  863. pch = 16, title = "Transcriptome groups", cex = 1.4)
  865. dev.off()
  867. ###############################################################################
  868. # 1d. Aldh1l1.pos.pulldowns vs RiboTag.bulk.cortexs
  869. # subset DESeq2 norm counts to only include Aldh1l1 and RiboTag.bulk.cortex samples
  871. #Updated clean_counts
  872. Aldh1l1_Vs_RiboTag.bulk.cortex_clean_counts <- counts_clean_filtered_rows_9_2_8_1[,c(14:16, 22:23)]
  874. #Updated traits
  875. Aldh1l1_Vs_RiboTag.bulk.cortex_traits <- traits_9_2_8_1[c(14:16, 22:23),]
  876. # Sample.ID.. Group
  877. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  878. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  879. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  880. # 23238R-04-01 CR1 RiboTag.bulk.cortex
  881. # 23238R-04-02 CR2 RiboTag.bulk.cortex
  883. colnames(Aldh1l1_Vs_RiboTag.bulk.cortex_clean_counts)
  884. # [1] "21047FL-134-01-41" "21047FL-134-01-42" "21047FL-134-01-43" "23238R-04-01" "23238R-04-02"
  886. #create DESeq data set (dds) matrix for the 3 Aldh1l1 and 3 Camk2a pos.pulldowns
  887. Aldh1l1_Vs_RiboTag.bulk.cortex_dds <- DESeqDataSetFromMatrix(countData = Aldh1l1_Vs_RiboTag.bulk.cortex_clean_counts, #filtered, unnorm counts is the input here
  888. colData = Aldh1l1_Vs_RiboTag.bulk.cortex_traits,
  889. design = ~Group)
  890. dim(Aldh1l1_Vs_RiboTag.bulk.cortex_dds)
  891. # [1] 16825 6
  893. # set the factor level
  894. Aldh1l1_Vs_RiboTag.bulk.cortex_dds$Group <- relevel(Aldh1l1_Vs_RiboTag.bulk.cortex_dds$Group, ref = "RiboTag.bulk.cortex")
  896. #Run DESeq
  897. Aldh1l1_Vs_RiboTag.bulk.cortex_dds<- DESeq(Aldh1l1_Vs_RiboTag.bulk.cortex_dds)
  900. # Create PCA plots
  901. pdf(file = "1d-1. Exp. 9-2 and 8-1_RNA-seq_Aldh1l1 vs RiboTag.bulk.cortexs_DESeq2_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  902. # perform DESeq2-basedPCA on rlog transformed data of 16 samples w/o neg pulldowns
  904. rldr <- rlog(Aldh1l1_Vs_RiboTag.bulk.cortex_dds, blind = TRUE)
  905. plotPCA(rldr, intgroup = "Group") # plots based on normalized counts
  907. pcaData <- plotPCA(rldr, intgroup = "Group", returnData=TRUE)
  908. percentVar <- round(100 * attr(pcaData, "percentVar"))
  909. pca3 <- ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  910. geom_point(size=3) +
  911. xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  912. ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  913. theme(panel.background = element_rect(fill = "transparent", color = "black"),
  914. legend.background = element_rect(fill = "transparent"),
  915. legend.key = element_rect(fill = "transparent")) +
  916. labs(color = "Transcriptome groups") +
  917. coord_fixed()
  919. dev.off()
  921. # perform limma:plotMDS-based PCA of 18 samples w/o neg pulldowns
  922. pdf(file = "1d-2. Exp. 9-2 and 8-1_RNA-seq_Aldh1l1 vs Camk2a_Limma_PCA_CR_06182024.pdf", width = 8, height = 8.5)
  923. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  924. heights = c(1.3,0.7,1.3), # Heights of the rows
  925. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  927. Grouping_vector_Aldh1l1VsCamk2a <- Aldh1l1VsCamk2a_traits$Group
  929. Aldh1l1VsCamk2a_norm_counts <- normalized_counts_2[,c(12:14, 20:22)]
  931. colnames(Aldh1l1VsCamk2a_norm_counts)
  932. # [1] "21047FL-134-01-41" "21047FL-134-01-42" "21047FL-134-01-43" "22082R-12-01" "22082R-12-02"
  933. # [6] "22082R-12-03"
  935. pch_values <- ifelse(Grouping_vector_Aldh1l1VsCamk2a== "Aldh1l1.pos.pulldown", 16,
  936. ifelse(Grouping_vector_Aldh1l1VsCamk2a== "RiboTag.bulk.cortex", 16, NA))
  939. pt_colors <- ifelse(Grouping_vector_Aldh1l1VsCamk2a== "Aldh1l1.pos.pulldown", "purple",
  940. ifelse(Grouping_vector_Aldh1l1VsCamk2a== "RiboTag.bulk.cortex", "dodgerblue", NA))
  942. legend_groups <- c("Aldh1l1.pos.pulldown", "RiboTag.bulk.cortex")
  944. plotMDS_Aldh1l1<- plotMDS(log2(Aldh1l1VsCamk2a_norm_counts), #top = 500,
  945. labels = NULL, pch = pch_values, col = pt_colors,
  946. cex = 3, dim.plot = c(1,2), gene.selection = "pairwise",
  947. xlab = NULL, ylab = NULL, plot = TRUE, var.explained = TRUE)
  948. mtext(side=3, text="MDS Plot for log2(mean mRNA counts)", cex=1, line=2)
  950. plot.new()
  951. legend("topright", legend = legend_groups, col = c("purple", "dodgerblue"),
  952. pch = 16, title = "Transcriptome groups", cex = 1.4)
  954. dev.off()
  957. ################################################################################
  958. ################################################################################
  959. # Step 2. Correlations of mean mRNA counts
  960. # a. Aldh1l1.bulk.cortex (n=3/group) vs Aldh1l1.RiboTag.bulk.cortex (n=2/group)
  961. # b. Aldh1l1.RiboTag.bulk.cortex (n=2/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  962. # c. Aldh1l1.pos.pulldown (n=3/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  964. # check for finite values
  965. any(is.finite(normalized_counts_2))
  966. # [1] TRUE zeroes are considered finite
  968. dim(normalized_counts_2)
  969. # [1] 17054 27 (4 neg pulldowns removed)
  971. # CR: Because there are zeroes in the DESeq2 normalized counts data matrix,
  972. # the log2 transformation changes these values to infinite values (-ifn).
  973. # A correlation coefficient cannot be calculated with infinite values, so
  974. # I am going to add a "pseudocount" of + 1 to all the counts values to
  975. # avoid taking a log of 0.
  977. norm_counts_pseudo <- normalized_counts_2 + 1
  978. any(is.finite(norm_counts_pseudo))
  980. norm_counts_log <- log2(norm_counts_pseudo)
  981. any(is.finite(norm_counts_log))
  983. dim(norm_counts_log)
  984. # [1] 17054 27 (4 neg pulldowns removed)
  986. # load library to generate plot
  987. library(ggplot2)
  989. ################################################################################
  990. # 2a. Aldh1l1.bulk.cortex (n=3/group) vs Aldh1l1.RiboTag.bulk.cortex (n=2/group)
  992. pdf(file = "2a-c. Exp. 9-2 and 8-1_RNA-seq_counts correlations_CR_06182024.pdf", width = 6, height = 6, family = "Arial")
  993. par(mfrow= c(2,1)) # centers plot on single page (I think)
  995. # subset into groups of interest and take the row mean
  996. Aldh1l1.bulk.cortex <- rowMeans(norm_counts_log[,c(4:6)])
  997. Aldh1l1.RiboTag.bulk.cortex <- rowMeans(norm_counts_log[,c(22:24)])
  999. # generate a linear correlation
  1000. correlation_2a <- cor(Aldh1l1.bulk.cortex, Aldh1l1.RiboTag.bulk.cortex)
  1001. correlation_2a
  1002. # [1] 0.9788304
  1004. # convert data to a data frame
  1005. df_2a <- data.frame(Aldh1l1.bulk.cortex = Aldh1l1.bulk.cortex, Aldh1l1.RiboTag.bulk.cortex = Aldh1l1.RiboTag.bulk.cortex)
  1007. # create correlation plot
  1008. correlation_plot_2a <- ggplot(df_2a, aes(x = Aldh1l1.bulk.cortex, y = Aldh1l1.RiboTag.bulk.cortex)) +
  1009. geom_point(color = "black", fill = "grey", pch = 21, alpha = 0.5, size = 3) +
  1010. labs(x = "Aldh1l1.bulk.cortex log2(mean mRNA counts)",
  1011. y = "Aldh1l1.RiboTag.bulk.cortex log2(mean mRNA counts)") +
  1012. ggtitle("2a. Aldh1l1.bulk.cortex vs Aldh1l1.RiboTag.bulk.cortex") +
  1013. theme_minimal() +
  1014. theme(
  1015. panel.background = element_blank(), # Remove panel background
  1016. panel.grid.major = element_blank(), # Remove major grid lines
  1017. panel.grid.minor = element_blank(), # Remove minor grid lines
  1018. axis.line = element_line(color = "black"), # Add black x and y axes
  1019. axis.ticks = element_line(color = "black"), # Add ticks to the axes
  1020. plot.title = element_text(hjust = 0.5) # Center the title
  1021. ) +
  1022. annotate("text", x = max(Aldh1l1.bulk.cortex), y = min(Aldh1l1.RiboTag.bulk.cortex),
  1023. label = paste("Correlation coefficient:", round(correlation_2a, 2)),
  1024. hjust = 1, vjust = 0, color = "black") +
  1025. annotate("text", x = min(df_2a$Aldh1l1.bulk.cortex), y = min(df_2a$Aldh1l1.RiboTag.bulk.cortex) - 2,
  1026. label = paste("\nNumber of genes:", nrow(df_2a)),
  1027. hjust = 0, vjust = 1, color = "black") + # Adjusted y coordinate
  1028. scale_x_continuous(breaks = seq(0, 20, by = 5), limits = c(0, 20)) +
  1029. scale_y_continuous(breaks = seq(0, 20, by = 5), limits = c(0, 20))
  1031. print(correlation_plot_2a)
  1033. ################################################################################
  1034. # 2b. Aldh1l1.RiboTag.bulk.cortex (n=2/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  1035. # subset into groups of interest and take the row mean
  1036. Aldh1l1.RiboTag.bulk.cortex <- rowMeans(norm_counts_log[,c(22:24)])
  1037. Aldh1l1.RiboTag.pos.pulldown <- rowMeans(norm_counts_log[,c(25:27)])
  1039. # convert data to a data frame
  1040. df_2b <- data.frame(Aldh1l1.RiboTag.bulk.cortex = Aldh1l1.RiboTag.bulk.cortex, Aldh1l1.RiboTag.pos.pulldown = Aldh1l1.RiboTag.pos.pulldown)
  1042. # calculate correlation coefficient
  1043. correlation_2b <- cor(Aldh1l1.RiboTag.bulk.cortex, Aldh1l1.RiboTag.pos.pulldown) # cor() is based on Pearson's correlation by default, can change if needed by adding method = "spearman" or others
  1044. correlation_2b
  1045. # [1] 0.9169184
  1047. correlation_plot_2b <- ggplot(df_2b, aes(x = Aldh1l1.RiboTag.bulk.cortex, y = Aldh1l1.RiboTag.pos.pulldown)) +
  1048. geom_point(color = "black", fill = "red4", pch = 21, alpha = 0.5, size = 3) +
  1049. labs(x = "Aldh1l1.RiboTag.bulk.cortex log2(mean mRNA counts)",
  1050. y = "Aldh1l1.RiboTag.pos.pulldown log2(mean mRNA counts)") +
  1051. ggtitle("2b. Aldh1l1.Ribotag bulk.cortex vs pulldown") +
  1052. theme_minimal() +
  1053. theme(
  1054. panel.background = element_blank(), # Remove panel background
  1055. panel.grid.major = element_blank(), # Remove major grid lines
  1056. panel.grid.minor = element_blank(), # Remove minor grid lines
  1057. axis.line = element_line(color = "black"), # Add black x and y axes
  1058. axis.ticks = element_line(color = "black"), # Add ticks to the axes
  1059. plot.title = element_text(hjust = 0.5) # Center the title
  1060. ) +
  1061. annotate("text", x = max(Aldh1l1.RiboTag.bulk.cortex), y = min(Aldh1l1.RiboTag.pos.pulldown),
  1062. label = paste("Correlation coefficient:", round(correlation_2b, 2)),
  1063. hjust = 1, vjust = 0, color = "black") +
  1064. annotate("text", x = min(df_2b$Aldh1l1.RiboTag.bulk.cortex), y = min(df_2b$Aldh1l1.RiboTag.pos.pulldown) - 2,
  1065. label = paste("Number of genes:", nrow(df_2a)), # # of genes are the same for each df
  1066. hjust = 0, vjust = 1, color = "black") + # Adjusted y coordinate
  1067. scale_x_continuous(breaks = seq(0, 20, by = 5), limits = c(0, 20)) +
  1068. scale_y_continuous(breaks = seq(0, 20, by = 5), limits = c(0, 20))
  1070. print(correlation_plot_2b)
  1072. ################################################################################
  1073. # 2c.Aldh1l1.pos.pulldown (n=3/group) vs Aldh1l1.RiboTag.pos.pulldown (n=3/group)
  1074. # subset into groups of interest and take the row mean
  1075. Aldh1l1.pos.pulldown <- rowMeans(norm_counts_log[,c(12:14)])
  1076. Aldh1l1.RiboTag.pos.pulldown <- rowMeans(norm_counts_log[,c(25:27)])
  1078. # convert data to a data frame
  1079. df_2c <- data.frame(Aldh1l1.pos.pulldown = Aldh1l1.pos.pulldown, Aldh1l1.RiboTag.pos.pulldown = Aldh1l1.RiboTag.pos.pulldown)
  1081. # calculate correlation coefficient
  1082. correlation_2c <- cor(Aldh1l1.pos.pulldown, Aldh1l1.RiboTag.pos.pulldown) # cor() is based on Pearson's correlation by default, can change if needed by adding method = "spearman" or others
  1083. correlation_2c
  1084. # [1] 0.9326028
  1086. correlation_plot_2c <- ggplot(df_2c, aes(x = Aldh1l1.pos.pulldown, y = Aldh1l1.RiboTag.pos.pulldown)) +
  1087. geom_point(color = "black", fill = "red2", pch = 21, alpha = 0.5, size = 3) +
  1088. labs(x = "Aldh1l1.pos.pulldown log2(mean mRNA counts)",
  1089. y = "Aldh1l1.RiboTag log2(mean mRNA counts)") +
  1090. ggtitle("2c. Aldh1l1 TurboID vs RiboTag pulldown RNA") +
  1091. theme_minimal() +
  1092. theme(
  1093. panel.background = element_blank(), # Remove panel background
  1094. panel.grid.major = element_blank(), # Remove major grid lines
  1095. panel.grid.minor = element_blank(), # Remove minor grid lines
  1096. axis.line = element_line(color = "black"), # Add black x and y axes
  1097. axis.ticks = element_line(color = "black"), # Add ticks to the axes
  1098. plot.title = element_text(hjust = 0.5) # Center the title
  1099. ) +
  1100. annotate("text", x = max(Aldh1l1.pos.pulldown), y = min(Aldh1l1.RiboTag.pos.pulldown),
  1101. label = paste("Correlation coefficient:", round(correlation_2c, 2)),
  1102. hjust = 1, vjust = 0, color = "black") +
  1103. annotate("text", x = min(df_2c$Aldh1l1.pos.pulldown), y = min(df_2c$Aldh1l1.RiboTag.pos.pulldown) - 2,
  1104. label = paste("Number of genes:", nrow(df_2a)), # # of genes are the same for each df
  1105. hjust = 1, vjust = 0, color = "black") + # Adjusted y coordinate
  1106. scale_x_continuous(breaks = seq(0, 20, by = 5), limits = c(0, 20)) +
  1107. scale_y_continuous(breaks = seq(0, 20, by = 5), limits = c(0, 20))
  1109. print(correlation_plot_2c)
  1111. dev.off()
  1114. ################################################################################
  1115. ################################################################################
  1116. # Step 3 Part I: Diffex volcanoes
  1118. # C2: Astrocyte-RiboTag pulldown vs bulk.cortex
  1119. # C4: Astrocyte-RiboTag vs astrocyte-TurboID pulldown
  1120. ################################################################################
  1121. #load Libraries
  1122. library(DESeq2)
  1123. library(tidyverse)
  1124. library(biomaRt)
  1125. library(pheatmap)
  1126. library(ggplot2)
  1127. library(ggrepel)
  1129. ################################################################################
  1130. # C2: Astrocyte-RiboTag pulldown vs bulk.cortex
  1131. pdf(file = "C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_06182024.pdf", height = 8, width = 8.5)
  1132. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  1133. heights = c(1.3,0.7,1.3), # Heights of the rows
  1134. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  1136. final_clean_counts <- counts_clean_filtered_rows_9_2_8_1
  1137. dim(final_clean_counts)
  1138. # [1] 17054 31
  1140. #Updated clean_counts for comparison 1
  1141. C2_clean_counts <- final_clean_counts[,c(24:26, 29:31)]
  1143. #Updated traits_9_2_8_1 for comparison 1
  1144. C2_traits_9_2_8_1 <- traits_9_2_8_1[c(24:26, 29:31),]
  1145. C2_traits_9_2_8_1
  1146. # Sample.ID.. Group
  1147. # 23238R-04-03 CR3 Aldh1l1.RiboTag.bulk.cortex
  1148. # 23238R-04-04 CR4 Aldh1l1.RiboTag.bulk.cortex
  1149. # 23238R-04-05 CR5 Aldh1l1.RiboTag.bulk.cortex
  1150. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1151. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1152. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1154. #Deseq2 data matrix is needed for downstream analysis
  1155. dds2 <- DESeqDataSetFromMatrix(countData = C2_clean_counts,
  1156. colData = C2_traits_9_2_8_1,
  1157. design = ~Group)
  1158. dim(dds2)
  1159. # [1] 17054 6
  1161. dds2
  1163. # set the factor level
  1164. dds2$Group <- relevel(dds2$Group, ref = "Aldh1l1.RiboTag.bulk.cortex")
  1166. #Run DESeq
  1167. dds2 <- DESeq(dds2)
  1169. #diff p adj value
  1170. res05 <- results(dds2, alpha = 0.05)
  1171. res05 #table w statistics
  1173. #diff ex analysis
  1174. summary(res05)
  1175. # out of 17051 with nonzero bulk.cortex read count
  1176. # adjusted p-value < 0.05
  1177. # LFC > 0 (up) : 5895, 35%
  1178. # LFC < 0 (down) : 5347, 31%
  1179. # outliers [1] : 0, 0%
  1180. # low counts [2] : 0, 0%
  1181. # (mean count < 0)
  1183. dim(res05)
  1184. # [1] 17054 6
  1186. #contrasts
  1187. resultsNames(dds2)
  1188. # [1] "Intercept"
  1189. # [2] "Group_Aldh1l1.RiboTag.pos.pulldown_vs_Aldh1l1.RiboTag.bulk.cortex"
  1191. #plot MA
  1192. #In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.
  1193. plotMA(res05, main = "Comparison 2: Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex")
  1195. #use the log transform on the data set
  1196. #The above will plot the most varied genes across all conditions
  1197. rld <- rlog(dds2, blind=TRUE)
  1198. de2 <- as.data.frame(res05)
  1199. de2$padj[is.na(de2$padj)] <- 1
  1201. write.csv(de2, file = "C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_06182024.csv")
  1203. var <- order(de2$padj, decreasing=F)
  1204. #topAdjpGenes <- which(var <= 20)
  1205. #allows you to hold your enrs ids and sample names in a df
  1206. #Most sig genes! Top 50
  1207. matrix <- assay(rld)[ var[1:50], ]
  1208. matrix <- matrix - rowMeans(matrix)
  1210. matrixx <- matrix
  1211. #thresh_val <- min(abs(max(matrix)),abs(min(matrix)))
  1212. thresh_val <- 0.5
  1213. matrixx[matrixx < -thresh_val] <- -thresh_val
  1214. matrixx[matrixx > thresh_val] <- thresh_val
  1216. annotation_data <- as.data.frame(colData(rld)[c("Group")])
  1217. pheatmap(matrixx, annotation_col=annotation_data, show_rownames = T,show_colnames = T, cluster_cols=TRUE, main = "Top 50 padj C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054")
  1219. #up and down list created to pull out diff ex genes
  1220. up <- res05[which(res05$log2FoldChange > 1 & res05$padj < .05),]
  1221. down <- res05[which(res05$log2FoldChange < -1 & res05$padj < .05),]
  1223. write.csv(up, file = "C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_Up_06182024.csv")
  1224. write.csv(down, file = "C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_Down_06182024.csv")
  1226. #The counts here need to match with the one counter above by the model!
  1227. up_genes <- read.csv("C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_Up_06182024.csv", header = TRUE, row.names = 1)
  1228. down_genes <- read.csv("C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_Down_06182024.csv", header = TRUE, row.names = 1)
  1230. #bulk.cortex df of up and down genes based on a criteria that were extracted!
  1231. up_down_genes <- rbind(up_genes,down_genes)
  1233. dim(up_down_genes)
  1234. # [1] 5512 6
  1236. #volcano plot of diffex with padj df
  1237. de2 <- read.csv("C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_06182024.csv")
  1239. de2$diffexpressed <- "NS"
  1240. de2$diffexpressed[de2$log2FoldChange > 1.0 & de2$padj < 0.05] <- "UP"
  1241. de2$diffexpressed[de2$log2FoldChange < -1.0 & de2$padj < 0.05] <- "DOWN"
  1243. de2$delabel <- NA
  1244. de2$delabel[de2$diffexpressed != "NS"] <- de2$X[de2$diffexpressed != "NS"]
  1246. ggplot(data=de2, aes(x=log2FoldChange, y=-log10(padj), col=diffexpressed, label=delabel)) +
  1247. geom_point() +
  1248. ggtitle(label = "C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054") +
  1249. theme_minimal() +
  1250. geom_text_repel(segment.colour = NA, max.overlaps = getOption("ggrepel.max.overlaps", default = 15)) +
  1251. scale_color_manual(values=c("black", "grey", "red4")) +
  1252. geom_vline(xintercept=c(-1.0, 1.0), linetype = "longdash", col="black") +
  1253. geom_hline(yintercept=-log10(0.05), linetype = "longdash", col="black") +
  1254. theme(legend.position="none") +
  1255. scale_y_continuous(name = "-log10(padj)",
  1256. limits = c(-5, 300)) +
  1257. scale_x_continuous(name = "log2(difference)
  1258. Aldh1l1.RiboTag.pulldown - Aldh1l1.RiboTag bulk.cortex RNA",
  1259. limits = c(-10, 10))
  1261. #Perform PCA on rlog transformed data
  1262. rldr <- rlog(dds2, blind = TRUE)
  1263. plotPCA(rldr, intgroup = "Group")
  1265. dev.off()
  1267. ################################################################################
  1268. ################################################################################
  1269. # C4: Astrocyte-RiboTag pulldown vs Astrocyte-TurboID pulldown
  1270. #use base pdf function to extract plots into a single pdf
  1271. pdf(file = "C4_Aldh1l1.TurboID vs Aldh1l1.RiboTag pulldowns_rowmeans_17054_Diffex_06182024.pdf", height = 8, width = 8.5)
  1272. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  1273. heights = c(1.3,0.7,1.3), # Heights of the rows
  1274. widths = c(0.65,0.65,0.7)) # Widths of the columns -- the distance to squash the second plot to the left, (because we do not duplicate legends)
  1276. final_clean_counts <- counts_clean_filtered_rows_9_2_8_1
  1277. dim(final_clean_counts)
  1278. # [1] 17054 31
  1280. #Updated clean_counts for comparison 1
  1281. C4_clean_counts <- final_clean_counts[,c(14:16, 29:31)]
  1283. #Updated traits_9_2_8_1 for comparison 1
  1284. C4_traits_9_2_8_1 <- traits_9_2_8_1[c(14:16, 29:31),]
  1285. C4_traits_9_2_8_1
  1286. # Sample.ID.. Group
  1287. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  1288. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  1289. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  1290. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1291. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1292. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1294. #Deseq2 data matrix is needed for downstream analysis
  1295. dds4 <- DESeqDataSetFromMatrix(countData = C4_clean_counts,
  1296. colData = C4_traits_9_2_8_1,
  1297. design = ~Group)
  1298. dim(dds4)
  1299. # [1] 17054 6
  1301. dds4
  1303. # set the factor level
  1304. dds4$Group <- relevel(dds4$Group, ref = "Aldh1l1.RiboTag.pos.pulldown")
  1306. #Run DESeq
  1307. dds4 <- DESeq(dds4)
  1309. #diff p adj value
  1310. res05 <- results(dds4, alpha = 0.05)
  1311. res05 #table w statistics
  1313. #diff ex analysis
  1314. summary(res05)
  1315. # out of 17054 with nonzero bulk.cortex read count
  1316. # adjusted p-value < 0.05
  1317. # LFC > 0 (up) : 4519, 26%
  1318. # LFC < 0 (down) : 4555, 27%
  1319. # outliers [1] : 0, 0%
  1320. # low counts [2] : 0, 0%
  1321. # (mean count < 1)
  1323. dim(res05)
  1324. # [1] 17054 6
  1326. #contrasts
  1327. resultsNames(dds4)
  1328. # [1] "Intercept"
  1329. # [2] "Group_Aldh1l1.pos.pulldown_vs_Aldh1l1.RiboTag.pos.pulldown"
  1331. #plot MA
  1332. #In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.
  1333. plotMA(res05, main = "Comparison 4: Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown")
  1335. #use the log transform on the data set
  1336. #The above will plot the most varied genes across all conditions
  1337. rld <- rlog(dds4, blind=TRUE)
  1338. de4 <- as.data.frame(res05)
  1339. de4$padj[is.na(de4$padj)] <- 1
  1341. write.csv(de4, file = "C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_06182024.csv")
  1343. var <- order(de4$padj, decreasing=F)
  1344. #topAdjpGenes <- which(var <= 20)
  1345. #allows you to hold your enrs ids and sample names in a df
  1346. #Most sig genes! Top 50
  1347. matrix <- assay(rld)[ var[1:50], ]
  1348. matrix <- matrix - rowMeans(matrix)
  1350. matrixx <- matrix
  1351. #thresh_val <- min(abs(max(matrix)),abs(min(matrix)))
  1352. thresh_val <- 0.5
  1353. matrixx[matrixx < -thresh_val] <- -thresh_val
  1354. matrixx[matrixx > thresh_val] <- thresh_val
  1356. annotation_data <- as.data.frame(colData(rld)[c("Group")])
  1357. pheatmap(matrixx, annotation_col=annotation_data, show_rownames = T,show_colnames = T, cluster_cols=TRUE, main = "Top 50 padj C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054")
  1359. #up and down list created to pull out diff ex genes
  1360. up <- res05[which(res05$log2FoldChange > 1 & res05$padj < .05),]
  1361. down <- res05[which(res05$log2FoldChange < -1 & res05$padj < .05),]
  1363. write.csv(up, file = "C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_Up_06182024.csv")
  1364. write.csv(down, file = "C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_Down_06182024.csv")
  1366. #The counts here need to match with the one counter above by the model!
  1367. up_genes <- read.csv("C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_Up_06182024.csv", header = TRUE, row.names = 1)
  1368. down_genes <- read.csv("C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_Down_06182024.csv", header = TRUE, row.names = 1)
  1370. #bulk.cortex df of up and down genes based on a criteria that were extracted!
  1371. up_down_genes <- rbind(up_genes,down_genes)
  1373. dim(up_down_genes)
  1374. # [1] 4362 6
  1376. #volcano plot of diffex with padj df
  1377. de4 <- read.csv("C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_06182024.csv")
  1379. de4$diffexpressed <- "NS"
  1380. de4$diffexpressed[de4$log2FoldChange > 1.0 & de4$padj < 0.05] <- "UP"
  1381. de4$diffexpressed[de4$log2FoldChange < -1.0 & de4$padj < 0.05] <- "DOWN"
  1383. de4$delabel <- NA
  1384. de4$delabel[de4$diffexpressed != "NS"] <- de4$X[de4$diffexpressed != "NS"]
  1386. ggplot(data=de4, aes(x=log2FoldChange, y=-log10(padj), col=diffexpressed, label=delabel)) +
  1387. geom_point() +
  1388. ggtitle(label = "C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054") +
  1389. theme_minimal() +
  1390. geom_text_repel(segment.colour = NA, max.overlaps = getOption("ggrepel.max.overlaps", default = 15)) +
  1391. scale_color_manual(values=c("red4", "grey", "purple4")) +
  1392. geom_vline(xintercept=c(-1.0, 1.0), linetype = "longdash", col="black") +
  1393. geom_hline(yintercept=-log10(0.05), linetype = "longdash", col="black") +
  1394. theme(legend.position="none") +
  1395. scale_y_continuous(name = "-log10(padj)",
  1396. limits = c(-5, 200)) +
  1397. scale_x_continuous(name = "log2(difference)
  1398. Aldh1l1.TurboID - Aldh1l1.RiboTag pulldown RNA",
  1399. limits = c(-30, 30))
  1401. #Perform PCA on rlog transformed data
  1402. rldr <- rlog(dds1, blind = TRUE)
  1403. plotPCA(rldr, intgroup = "Group")
  1405. dev.off()
  1407. ################################################################################
  1408. ################################################################################
  1409. # Step 3 Part II: Diffex GSEA
  1410. ################################################################################
  1411. # Gene set erichment analysis with Eric Dammer's GOparallel pipeline
  1412. # https://github.com/edammer/GOparallel
  1413. # performs hypergeometric overlap with Fisher Exact Test for enrichment p<0.05 and 5 minimum genes per ontology
  1414. # currently outputs files that can be used as GO-Elite input, not read by this script...
  1415. # outputs Z score barplots for each input list (module, or DiffEx Protein list) across 6 ontology types in the larger GMT files available from the Bader Lab website.
  1417. # Parameters for GOparallel GSEA on volcano-defined DEG lists
  1418. source("GOparallel-FET.R")
  1420. ## Parameters saved to variables in memory.
  1421. filePath <- getwd()
  1422. #Folder that (may) contain the input file specified above, and which will contain the outFilename project Folder.
  1423. modulesInMemory=FALSE
  1424. ANOVAgroups=TRUE #change to FALSE if using inputFile
  1425. #if true, modulesInMemory ignored. Volcano pipeline code should already have been run!
  1426. #inputFile will be ignored
  1427. maxBarsPerOntology=6
  1429. panelDimensions=c(2,2) #dimensions of the individual parblots within a page of the main barplot PDF output
  1430. pageDimensions=c(8.5,11) #main barplot PDF output page dimensions, in inches
  1431. parallelThreads=7 # need at least 2 threads to run
  1432. removeRedundantGOterms=TRUE
  1433. #if true, the 3 GO ontology types are collapased into a minimal set of less redundant terms using the below OBO file
  1434. cocluster=FALSE
  1435. #If TRUE, output PDF of signed Zscore coclustering on GO:cellular component terms (useful for WGCNA modules)
  1437. # "Clean up" to run before GOPar
  1438. rm(ANOVAout)
  1440. ###############################################################################
  1441. # Relevant diffex comparisons:
  1442. # C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_06182024.csv
  1443. # C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_06182024.csv
  1445. ###############################################################################
  1446. # Comparsion 2: Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex
  1447. # organize DESeq2 diffex output file to mirror LFQ-MS ANOVAout file
  1448. C2_ANOVA<-read.csv("C2_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_06182024.csv", header=TRUE, row.names=1)
  1449. C2_ANOVAout<-C2_ANOVA[,c(4,1,5,2)]
  1450. colnames(C2_ANOVAout)[c(3,4)]<-c("Aldh1l1.RiboTag.pull-Aldh1l1.RiboTag.bulk.cortex","diff Aldh1l1.RiboTag.pull.Aldh1l1.RiboTag.bulk.cortex")
  1451. write.csv(C2_ANOVAout, file = "C2_RNA_ANOVAout_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_rowmeans_17054_Diffex_06182024.csv")
  1453. # GOpar
  1454. flip=c(0)
  1455. outFilename <- "C2_RNA_Aldh1l1.RiboTag.pulldown vs Aldh1l1.RiboTag.bulk.cortex_pval_DEGs-GO"
  1456. GMTdatabaseFile="Mouse_GO_AllPathways_with_GO_iea_November_07_2023_symbol.gmt"
  1457. GOparallel(C2_ANOVAout) # parameters are set in global environment as above; if not set, the function falls back to defaults and looks for all inputs available.
  1459. ###############################################################################
  1460. # Comparison 4: Astrocyte-TurboID pulldown vs Astrocyte-RiboTag pulldown
  1462. # organize DESeq2 diffex output file to mirror LFQ-MS ANOVAout file
  1463. C4_ANOVA<-read.csv("C4_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_06182024.csv", header=TRUE, row.names=1)
  1464. C4_ANOVAout<-C4_ANOVA[,c(4,1,5,2)]
  1465. colnames(C4_ANOVAout)[c(3,4)]<-c("Aldh1l1.TurboID-Aldh1l1.RiboTag","diff Aldh1l1.TurboID.Aldh1l1.Ribotag")
  1466. write.csv(C4_ANOVAout, file = "C4_RNA_ANOVAout_Aldh1l1.pos.pulldown vs Aldh1l1.RiboTag.pos.pulldown_rowmeans_17054_Diffex_06182024.csv")
  1468. # GOpar
  1469. flip=c(0)
  1470. outFilename <- "C4_RNA_Aldh1l1.TurboID vs Aldh1l1.RiboTag pulldowns_pval_DEGs-GO"
  1471. GMTdatabaseFile="Mouse_GO_AllPathways_with_GO_iea_November_07_2023_symbol.gmt"
  1472. GOparallel(C4_ANOVAout) # parameters are set in global environment as above; if not set, the function falls back to defaults and looks for all inputs available.
  1475. ####################################################
  1476. # heatmap of bulk vs pulldown RNA
  1477. #Load necessary libraries
  1478. library(pheatmap)
  1479. library(dplyr)
  1480. library(grDevices) # for colorRampPalette
  1482. # Load top 50 cell-type gene list from Barres lab
  1483. top50_cell_type_genes_Sloan <- read.csv("mouse_cell_markers_Sloan_03062023.csv", header = TRUE)
  1485. # Create a mapping of genes to cell types
  1486. gene_to_cell_type <- stack(top50_cell_type_genes_Sloan)
  1487. colnames(gene_to_cell_type) <- c("Gene", "CellType")
  1489. # Load expression data
  1490. final_clean_counts <- read.csv("1b.Exp. 9-2 and 8-1_RNA-seq_clean_counts_diffex filter-17054x31_CR_06182024.csv", row.names = 1, check.names=FALSE)
  1492. dim(final_clean_counts)
  1493. # [1] 17054 31
  1495. # Load traits data
  1496. traits_9_2_8_1 <- read.csv("Exp. 9-2 and 8-1 RNA traits.csv", row.names = 1)
  1497. dim(traits_9_2_8_1)
  1498. # [1] 31 2
  1500. ################################################################################
  1501. # Comparison 2: Aldh1l1.RiboTag.pos.pulldown vs Aldh1l1.RiboTag.bulk.cortex
  1502. # subset counts data
  1504. # Subset counts data
  1505. Aldh1l1.RiboTag_pullVsBulk_clean_counts <- final_clean_counts[,c(24:26, 29:31)]
  1507. # Updated traits_9_2_8_1
  1508. Aldh1l1.RiboTag_pullVsBulk_traits <- traits_9_2_8_1[c(24:26, 29:31),]
  1509. Aldh1l1.RiboTag_pullVsBulk_traits
  1510. # Sample.ID.. Group
  1511. # 23238R-04-03 CR3 Aldh1l1.RiboTag.bulk.cortex
  1512. # 23238R-04-04 CR4 Aldh1l1.RiboTag.bulk.cortex
  1513. # 23238R-04-05 CR5 Aldh1l1.RiboTag.bulk.cortex
  1514. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1515. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1516. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1518. # Check if sample names match
  1519. if(!all(row.names(Aldh1l1.RiboTag_pullVsBulk_traits) %in% colnames(Aldh1l1.RiboTag_pullVsBulk_clean_counts))) {
  1520. stop("Sample names in traits file do not match column names in expression data.")
  1521. }
  1523. group_info <- Aldh1l1.RiboTag_pullVsBulk_traits$Group
  1524. names(group_info) <- row.names(Aldh1l1.RiboTag_pullVsBulk_traits)
  1526. annotation_col <- data.frame(Group = group_info)
  1527. rownames(annotation_col) <- names(group_info)
  1529. # Filter for astrocyte genes
  1530. astrocyte_genes <- gene_to_cell_type %>% filter(CellType == "Astrocyte") %>% pull(Gene)
  1532. # Filter counts for astrocyte genes
  1533. Aldh1l1.RiboTag_pullVsBulk_filtered_data <- Aldh1l1.RiboTag_pullVsBulk_clean_counts[rownames(Aldh1l1.RiboTag_pullVsBulk_clean_counts) %in% astrocyte_genes, ]
  1535. # Calculate row z-scores
  1536. row_z_scores <- t(apply(Aldh1l1.RiboTag_pullVsBulk_filtered_data, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1537. colnames(row_z_scores) <- colnames(Aldh1l1.RiboTag_pullVsBulk_filtered_data)
  1539. # Adjust row names to include cell type information
  1540. new_row_names <- paste(rownames(row_z_scores), "Astrocyte", sep = " | ")
  1541. rownames(row_z_scores) <- new_row_names
  1543. pdf(file = "C2_Exp. 8-1_RNA-seq_Aldh1l1.RiboTag_pullVsBulk_Zhang_heatmap_07022024.pdf", width = 8, height = 8)
  1544. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  1545. heights = c(1,1,1), # Equal heights for all rows
  1546. widths = c(1,1,1)) # Equal widths for all columns
  1548. # Define a custom color palette
  1549. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1551. # Create the heatmap with annotation and custom colors
  1552. RiboTag_astrocyte_heatmap <- pheatmap(row_z_scores,
  1553. cluster_rows = FALSE,
  1554. cluster_cols = FALSE,
  1555. annotation_col = annotation_col,
  1556. show_rownames = TRUE,
  1557. show_colnames = TRUE,
  1558. color = custom_colors)
  1559. ###############################################################################
  1560. # Filter for neuronal genes
  1561. neuronal_genes <- gene_to_cell_type %>% filter(CellType == "Neuron") %>% pull(Gene)
  1563. # Filter counts for astrocyte genes
  1564. Aldh1l1.RiboTag_pullVsBulk_filtered_data <- Aldh1l1.RiboTag_pullVsBulk_clean_counts[rownames(Aldh1l1.RiboTag_pullVsBulk_clean_counts) %in% neuronal_genes, ]
  1566. # Calculate row z-scores
  1567. row_z_scores <- t(apply(Aldh1l1.RiboTag_pullVsBulk_filtered_data, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1568. colnames(row_z_scores) <- colnames(Aldh1l1.RiboTag_pullVsBulk_filtered_data)
  1570. # Adjust row names to include cell type information
  1571. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  1572. rownames(row_z_scores) <- new_row_names
  1574. # Define a custom color palette
  1575. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1577. # Create the heatmap with annotation and custom colors
  1578. RiboTag_neuron_heatmap <- pheatmap(row_z_scores,
  1579. cluster_rows = FALSE,
  1580. cluster_cols = FALSE,
  1581. annotation_col = annotation_col,
  1582. show_rownames = TRUE,
  1583. show_colnames = TRUE,
  1584. color = custom_colors)
  1585. dev.off()
  1586. ################################################################################
  1587. ################################################################################
  1588. # Comparison 4: Astrocyte-TurboID pulldown vs Astrocyte-RiboTag pulldown
  1590. # Subset counts data
  1591. Aldh1l1_TurboVsRibo_clean_counts <- final_clean_counts[,c(4:6, 29:31)]
  1593. # Updated traits_9_2_8_1
  1594. Aldh1l1_TurboVsRibo_traits <- traits_9_2_8_1[c(14:16, 29:31),]
  1595. # Sample.ID.. Group
  1596. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  1597. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  1598. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  1599. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1600. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1601. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1603. # Check if sample names match
  1604. if(!all(row.names(Aldh1l1_TurboVsRibo_traits) %in% colnames(Aldh1l1_TurboVsRibo_clean_counts))) {
  1605. stop("Sample names in traits file do not match column names in expression data.")
  1606. }
  1608. group_info <- Aldh1l1_TurboVsRibo_traits$Group
  1609. names(group_info) <- row.names(Aldh1l1_TurboVsRibo_traits)
  1611. annotation_col <- data.frame(Group = group_info)
  1612. rownames(annotation_col) <- names(group_info)
  1614. # Filter for astrocyte genes
  1615. astrocyte_genes <- gene_to_cell_type %>% filter(CellType == "Astrocyte") %>% pull(Gene)
  1617. # Filter counts for astrocyte genes
  1618. Aldh1l1_TurboVsRibo_filtered_data <- Aldh1l1_TurboVsRibo_clean_counts[rownames(Aldh1l1_TurboVsRibo_clean_counts) %in% astrocyte_genes, ]
  1620. # Calculate row z-scores
  1621. row_z_scores <- t(apply(Aldh1l1_TurboVsRibo_filtered_data, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1622. colnames(row_z_scores) <- colnames(Aldh1l1_TurboVsRibo_filtered_data)
  1624. # Adjust row names to include cell type information
  1625. new_row_names <- paste(rownames(row_z_scores), "Astrocyte", sep = " | ")
  1626. rownames(row_z_scores) <- new_row_names
  1628. pdf(file = "C4_Exp. 8-1_RNA-seq_Aldh1l1_TurboVsRibo_pulls_Zhang_heatmap_07022024.pdf", width = 8, height = 8)
  1629. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  1630. heights = c(1,1,1), # Equal heights for all rows
  1631. widths = c(1,1,1)) # Equal widths for all columns
  1633. # Define a custom color palette
  1634. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1636. # Create the heatmap with annotation and custom colors
  1637. Aldh1l1_pulls_astrocyte_heatmap <- pheatmap(row_z_scores,
  1638. cluster_rows = FALSE,
  1639. cluster_cols = FALSE,
  1640. annotation_col = annotation_col,
  1641. show_rownames = TRUE,
  1642. show_colnames = TRUE,
  1643. color = custom_colors)
  1644. ###############################################################################
  1645. # Filter for neuronal genes
  1646. neuronal_genes <- gene_to_cell_type %>% filter(CellType == "Neuron") %>% pull(Gene)
  1648. # Filter counts for astrocyte genes
  1649. Aldh1l1_TurboVsRibo_filtered_data <- Aldh1l1_TurboVsRibo_clean_counts[rownames(Aldh1l1_TurboVsRibo_clean_counts) %in% neuronal_genes, ]
  1651. # Calculate row z-scores
  1652. row_z_scores <- t(apply(Aldh1l1_TurboVsRibo_filtered_data, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1653. colnames(row_z_scores) <- colnames(Aldh1l1_TurboVsRibo_filtered_data)
  1655. # Adjust row names to include cell type information
  1656. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  1657. rownames(row_z_scores) <- new_row_names
  1659. # Define a custom color palette
  1660. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1662. # Create the heatmap with annotation and custom colors
  1663. Aldh1l1_pulls_neuron_heatmap <- pheatmap(row_z_scores,
  1664. cluster_rows = FALSE,
  1665. cluster_cols = FALSE,
  1666. annotation_col = annotation_col,
  1667. show_rownames = TRUE,
  1668. show_colnames = TRUE,
  1669. color = custom_colors)
  1670. dev.off()
  1672. ################################################################################
  1673. ##################################################################################
  1674. # Filtering for cell-type markers from Zhang et al., 2014 and Sharma et al., 2015
  1676. # load Zhang and Sharma cell-type marker union list from Johnson et al. 2022 Nat Neurosci
  1677. Zhang_Sharma_markers <- read.csv("Johnson et al.. 2022 Nat Neurosci Sharma and Zhang marker union list_mouse_from Eric.csv", header = TRUE)
  1679. # Create a mapping of genes to cell types
  1680. gene_to_cell_type_2 <- stack(Zhang_Sharma_markers)
  1681. colnames(gene_to_cell_type_2) <- c("Gene", "CellType")
  1683. ################################################################################
  1684. # Comparison 2: Aldh1l1.RiboTag.pos.pulldown vs Aldh1l1.RiboTag.bulk.cortex
  1685. # subset counts data
  1687. # Subset counts data
  1688. Aldh1l1.RiboTag_pullVsBulk_clean_counts <- final_clean_counts[,c(24:26, 29:31)]
  1690. # Updated traits_9_2_8_1
  1691. Aldh1l1.RiboTag_pullVsBulk_traits <- traits_9_2_8_1[c(24:26, 29:31),]
  1692. Aldh1l1.RiboTag_pullVsBulk_traits
  1693. # Sample.ID.. Group
  1694. # 23238R-04-03 CR3 Aldh1l1.RiboTag.bulk.cortex
  1695. # 23238R-04-04 CR4 Aldh1l1.RiboTag.bulk.cortex
  1696. # 23238R-04-05 CR5 Aldh1l1.RiboTag.bulk.cortex
  1697. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1698. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1699. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1701. # Check if sample names match
  1702. if(!all(row.names(Aldh1l1.RiboTag_pullVsBulk_traits) %in% colnames(Aldh1l1.RiboTag_pullVsBulk_clean_counts))) {
  1703. stop("Sample names in traits file do not match column names in expression data.")
  1704. }
  1706. group_info <- Aldh1l1.RiboTag_pullVsBulk_traits$Group
  1707. names(group_info) <- row.names(Aldh1l1.RiboTag_pullVsBulk_traits)
  1709. annotation_col <- data.frame(Group = group_info)
  1710. rownames(annotation_col) <- names(group_info)
  1712. # Filter for astrocyte genes
  1713. astrocyte_genes_2 <- gene_to_cell_type_2 %>% filter(CellType == "Astrocytes") %>% pull(Gene)
  1715. # Filter counts for astrocyte genes
  1716. Aldh1l1.RiboTag_pullVsBulk_filtered_data_2 <- Aldh1l1.RiboTag_pullVsBulk_clean_counts[rownames(Aldh1l1.RiboTag_pullVsBulk_clean_counts) %in% astrocyte_genes_2, ]
  1718. # Calculate row z-scores
  1719. row_z_scores <- t(apply(Aldh1l1.RiboTag_pullVsBulk_filtered_data_2, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1720. colnames(row_z_scores) <- colnames(Aldh1l1.RiboTag_pullVsBulk_filtered_data_2)
  1722. # Adjust row names to include cell type information
  1723. new_row_names <- paste(rownames(row_z_scores), "Astrocyte", sep = " | ")
  1724. rownames(row_z_scores) <- new_row_names
  1726. sum(is.na(row_z_scores))
  1727. # [1] 12
  1728. sum(is.nan(row_z_scores))
  1729. # [1] 12
  1730. sum(is.infinite(row_z_scores))
  1731. # [1] 0
  1733. # Identify rows with NA, NaN values
  1734. rows_with_na <- apply(row_z_scores, 1, function(x) any(is.na(x)))
  1735. rows_with_nan <- apply(row_z_scores, 1, function(x) any(is.nan(x)))
  1737. # Print the indices of rows with NA values
  1738. which(rows_with_na)
  1739. # Eif2s3y Kdm5d
  1740. # 320 525
  1741. which(rows_with_nan)
  1742. # Eif2s3y Kdm5d
  1743. # 320 525
  1745. #rows_to_remove <- c(320, 525)
  1746. #row_z_scores_clean <- row_z_scores[-rows_to_remove, ]
  1747. row_z_scores_clean <- na.omit(row_z_scores)
  1748. rows_with_na <- apply(row_z_scores_clean, 1, function(x) any(is.na(x)))
  1749. which(rows_with_na)
  1750. # named integer(0)
  1752. pdf(file = "C2_Exp. 8-1_RNA-seq_Aldh1l1.RiboTag_pullVsBulk_Zhang and Sharma_heatmaps_07112024.pdf", width = 8, height = 8)
  1753. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  1754. heights = c(1,1,1), # Equal heights for all rows
  1755. widths = c(1,1,1)) # Equal widths for all columns
  1757. # Define a custom color palette
  1758. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1760. # Create the heatmap with annotation and custom colors
  1761. RiboTag_astrocyte_heatmap_2 <- pheatmap(row_z_scores_clean,
  1762. cluster_rows = TRUE,
  1763. cluster_cols = TRUE,
  1764. annotation_col = annotation_col,
  1765. show_rownames = FALSE,
  1766. show_colnames = TRUE,
  1767. color = custom_colors,
  1768. main = "Top astrocytic markers from Zhang and Sharma")
  1769. ###############################################################################
  1770. # Filter for neuronal genes
  1771. neuronal_genes_2 <- gene_to_cell_type_2 %>% filter(CellType == "Neurons") %>% pull(Gene)
  1773. # Filter counts for astrocyte genes
  1774. Aldh1l1.RiboTag_pullVsBulk_filtered_data_2 <- Aldh1l1.RiboTag_pullVsBulk_clean_counts[rownames(Aldh1l1.RiboTag_pullVsBulk_clean_counts) %in% neuronal_genes_2, ]
  1776. # Calculate row z-scores
  1777. row_z_scores <- t(apply(Aldh1l1.RiboTag_pullVsBulk_filtered_data_2, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1778. colnames(row_z_scores) <- colnames(Aldh1l1.RiboTag_pullVsBulk_filtered_data_2)
  1780. # Adjust row names to include cell type information
  1781. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  1782. rownames(row_z_scores) <- new_row_names
  1784. # Define a custom color palette
  1785. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1787. # Create the heatmap with annotation and custom colors
  1788. RiboTag_neuron_heatmap_2 <- pheatmap(row_z_scores,
  1789. cluster_rows = TRUE,
  1790. cluster_cols = TRUE,
  1791. annotation_col = annotation_col,
  1792. show_rownames = FALSE,
  1793. show_colnames = TRUE,
  1794. color = custom_colors,
  1795. main = "Top neuronal markers from Zhang and Sharma")
  1796. dev.off()
  1798. ################################################################################
  1799. ###############################################################################
  1800. # Filter for neuronal genes
  1801. neuronal_genes <- gene_to_cell_type %>% filter(CellType == "Neuron") %>% pull(Gene)
  1803. # Filter counts for astrocyte genes
  1804. TurboVsRibo_bulk.cortexs_filtered_data <- TurboVsRibo_bulk.cortexs_clean_counts[rownames(TurboVsRibo_bulk.cortexs_clean_counts) %in% neuronal_genes, ]
  1806. # Calculate row z-scores
  1807. row_z_scores <- t(apply(TurboVsRibo_bulk.cortexs_filtered_data, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1808. colnames(row_z_scores) <- colnames(TurboVsRibo_bulk.cortexs_filtered_data)
  1810. # Adjust row names to include cell type information
  1811. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  1812. rownames(row_z_scores) <- new_row_names
  1814. # Define a custom color palette
  1815. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1817. # Create the heatmap with annotation and custom colors
  1818. Aldh1l1_bulk.cortexs_neuron_heatmap <- pheatmap(row_z_scores,
  1819. cluster_rows = FALSE,
  1820. cluster_cols = FALSE,
  1821. annotation_col = annotation_col,
  1822. show_rownames = TRUE,
  1823. show_colnames = TRUE,
  1824. color = custom_colors)
  1825. dev.off()
  1827. # Subset counts data
  1828. Aldh1l1_TurboVsRibo_clean_counts <- final_clean_counts[,c(14:16, 29:31)]
  1830. # Updated traits_9_2_8_1
  1831. Aldh1l1_TurboVsRibo_traits <- traits_9_2_8_1[c(14:16, 29:31),]
  1832. Aldh1l1_TurboVsRibo_traits
  1833. # Sample.ID.. Group
  1834. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  1835. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  1836. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  1837. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1838. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1839. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1841. # Check if sample names match
  1842. if(!all(row.names(Aldh1l1_TurboVsRibo_traits) %in% colnames(Aldh1l1_TurboVsRibo_clean_counts))) {
  1843. stop("Sample names in traits file do not match column names in expression data.")
  1844. }
  1846. group_info <- Aldh1l1_TurboVsRibo_traits$Group
  1847. names(group_info) <- row.names(Aldh1l1_TurboVsRibo_traits)
  1849. annotation_col <- data.frame(Group = group_info)
  1850. rownames(annotation_col) <- names(group_info)
  1852. # Filter for astrocyte genes
  1853. astrocyte_genes_2 <- gene_to_cell_type_2 %>% filter(CellType == "Astrocytes") %>% pull(Gene)
  1855. # Filter counts for astrocyte genes
  1856. Aldh1l1_TurboVsRibo_filtered_data_2 <- Aldh1l1_TurboVsRibo_clean_counts[rownames(Aldh1l1_TurboVsRibo_clean_counts) %in% astrocyte_genes_2, ]
  1858. # Calculate row z-scores
  1859. row_z_scores <- t(apply(Aldh1l1_TurboVsRibo_filtered_data_2, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1860. colnames(row_z_scores) <- colnames(Aldh1l1_TurboVsRibo_filtered_data_2)
  1862. # Adjust row names to include cell type information
  1863. new_row_names <- paste(rownames(row_z_scores), "Astrocyte", sep = " | ")
  1864. rownames(row_z_scores) <- new_row_names
  1866. pdf(file = "C4_Exp. 8-1_RNA-seq_Aldh1l1_TurboVsRibo_pulls_Zhang and Sharma_heatmaps_07112024.pdf", width = 8, height = 8)
  1867. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  1868. heights = c(1,1,1), # Equal heights for all rows
  1869. widths = c(1,1,1)) # Equal widths for all columns
  1871. # Define a custom color palette
  1872. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1874. # Create the heatmap with annotation and custom colors
  1875. Aldh1l1_pulls_astrocyte_heatmap_2 <- pheatmap(row_z_scores,
  1876. cluster_rows = TRUE,
  1877. cluster_cols = TRUE,
  1878. annotation_col = annotation_col,
  1879. show_rownames = FALSE,
  1880. show_colnames = TRUE,
  1881. color = custom_colors,
  1882. main = "Top astrocytic markers from Zhang and Sharma")
  1883. ###############################################################################
  1884. # Filter for neuronal genes
  1885. neuronal_genes_2 <- gene_to_cell_type_2 %>% filter(CellType == "Neurons") %>% pull(Gene)
  1887. # Filter counts for astrocyte genes
  1888. Aldh1l1_TurboVsRibo_filtered_data_2 <- Aldh1l1_TurboVsRibo_clean_counts[rownames(Aldh1l1_TurboVsRibo_clean_counts) %in% neuronal_genes_2, ]
  1890. # Calculate row z-scores
  1891. row_z_scores <- t(apply(Aldh1l1_TurboVsRibo_filtered_data_2, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1892. colnames(row_z_scores) <- colnames(Aldh1l1_TurboVsRibo_filtered_data_2)
  1894. # Adjust row names to include cell type information
  1895. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  1896. rownames(row_z_scores) <- new_row_names
  1898. # Define a custom color palette
  1899. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  1901. # Create the heatmap with annotation and custom colors
  1902. Aldh1l1_pulls_neuron_heatmap_2 <- pheatmap(row_z_scores,
  1903. cluster_rows = TRUE,
  1904. cluster_cols = TRUE,
  1905. annotation_col = annotation_col,
  1906. show_rownames = FALSE,
  1907. show_colnames = TRUE,
  1908. color = custom_colors,
  1909. main = "Top neuronal markers from Zhang and Sharma")
  1910. dev.off()
  1912. ################################################################################
  1913. #################################################################################
  1914. # load Seyfried et al., 2017 Zhang only list
  1915. Zhang_markers <- read.csv("Seyfried et al., 2017 Cell Syst S4 Table (Zhang list).csv", header = TRUE)
  1916. gene_to_cell_type_3 <- stack(Zhang_markers)
  1917. colnames(gene_to_cell_type_3) <- c("Gene", "CellType")
  1919. # Load expression data
  1920. final_clean_counts <- read.csv("1b.Exp. 9-2 and 8-1_RNA-seq_clean_counts_diffex filter-17054x31_CR_06182024.csv", row.names = 1, check.names=FALSE)
  1922. dim(final_clean_counts)
  1923. # [1] 17054 31
  1925. # Load traits data
  1926. traits_9_2_8_1 <- read.csv("Exp. 9-2 and 8-1 RNA traits.csv", row.names = 1)
  1927. dim(traits_9_2_8_1)
  1928. # [1] 31 2
  1930. ##################################################################################
  1931. ##################################################################################
  1932. # Comparison 2: Aldh1l1.RiboTag.pos.pulldown vs Aldh1l1.RiboTag.bulk.cortex
  1934. # Subset counts data
  1935. Aldh1l1.RiboTag_pullVsBulk_clean_counts <- final_clean_counts[,c(24:26, 29:31)]
  1937. # Updated traits_9_2_8_1
  1938. Aldh1l1.RiboTag_pullVsBulk_traits <- traits_9_2_8_1[c(24:26, 29:31),]
  1939. Aldh1l1.RiboTag_pullVsBulk_traits
  1940. # Sample.ID.. Group
  1941. # 23238R-04-03 CR3 Aldh1l1.RiboTag.bulk.cortex
  1942. # 23238R-04-04 CR4 Aldh1l1.RiboTag.bulk.cortex
  1943. # 23238R-04-05 CR5 Aldh1l1.RiboTag.bulk.cortex
  1944. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  1945. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  1946. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  1948. # Check if sample names match
  1949. if(!all(row.names(Aldh1l1.RiboTag_pullVsBulk_traits) %in% colnames(Aldh1l1.RiboTag_pullVsBulk_clean_counts))) {
  1950. stop("Sample names in traits file do not match column names in expression data.")
  1951. }
  1953. group_info <- Aldh1l1.RiboTag_pullVsBulk_traits$Group
  1954. names(group_info) <- row.names(Aldh1l1.RiboTag_pullVsBulk_traits)
  1956. annotation_col <- data.frame(Group = group_info)
  1957. rownames(annotation_col) <- names(group_info)
  1959. # Filter for Astrocytes genes
  1960. Astrocyte_genes_3 <- gene_to_cell_type_3 %>% filter(CellType == "Astrocytes") %>% pull(Gene)
  1962. # Filter counts for Astrocytes genes
  1963. Aldh1l1.RiboTag_pullVsBulk_filtered_data_3 <- Aldh1l1.RiboTag_pullVsBulk_clean_counts[rownames(Aldh1l1.RiboTag_pullVsBulk_clean_counts) %in% Astrocyte_genes_3, ]
  1965. # Calculate row z-scores
  1966. row_z_scores <- t(apply(Aldh1l1.RiboTag_pullVsBulk_filtered_data_3, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  1967. colnames(row_z_scores) <- colnames(Aldh1l1.RiboTag_pullVsBulk_filtered_data_3)
  1969. # Adjust row names to include cell type information
  1970. new_row_names <- paste(rownames(row_z_scores), "Astrocytes", sep = " | ")
  1971. rownames(row_z_scores) <- new_row_names
  1973. sum(is.na(row_z_scores))
  1974. # [1] 12
  1975. sum(is.nan(row_z_scores))
  1976. # [1] 12
  1977. sum(is.infinite(row_z_scores))
  1978. # [1] 0
  1980. # Identify rows with NA, NaN values
  1981. rows_with_na <- apply(row_z_scores, 1, function(x) any(is.na(x)))
  1982. rows_with_nan <- apply(row_z_scores, 1, function(x) any(is.nan(x)))
  1984. # Print the indices of rows with NA values
  1985. which(rows_with_na)
  1986. # Eif2s3y Kdm5d
  1987. # 320 525
  1988. which(rows_with_nan)
  1989. # Eif2s3y Kdm5d
  1990. # 320 525
  1992. #rows_to_remove <- c(320, 525)
  1993. #row_z_scores_clean <- row_z_scores[-rows_to_remove, ]
  1994. row_z_scores_clean <- na.omit(row_z_scores)
  1995. rows_with_na <- apply(row_z_scores_clean, 1, function(x) any(is.na(x)))
  1996. which(rows_with_na)
  1997. # named integer(0)
  1999. pdf(file = "C2_Exp. 8-1_RNA-seq_Aldh1l1.RiboTag_pullVsBulk_Zhang all_heatmap_07122024.pdf", width = 8, height = 8)
  2000. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  2001. heights = c(1,1,1), # Equal heights for all rows
  2002. widths = c(1,1,1)) # Equal widths for all columns
  2004. # Define a custom color palette
  2005. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  2007. # Create the heatmap with annotation and custom colors
  2008. RiboTag_Astrocytes_heatmap_3 <- pheatmap(row_z_scores_clean,
  2009. cluster_rows = TRUE,
  2010. cluster_cols = TRUE,
  2011. annotation_col = annotation_col,
  2012. show_rownames = FALSE,
  2013. show_colnames = TRUE,
  2014. color = custom_colors,
  2015. main = "Top astrocytic markers from Zhang")
  2016. ###############################################################################
  2017. # Filter for neuronal genes
  2018. neuronal_genes_3 <- gene_to_cell_type_3 %>% filter(CellType == "Neurons") %>% pull(Gene)
  2020. # Filter counts for Astrocytes genes
  2021. Aldh1l1.RiboTag_pullVsBulk_filtered_data_3 <- Aldh1l1.RiboTag_pullVsBulk_clean_counts[rownames(Aldh1l1.RiboTag_pullVsBulk_clean_counts) %in% neuronal_genes_3, ]
  2023. # Calculate row z-scores
  2024. row_z_scores <- t(apply(Aldh1l1.RiboTag_pullVsBulk_filtered_data_3, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  2025. colnames(row_z_scores) <- colnames(Aldh1l1.RiboTag_pullVsBulk_filtered_data_3)
  2027. # Adjust row names to include cell type information
  2028. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  2029. rownames(row_z_scores) <- new_row_names
  2031. # Define a custom color palette
  2032. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  2034. # Create the heatmap with annotation and custom colors
  2035. RiboTag_neuron_heatmap_3 <- pheatmap(row_z_scores,
  2036. cluster_rows = TRUE,
  2037. cluster_cols = TRUE,
  2038. annotation_col = annotation_col,
  2039. show_rownames = FALSE,
  2040. show_colnames = TRUE,
  2041. color = custom_colors,
  2042. main = "Top neuronal markers from Zhang")
  2043. dev.off()
  2045. ################################################################################
  2046. # Comparison 4: Astrocyte-TurboID pulldown vs Astrocyte-RiboTag pulldown
  2048. # Subset counts data
  2049. Aldh1l1_TurboVsRibo_clean_counts <- final_clean_counts[,c(14:16, 29:31)]
  2051. # Updated traits_9_2_8_1
  2052. Aldh1l1_TurboVsRibo_traits <- traits_9_2_8_1[c(14:16, 29:31),]
  2053. Aldh1l1_TurboVsRibo_traits
  2054. # Sample.ID.. Group
  2055. # 21047FL-134-01-41 15-RNA 30L Aldh1l1.pos.pulldown
  2056. # 21047FL-134-01-42 16-RNA 31B Aldh1l1.pos.pulldown
  2057. # 21047FL-134-01-43 17-RNA 24L Aldh1l1.pos.pulldown
  2058. # 23238R-04-08 CR8 Aldh1l1.RiboTag.pos.pulldown
  2059. # 23238R-04-09 CR9 Aldh1l1.RiboTag.pos.pulldown
  2060. # 23238R-04-10 CR10 Aldh1l1.RiboTag.pos.pulldown
  2062. # Check if sample names match
  2063. if(!all(row.names(Aldh1l1_TurboVsRibo_traits) %in% colnames(Aldh1l1_TurboVsRibo_clean_counts))) {
  2064. stop("Sample names in traits file do not match column names in expression data.")
  2065. }
  2067. group_info <- Aldh1l1_TurboVsRibo_traits$Group
  2068. names(group_info) <- row.names(Aldh1l1_TurboVsRibo_traits)
  2070. annotation_col <- data.frame(Group = group_info)
  2071. rownames(annotation_col) <- names(group_info)
  2073. # Filter for Astrocytes genes
  2074. Astrocyte_genes_3 <- gene_to_cell_type_3 %>% filter(CellType == "Astrocytes") %>% pull(Gene)
  2076. # Filter counts for Astrocytes genes
  2077. Aldh1l1_TurboVsRibo_filtered_data_3 <- Aldh1l1_TurboVsRibo_clean_counts[rownames(Aldh1l1_TurboVsRibo_clean_counts) %in% Astrocyte_genes_3, ]
  2079. # Calculate row z-scores
  2080. row_z_scores <- t(apply(Aldh1l1_TurboVsRibo_filtered_data_3, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  2081. colnames(row_z_scores) <- colnames(Aldh1l1_TurboVsRibo_filtered_data_3)
  2083. # Adjust row names to include cell type information
  2084. new_row_names <- paste(rownames(row_z_scores), "Astrocytes", sep = " | ")
  2085. rownames(row_z_scores) <- new_row_names
  2087. pdf(file = "C4_Exp. 8-1_RNA-seq_Aldh1l1_TurboVsRibo_pulls_Zhang all_heatmaps_07122024.pdf", width = 8, height = 8)
  2088. layout(matrix(c(1,1,2, 1,1,3, 4,5,6), nrow = 3, ncol = 3, byrow=TRUE),
  2089. heights = c(1,1,1), # Equal heights for all rows
  2090. widths = c(1,1,1)) # Equal widths for all columns
  2092. # Define a custom color palette
  2093. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  2095. # Create the heatmap with annotation and custom colors
  2096. Aldh1l1_pulls_Astrocytes_heatmap_3 <- pheatmap(row_z_scores,
  2097. cluster_rows = TRUE,
  2098. cluster_cols = TRUE,
  2099. annotation_col = annotation_col,
  2100. show_rownames = FALSE,
  2101. show_colnames = TRUE,
  2102. color = custom_colors,
  2103. main = "Top astrocytic markers from Zhang")
  2104. ###############################################################################
  2105. # Filter for neuronal genes
  2106. neuronal_genes_3 <- gene_to_cell_type_3 %>% filter(CellType == "Neurons") %>% pull(Gene)
  2108. # Filter counts for Astrocytes genes
  2109. Aldh1l1_TurboVsRibo_filtered_data_3 <- Aldh1l1_TurboVsRibo_clean_counts[rownames(Aldh1l1_TurboVsRibo_clean_counts) %in% neuronal_genes_3, ]
  2111. # Calculate row z-scores
  2112. row_z_scores <- t(apply(Aldh1l1_TurboVsRibo_filtered_data_3, 1, function(x) scale(x, center = TRUE, scale = TRUE)))
  2113. colnames(row_z_scores) <- colnames(Aldh1l1_TurboVsRibo_filtered_data_3)
  2115. # Adjust row names to include cell type information
  2116. new_row_names <- paste(rownames(row_z_scores), "Neuron", sep = " | ")
  2117. rownames(row_z_scores) <- new_row_names
  2119. # Define a custom color palette
  2120. custom_colors <- colorRampPalette(c("lightblue", "white", "darkblue"))(100)
  2122. # Create the heatmap with annotation and custom colors
  2123. Aldh1l1_pulls_neuron_heatmap_3 <- pheatmap(row_z_scores,
  2124. cluster_rows = TRUE,
  2125. cluster_cols = TRUE,
  2126. annotation_col = annotation_col,
  2127. show_rownames = FALSE,
  2128. show_colnames = TRUE,
  2129. color = custom_colors,
  2130. main = "Top neuronal markers from Zhang")
  2131. dev.off()
  2133. ################################################################################
  2134. # end of data analysis and visualization

Astrocyte-RiboTag_mRNA seq analysis.R at commit 4303bf0, under MIT · at the source

Overview

Authors: Christina C Ramelow1,2,3, Eric B Dammer2,4, Hailian Xiao1,2, Lihong Cheng1,5, Prateek Kumar6, Claudia Espinosa-Garcia6, Maureen M Sampson3, Dilpreet Kour6, Ruth S Nelson6, Sneha Malepati6, Rashmi Kumari6, Wooyoung Eric Jang6, Qi Guo2,4, Pritha Bagchi7, Duc M Duong2,7, Nicholas T Seyfried1,2,4,7, Steven A Sloan3, Srikant Rangaraju1,5
ORCID iDs: Christina C Ramelow, Eric B Dammer, Prateek Kumar, Claudia Espinosa-Garcia, Maureen M Sampson, Dilpreet Kour, Sneha Malepati, Wooyoung Eric Jang, Qi Guo, Pritha Bagchi, Nicholas T Seyfried, Steven A Sloan, Srikant Rangaraju
  1. Department of Neurology, Emory University, Atlanta, GA USA
  2. Center for Neurodegenerative Disease, Emory University, Atlanta, GA USA
  3. Department of Human Genetics, Emory University, Atlanta, GA USA
  4. Department of Biochemistry, Emory University, Atlanta, GA USA
  5. Department of Pharmacology and Chemical Biology, Emory University, Atlanta, GA USA
  6. Department of Neurology, Yale University, New Haven, CT USA
  7. Emory Integrated Proteomics Core, Emory University, Atlanta, GA USA
Institutions: Emory University (United States); Yale University (United States)
Journal: Nature communications, volume 17, issue 1, article 4636
Dates: received 31 January 2025; accepted 10 March 2026; published online 31 March 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-71098-4 · PMID 41917013 · PMCID PMC13199526 · OpenAlex W7143899533
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: genetics / omics (modality), human (organism), mouse (organism), cellular / molecular (subfield)
Methods: Statistics, Evoked potentials
Keywords: Molecular neuroscience, Gene expression analysis, Proteomic analysis, Cellular neuroscience, Bioinformatics
MeSH: Gene Expression Profiling*, Neurons*, Proteome*, Proteomics*, Transcriptome*, Animals, Astrocytes, HEK293 Cells, Humans, Lipopolysaccharides, Mice, Microglia, Multiomics, RNA, Messenger (* major topic)
Topic: Biotin and Related Studies (Cell Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: U.S. Department of Health &amp; Human Services | NIH | National Institute on Aging (R01AG075820, 1F31AG079597-01A1); NIA NIH HHS (R01 AG071587, R01 AG075820); U.S. Department of Health &amp; Human Services | NIH | National Institute of Mental Health (NIMH) (R01MH125956)
Citations: cited by 3 papers (Europe PMC); 155 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 5 matches between paragraphs and lines of code.

cramelow/SPARO

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 4303bf04063009a22d82e7b276f59d38f83692dc, 5 March 2026
Languages: R (7)
Size: 8 files, 7 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: license file
Not found: README, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (7 files), limma (7 files), tidyverse (7 files), DESeq2 (4 files), pheatmap (4 files), WGCNA (3 files), data.table (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
8 files
At the source: github.com/cramelow/SPARO

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:

Read it in the paper: doi.org/10.1038/s41467-026-71098-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;
  • 7 scripts, each with its path and the digest of its content;
  • 5 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/s41467-026-71098-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, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 18 authors, 5 keywords, 14 MeSH terms, 3 funders, 154 references.

Cite

This paper

Ramelow, C. C., Dammer, E. B., Xiao, H., Cheng, L., Kumar, P., Espinosa-Garcia, C., Sampson, M. M., Kour, D., Nelson, R. S., Malepati, S., Kumari, R., Jang, W. E., Guo, Q., Bagchi, P., Duong, D. M., Seyfried, N. T., Sloan, S. A., & Rangaraju, S. (2026). Simultaneous profiling of native-state proteomes and transcriptomes of neural cell types using proximity labeling. Nature communications, 17(1), 4636. https://doi.org/10.1038/s41467-026-71098-4

BibTeX

@article{ramelow2026simultaneous,
author = {Ramelow, Christina C and Dammer, Eric B and Xiao, Hailian and Cheng, Lihong and Kumar, Prateek and Espinosa-Garcia, Claudia and Sampson, Maureen M and Kour, Dilpreet and Nelson, Ruth S and Malepati, Sneha and Kumari, Rashmi and Jang, Wooyoung Eric and Guo, Qi and Bagchi, Pritha and Duong, Duc M and Seyfried, Nicholas T and Sloan, Steven A and Rangaraju, Srikant},
title = {{Simultaneous profiling of native-state proteomes and transcriptomes of neural cell types using proximity labeling}},
journal = {Nature communications},
year = {2026},
month = mar,
volume = {17},
number = {1},
pages = {4636},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-71098-4},
url = {https://doi.org/10.1038/s41467-026-71098-4},
pmid = {41917013},
pmcid = {PMC13199526}
}

RIS

TY - JOUR
AU - Ramelow, Christina C
AU - Dammer, Eric B
AU - Xiao, Hailian
AU - Cheng, Lihong
AU - Kumar, Prateek
AU - Espinosa-Garcia, Claudia
AU - Sampson, Maureen M
AU - Kour, Dilpreet
AU - Nelson, Ruth S
AU - Malepati, Sneha
AU - Kumari, Rashmi
AU - Jang, Wooyoung Eric
AU - Guo, Qi
AU - Bagchi, Pritha
AU - Duong, Duc M
AU - Seyfried, Nicholas T
AU - Sloan, Steven A
AU - Rangaraju, Srikant
TI - Simultaneous profiling of native-state proteomes and transcriptomes of neural cell types using proximity labeling
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/03/31
VL - 17
IS - 1
SP - 4636
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-71098-4
UR - https://doi.org/10.1038/s41467-026-71098-4
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-71098-4",
"type": "article-journal",
"title": "Simultaneous profiling of native-state proteomes and transcriptomes of neural cell types using proximity labeling",
"container-title": "Nature communications",
"author": [
{
"family": "Ramelow",
"given": "Christina C"
},
{
"family": "Dammer",
"given": "Eric B"
},
{
"family": "Xiao",
"given": "Hailian"
},
{
"family": "Cheng",
"given": "Lihong"
},
{
"family": "Kumar",
"given": "Prateek"
},
{
"family": "Espinosa-Garcia",
"given": "Claudia"
},
{
"family": "Sampson",
"given": "Maureen M"
},
{
"family": "Kour",
"given": "Dilpreet"
},
{
"family": "Nelson",
"given": "Ruth S"
},
{
"family": "Malepati",
"given": "Sneha"
},
{
"family": "Kumari",
"given": "Rashmi"
},
{
"family": "Jang",
"given": "Wooyoung Eric"
},
{
"family": "Guo",
"given": "Qi"
},
{
"family": "Bagchi",
"given": "Pritha"
},
{
"family": "Duong",
"given": "Duc M"
},
{
"family": "Seyfried",
"given": "Nicholas T"
},
{
"family": "Sloan",
"given": "Steven A"
},
{
"family": "Rangaraju",
"given": "Srikant"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "4636",
"DOI": "10.1038/s41467-026-71098-4",
"PMID": "41917013",
"PMCID": "PMC13199526",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-71098-4",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
31
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41467-024-47028-7 [code]
Native-state proteomics of Parvalbumin interneurons identifies unique molecular signatures and vulnerabilities to early Alzheimer's pathology.
Journal: Nature communications
In common: WGCNA, ggplot2, genetics / omics, mouse, cellular / molecular, 13 references, 6 authors
[2] doi:10.1016/j.mcpro.2025.101458 [code]
Proximity Labeling of the Tau Repeat Domain Enriches RNA-Binding Proteins That Are Altered in Alzheimer's Disease and Related Tauopathies.
Journal: Molecular & cellular proteomics : MCP
In common: WGCNA, ggplot2, genetics / omics, cellular / molecular, 10 references, 3 authors
[3] doi:10.1016/j.mcpro.2023.100546 [code]
Cellular Proteomic Profiling Using Proximity Labeling by TurboID-NES in Microglial and Neuronal Cell Lines.
Journal: Molecular & cellular proteomics : MCP
In common: tidyverse, genetics / omics, mouse, cellular / molecular, 13 references, 2 authors
[4] doi:10.1038/s41467-024-52297-3 [code]
Entorhinal cortex vulnerability to human APP expression promotes hyperexcitability and tau pathology.
Journal: Nature communications
In common: mouse, 4 references, 4 authors
[5] doi:10.1016/j.mcpro.2023.100678 [code]
Identification of State-Specific Proteomic and Transcriptomic Signatures of Microglia-Derived Extracellular Vesicles.
Journal: Molecular & cellular proteomics : MCP
In common: WGCNA, DESeq2, pheatmap, 2 other tools, genetics / omics, mouse, cellular / molecular, 3 references, 2 authors
[6] doi:10.1016/j.celrep.2022.111474 [code]
Activity-dependent translation dynamically alters the proteome of the perisynaptic astrocyte process.
Journal: Cell reports
In common: limma, data.table, ggplot2, 1 other tool, genetics / omics, cellular / molecular, 12 references
[7] doi:10.1016/j.mcpro.2022.100422
Deciphering Spatial Protein-Protein Interactions in Brain Using Proximity Labeling.
Journal: Molecular & cellular proteomics : MCP
In common: genetics / omics, cellular / molecular, 12 references
[8] doi:10.1038/s41467-025-63119-5 [code]
An integrated transcriptomic and proteomic map of the mouse hippocampus at synaptic resolution.
Journal: Nature communications
In common: DESeq2, pheatmap, data.table, 2 other tools, genetics / omics, mouse, cellular / molecular, 11 references
[9] doi:10.1038/s41467-025-64328-8 [code]
Integrative brain omics approach highlights sn-1 lysophosphatidylethanolamine in Alzheimer's dementia.
Journal: Nature communications
In common: WGCNA, limma, pheatmap, 3 other tools, genetics / omics, 2 references, 2 authors
[10] doi:10.1186/s13024-023-00638-z [code]
Quantitative proteomics of cerebrospinal fluid from African Americans and Caucasians reveals shared and divergent changes in Alzheimer's disease.
Journal: Molecular neurodegeneration
In common: WGCNA, ggplot2, tidyverse, genetics / omics, 5 references, 2 authors

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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1038/s41467-026-71098-4.svg)](https://openscicode.org/badge/doi/10.1038/s41467-026-71098-4)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1038/s41467-026-71098-4.svg
   :target: https://openscicode.org/badge/doi/10.1038/s41467-026-71098-4
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1038/s41467-026-71098-4"><img src="https://openscicode.org/badge/doi/10.1038/s41467-026-71098-4.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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