OSCR

Sex specific effects of adoptive Tregs transfer on the brain and periphery in maternal immune activation offspring rescuing immune dysregulation.

Code ↔ Paper

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

The 8 matches
  1. [1] § Materials and methods › RNA-sequencing preparation and analysis ↔ Master_QuantSeqAnalysis.sh, lines 166–206 · score 0.98 · alignIntronMax, alignMatesGapMax, alignSJDBoverhangMin, outFilterMismatchNmax, STAR, mm10
  2. [2] § Materials and methods › RNA-sequencing preparation and analysis ↔ LimmaVoom_DE_RNAseq_FC_sexlitter_removeP30.R, lines 300–346 · score 0.88 · calcNormFactors, contrasts.fit, library composition, eBayes, design matrix, Voom
  3. [3] § Materials and methods › RNA-sequencing preparation and analysis ↔ LimmaVoom_DE_RNAseq_CB_sexlitter.R, lines 341–387 · score 0.88 · calcNormFactors, contrasts.fit, library composition, eBayes, design matrix, Voom
  4. [4] § Materials and methods › Weighted gene co-expression network analysis ↔ WGCNA_allconditions.R, lines 423–478 · score 0.79 · soft thresholding power, co expression, TOM, adjacency, topological, fit
  5. [5] § Materials and methods › Weighted gene co-expression network analysis ↔ WGCNA_allconditions.R, lines 1436–1475 · score 0.70 · Gene module membership, Hub gene, log2RPKM, WGCNA
  6. [6] § Materials and methods › RNA-sequencing preparation and analysis ↔ DEGcomparisons.R, lines 476–522 · score 0.55 · GO term enrichments, compareCluster, cutoff, ClusterProfiler, genes
  7. [7] § Materials and methods › RNA-sequencing preparation and analysis ↔ LimmaVoom_DE_RNAseq_CB_sexlitter.R, lines 1039–1085 · score 0.55 · GO term enrichments, compareCluster, cutoff, ClusterProfiler, genes
  8. [8] § Materials and methods › Weighted gene co-expression network analysis ↔ WGCNA_allconditions.R, lines 326–328 · score 0.51 · median absolute deviation, WGCNA

Paper

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

The paper is loaded when this pane is shown.

The authors' code

R · 2,181 lines · 83 KB · no license · 3 matches

  1. #####################################################################################
  2. # Gene expression from human brain across development: brainspan
  3. ####################################################################################
  4. #BiocManager::install("edgeR")
  5. library(edgeR)
  7. #source("https://bioconductor.org/biocLite.R")
  8. #biocLite("limma")
  9. library(limma)
  10. library(dplyr)
  11. library(tidyr)
  12. library(ggplot2)
  13. library(cowplot)
  14. library(gplots)
  17. #source("http://bioconductor.org/biocLite.R")
  18. #biocLite("biomaRt")
  19. library(biomaRt)
  20. library(WGCNA);
  21. # The following setting is important, do not omit.
  22. options(stringsAsFactors = FALSE)
  24. ######################################################################
  25. # Basic function to convert mouse to human gene names
  26. human = useMart(biomart="ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl")
  27. mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
  29. # convertMouseGeneList <- function(x){
  30. # #for hg38:
  31. # genesV2 = getLDS(attributes = c("mgi_symbol","ensembl_gene_id"), filters = "mgi_symbol", values = x , mart = mouse, attributesL = c("hgnc_symbol","ensembl_gene_id","entrezgene"), martL = human, uniqueRows=T)
  32. # humanx <- unique(genesV2[,])
  33. #
  34. # # Print the first 6 genes found to the screen
  35. # print(head(humanx))
  36. # return(humanx)
  37. # }
  39. ######################################################################
  40. #Change to your working directory
  41. setwd("/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/")
  42. ######################################################################
  43. ############################################################################################################################
  44. ############################################################################################################################
  45. ############################################################################################################################
  46. #Initial data processing, sample cleaning and differential expression analysis
  47. ############################################################################################################################
  48. ############################################################################################################################
  49. ############################################################################################################################
  51. ############################################################################################################################
  52. #make RPKM data
  53. ############################################################################################################################
  55. #loop for combining GeneID and counts for each txt file
  56. path ="/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/counts"
  58. setwd(path)
  60. # out.file<-""
  61. # group<-""
  62. # file.names <- dir(path, pattern ='*collapsedUMI.counts.txt')
  64. # for(i in 1:length(file.names)){
  65. # file <- read.table(file.names[i],skip=1, header=TRUE, sep="\t", stringsAsFactors=FALSE)
  66. # counts <- data.frame(file$Geneid,file[,7])
  67. # #extract file name without count.txt
  68. # m <- as.data.frame(file.names[i])
  69. # names(m) <- c("split")
  70. # m$split <- as.character(m$split)
  71. # m <-tidyr::separate(m,split,into= c("name","stuff","stuff2","stuff3"),sep="\\.")
  72. #
  73. # print(m$name)
  74. # #replace column headings with names
  75. # names(counts) <- c("Geneid",m$name)
  76. #
  77. # #write to output file
  78. # out.file <- counts
  79. # write.table(out.file, file =paste(m$name, ".out.txt", sep=""), sep="\t",quote=F)
  80. # }
  83. #take files from loop above as new input
  84. input.files <-dir(path, pattern ='*out.txt')
  85. input.files <- as.data.frame(input.files)
  86. names(input.files) <- c("filenames")
  88. #get treatments names
  89. input.files$Sample.ID <- gsub(".out.txt","",input.files$filenames)
  91. #read in info
  92. Info3 <- read.csv("/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/MasterExperimentInfo.csv")
  94. #remove PA92FC30 based on MSDS
  95. Info3 <- Info3 %>% filter(Sample.ID !="PA92FC30")
  98. setwd("/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/counts")
  99. path <- c("/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/counts")
  100. RG <- readDGE(Info3$filenames, group=Info3$Group)
  101. colnames(RG$counts) <- gsub(".out","",colnames(RG$counts))
  103. #save
  104. path ="/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/WGCNA"
  105. setwd(path)
  106. save(RG,Info3, file="RNAseqDEGlist.RData")
  108. # number of unique transcripts
  109. dim(RG)
  110. #53801 177
  112. ###########################################################
  113. # Filter for one count per million in at least 1/4 of libraries = 130/4 = 32
  114. allsamples <- table(Info3$Dam.Treatment.Group,Info3$Pup.Treatment.Group,Info3$Sex.x, Info3$Brain.Region)
  115. allsamples
  117. #smallest group is n=6
  118. keep <- rowSums(cpm(RG)>1)>=6
  119. RGcpm <- RG[keep,]
  120. dim(RGcpm)
  121. # 21283 177
  123. geneid <- rownames(RGcpm) #ensemble IDs for biomart
  125. ##################################################################################
  126. #read in gene info
  127. ##################################################################################
  129. #source("http://bioconductor.org/biocLite.R")
  130. #biocLite("biomaRt")
  131. library(biomaRt)
  133. #mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
  134. #us_mart <- useEnsembl(biomart = "ensembl", mirror = "uswest")
  135. us_mart <- useEnsembl(biomart = "ensembl", mirror = "uswest",dataset = "mmusculus_gene_ensembl")
  137. #head(listAttributes(human),100)
  138. #attributes <- listAttributes(us_mart)
  139. #attributes[grep("exon", attributes$name),]
  141. genes <- getBM(attributes = c("ensembl_gene_id","chromosome_name","start_position","end_position","mgi_symbol","external_gene_name","description"),
  142. filter= "ensembl_gene_id",
  143. values = rownames(RGcpm),
  144. mart = us_mart)
  148. genes$genelength <- abs(genes$end_position - genes$start_position)
  150. #remove duplicates:keeps only first entry
  151. genes <- genes[!duplicated(genes$ensembl_gene_id),]
  153. #match order: df[match(target, df$name),]
  154. genes <- genes[match(rownames(RGcpm),genes$ensembl_gene_id),]
  156. ##########################################################################
  157. #add gene info to DEGlist object
  158. ##########################################################################
  160. #add into DEGlist
  161. RGcpm$genes <- genes
  163. #save
  164. save(RG,RGcpm,genes, file="DEGlist.RData")
  166. ###########################################################
  167. #filtering plot
  168. ###########################################################
  169. library(RColorBrewer)
  170. nsamples <- ncol(RGcpm)
  172. colourCount = nsamples
  173. getPalette = colorRampPalette(brewer.pal(9, "Set1"))
  174. fill=getPalette(colourCount)
  176. #plot:
  177. pdf('FilteringCPM_plots.pdf')
  178. par(mfrow=c(1,2))
  180. #prefilter:
  181. lcpm <- cpm(RG, log=TRUE, prior.count=2)
  184. plot(density(lcpm[,1]), col=fill[1], lwd=2, ylim=c(0,0.5), las=2,
  185. main="", xlab="")
  186. title(main="A. Raw data", xlab="Log-cpm")
  187. abline(v=0, lty=3)
  188. for (i in 2:nsamples){
  189. den <- density(lcpm[,i])
  190. lines(den$x, den$y, col=fill[i], lwd=2)
  191. }
  192. #legend("topright", Samples, text.col=fill, bty="n")
  194. #filtered data
  195. #og-CPM of zero threshold (equivalent to a CPM value of 1) used in the filtering ste
  196. lcpm <- cpm(RGcpm, log=TRUE, prior.count=2)
  197. plot(density(lcpm[,1]), col=fill[1], lwd=2, ylim=c(0,0.5), las=2,
  198. main="", xlab="")
  199. title(main="B. Filtered data", xlab="Log-cpm")
  200. abline(v=0, lty=3)
  201. for (i in 2:nsamples){
  202. den <- density(lcpm[,i])
  203. lines(den$x, den$y, col=fill[i], lwd=2)
  204. }
  205. #legend("topright", Samples, text.col=fill, bty="n")
  206. dev.off()
  209. ###########################################################
  210. ###########################################################
  211. #reset library sizes
  212. RGcpm$samples$lib.size <- colSums(RGcpm$counts)
  215. #plot library sizes
  216. pdf('LibrarySizes.pdf',w=30,h=8)
  217. barplot(RGcpm$samples$lib.size,names=colnames(RGcpm),las=2)
  218. # Add a title to the plot
  219. title("Barplot of library sizes")
  220. dev.off()
  223. # Get log2 counts per million
  224. logcounts <- cpm(RGcpm,log=TRUE)
  225. # Check distributions of samples using boxplots
  226. pdf('NonNormalizedLogCPM.pdf',w=30,h=10)
  227. boxplot(logcounts, xlab="", ylab="Log2 counts per million",las=2)
  228. # Let's add a blue horizontal line that corresponds to the median logCPM
  229. abline(h=median(logcounts),col="blue")
  230. title("Boxplots of logCPMs (unnormalised)")
  231. dev.off()
  233. ##########################################################################
  234. #Diagnosis * treatment + sex with subjects correlation removed
  235. ##################################################################################
  236. #make design matrix
  237. Info3$Sex <- factor(Info3$Sex.x, levels=c("M","F"))
  239. Info3$Dam.Treatment.Group <- factor(Info3$Dam.Treatment.Group,levels = c("Saline","PolyIC"))
  241. Info3$Pup.Treatment.Group <- factor(Info3$Pup.Treatment.Group,levels = c("Saline","Treg"))
  243. Info3$Brain.Region <- factor(Info3$Brain.Region)
  245. Info3$Group <- factor(Info3$Group)
  247. #set design matrix
  248. design <- model.matrix(~0+Group, data=Info3) #if using a 0 intercept must set up contrasts
  250. colnames(design)
  253. ##################################################################################
  254. #Normalization TMM and Voom
  255. ##################################################################################
  256. #The two steps refer to different aspects of normalization. CPM "normalization" accounts for library size differences between samples, and produces normalized values that can be compared on an absolute scale (e.g., for filtering). TMM normalization accounts for composition bias, and computes normalization factors for comparing between libraries on a relative scale. CPM normalization doesn't account for composition bias, and TMM normalization doesn't produce normalized values. Thus, you need both steps in the analysis pipeline. This isn't a problem, as the two steps aren't really redundant.
  257. #TMM normalization for library composition
  258. DGE=calcNormFactors(RGcpm,method =c("TMM"))
  260. pdf('VoomTrend.pdf',w=6,h=4)
  261. v=voom(DGE,design,plot=T)
  262. dev.off()
  264. corfit <- duplicateCorrelation(v, design, block=Info3$Mouse.Number)
  266. fit <- lmFit(v, design, block = Info3$Mouse.Number, correlation = corfit$consensus)
  269. #compute RPKM off adjusted library sizes
  270. RPKM <- rpkm(DGE,gene.length =DGE$genes$genelength,normalized.lib.sizes = TRUE, log=F)
  272. #save
  273. save(RG,RGcpm,Info3,design,genes, RPKM,v,file="DEGlist.RData")
  276. ############################################################################################################################
  277. ############################################################################################################################
  278. ############################################################################################################################
  279. #WGCNA
  280. ############################################################################################################################
  281. ############################################################################################################################
  282. ############################################################################################################################
  284. #make WGCNA table
  285. #remove extra columns:
  286. DF <- RPKM
  288. #make gene matrix with genes as columns and rows as samples
  289. DF2 <- t(DF)
  290. DF2 <- as.data.frame(DF2)
  291. write.csv(DF2,"WGCNA_RPKMinputmatrix.csv")
  293. ################## process data ########################
  294. #detect genes with missing values or samples with missing values
  295. gsg = goodSamplesGenes(DF2, verbose = 3);
  296. gsg$allOK # TRUE
  298. #remove the genes with too many missing values
  299. if (!gsg$allOK)
  300. {
  301. # Optionally, print the gene and sample names that were removed:
  302. if (sum(!gsg$goodGenes)>0)
  303. printFlush(paste("Removing genes:", paste(names(DF2)[!gsg$goodGenes], collapse = ", ")));
  304. if (sum(!gsg$goodSamples)>0)
  305. printFlush(paste("Removing samples:", paste(rownames(DF2)[!gsg$goodSamples], collapse = ", ")));
  306. # Remove the offending genes and samples from the data:
  307. DF2 = DF2[gsg$goodSamples, gsg$goodGenes]
  308. }
  310. gsg = goodSamplesGenes(DF2, verbose = 3);
  311. gsg$allOK # TRUE
  314. #gene with RPKM value of 2 or higher in at least one sample.
  315. #max(col)>=2
  316. #https://stackoverflow.com/questions/24212739/how-to-find-the-highest-value-of-a-column-in-a-data-frame-in-r/24212879
  317. colMax <- function(data) sapply(data, max, na.rm = TRUE)
  319. DF3 <- DF2[, colMax(DF2) >= 2]
  320. dim(DF2)
  321. dim(DF3)
  323. #log2 transform:
  324. Log2DF3 <- log2(DF3+1)
  326. # median absolute deviation
  327. #remove if = 0 https://support.bioconductor.org/p/65124/
  328. colMad <- function(data) sapply(data, mad, na.rm = TRUE)
  330. DF4 = Log2DF3[, colMad(Log2DF3) != 0]
  331. dim(DF4) #6464 genes left
  333. write.csv(DF4,"WGCNA_log2RPKMinputmatrix.csv")
  337. #=====================================================================================
  338. #Next we cluster the samples (in contrast to clustering genes that will come later) to see if there are any obvious outliers.
  339. sampleTree = hclust(dist(DF4), method = "average");
  340. # Plot the sample tree: Open a graphic output window of size 12 by 9 inches
  341. # The user should change the dimensions if the window is too large or too small.
  342. sizeGrWindow(12,9)
  345. pdf(file = "sampleClustering_preCut.pdf", width = 12, height = 9);
  346. par(cex = 0.6);
  347. par(mar = c(0,4,2,0))
  348. plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
  349. cex.axis = 1.5, cex.main = 2)
  351. # Plot a line to show the cut
  352. #remove two samples E18 tube 8 and P60 tube 50
  353. abline(h = 160, col = "red");
  355. dev.off()
  357. # Determine cluster under the line
  358. clust = cutreeStatic(sampleTree, cutHeight = 160, minSize = 10)
  359. table(clust)
  360. # clust 1 contains the samples we want to keep. > all in this case
  361. keepSamples = (clust==1)
  362. datExpr = DF4[keepSamples, ]
  363. nGenes = ncol(datExpr)
  364. nSamples = nrow(datExpr)
  366. #=====================================================================================
  367. #phenotype data:
  368. #Read in experient information sheet
  369. names(Info3)
  370. info <- Info3
  372. #filter to include only seq files:
  373. info$Sample.ID <- as.character(info$Sample.ID)
  375. infowt <- info[which(info$Sample.ID %in% rownames(datExpr)),]
  377. table(infowt$Sample.ID %in% rownames(datExpr)) #177 samples
  380. #match sample names and trait data:
  381. Samples = rownames(datExpr)
  382. traitRows = match(Samples, infowt$Sample.ID)
  383. datTraits = infowt[traitRows, ]
  385. datTraits<- as.data.frame(datTraits)
  386. datTraits$Sample.ID
  388. #change to numeric: sex #F =2 = red, M= 1 = white
  389. datTraits$sex2 <- as.numeric(as.factor(datTraits$Sex))
  391. #change to numeric: dam treatment 1=saline, 2=polyIC
  392. datTraits$Dam.Treatment.Group2 <- as.numeric(datTraits$Dam.Treatment.Group)
  394. #change to numeric: pup treatment 1=saline, 2=treg
  395. datTraits$Pup.Treatment.Group2 <- as.numeric(datTraits$Pup.Treatment.Group)
  397. #change to numeric: 1=Cerebellum, 2= Frontal, 3= Cortex Hippocampus
  398. datTraits$Brain.Region2 <- as.numeric(datTraits$Brain.Region)
  400. #numeric traits only
  401. numerictraits <- datTraits[,c(20,21,19,22)]
  402. rownames(numerictraits) <- datTraits$Sample.ID
  404. #We now have the expression data in the variable datExpr, and the corresponding clinical traits in the variable datTraits. Before we continue with network construction and module detection, we visualize how the clinical traits relate to the sample dendrogram.
  405. # Re-cluster samples
  406. sampleTree2 = hclust(dist(datExpr), method = "average")
  407. # Convert traits to a color representation: white means low, red means high, grey means missing entry
  408. traitColors1 = numbers2colors(numerictraits$Dam.Treatment.Group2,signed = FALSE);
  409. traitColors2 = numbers2colors(numerictraits$Pup.Treatment.Group2,signed = FALSE);
  410. traitColors3 = numbers2colors(numerictraits$sex2,signed = FALSE);
  411. traitColors <-cbind(traitColors1,traitColors2,traitColors3)
  412. colnames(traitColors) <- c("Dam treatment","Pup treatment","sex")
  413. # Plot the sample dendrogram and the colors underneath.
  415. pdf(file = "sampleClustering_withtraits.pdf", width = 12, height = 9)
  416. plotDendroAndColors(sampleTree2, traitColors, groupLabels = colnames(traitColors),
  417. main = "Sample dendrogram and trait heatmap")
  418. dev.off()
  420. save(numerictraits,datTraits,datExpr, file = "WTWGCNAnetworkConstruction-inputdata.RData")
  423. #=====================================================================================
  424. #choice of the soft thresholding power β to which co-expression similarity is raised to calculate adjacency
  425. # Choose a set of soft-thresholding powers
  426. powers = c(c(1:30))
  427. # Call the network topology analysis function
  428. sft = pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed",corFnc="bicor" ,verbose = 5)
  429. # Plot the results:
  430. sizeGrWindow(9, 5)
  431. par(mfrow = c(1,2));
  432. cex1 = 0.9;
  433. # Scale-free topology fit index as a function of the soft-thresholding power
  434. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  435. xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
  436. main = paste("Scale independence"));
  437. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  438. labels=powers,cex=cex1,col="blue");
  439. # this line corresponds to using an R^2 cut-off of h
  440. abline(h=0.80,col="black")
  441. # Mean connectivity as a function of the soft-thresholding power
  442. plot(sft$fitIndices[,1], sft$fitIndices[,5],
  443. xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
  444. main = paste("Mean connectivity"))
  445. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="blue")
  447. write.csv(sft, "SoftThresholding_powercalc.csv")
  448. #lowest power for which the scale-free topology fit index of > 0.8 ->5
  450. #=====================================================================================
  451. #Run this section on an external server, save the network data and load into R studio
  452. #=====================================================================================
  453. #Co-expression similarity and adjacency
  454. library(WGCNA)
  455. allowWGCNAThreads()
  457. # The following setting is important, do not omit.
  458. options(stringsAsFactors = FALSE)
  460. #load data
  464. Allnet = blockwiseModules(datExpr, maxBlockSize = 50000,
  465. power = 5, TOMType = "signed", minModuleSize = 20,
  466. networkType = "signed",
  467. corType = "bicor", #biweight midcorrelation
  468. maxPOutliers = 0.05, #forces bicor to never regard more than the specified proportion of samples as outliers.
  469. reassignThreshold = 0,
  470. numericLabels = TRUE,
  471. saveTOMs = FALSE,
  472. nThreads = 12,
  473. #saveTOMFileBase = "TOM-BilboMGExpression",
  474. verbose = 3)
  476. save(Allnet, file = "NetworkConstruction-auto_WTWGCNA.RData")
  477. #=====================================================================================
  478. table(Allnet$colors)
  480. # Convert labels to colors for plotting
  481. AllnetModuleColors = labels2colors(Allnet$colors)
  484. AllnetTree = Allnet$dendrograms[[1]]
  486. #plot the gene dendrogram and the corresponding module colors
  487. sizeGrWindow(8,6);
  489. pdf(file = "AllnetDendrogram_allsamples.pdf", wi = 8, he = 6)
  490. plotDendroAndColors(AllnetTree, AllnetModuleColors,
  491. "Module colors",
  492. dendroLabels = FALSE, hang = 0.03,
  493. addGuide = TRUE, guideHang = 0.05,
  494. main = "All Samples Cluster Dendrogram")
  495. dev.off()
  498. #=====================================================================================
  500. # Calculate eigengenes
  501. MEList = moduleEigengenes(datExpr, colors = AllnetModuleColors)
  502. MEs = MEList$eigengenes
  503. # Calculate dissimilarity of module eigengenes
  504. MEDiss = 1-cor(MEs);
  505. # Cluster module eigengenes
  506. METree = hclust(as.dist(MEDiss), method = "average");
  507. # Plot the result
  508. sizeGrWindow(7, 6)
  510. pdf(file = "Pretrim_allnet_dendrogram.pdf", wi = 9, he = 6)
  511. plot(METree, main = "Clustering of module eigengenes",
  512. xlab = "", sub = "")
  513. # Plot the cut line into the dendrogram
  514. abline(h=0.6, col = "red")
  516. dev.off()
  520. #=====================================================================================
  521. #merge down to fewer modules based on plot height
  522. MEDissThres = 0.1
  524. # Call an automatic merging function
  525. merge = mergeCloseModules(datExpr, AllnetModuleColors, cutHeight = MEDissThres, verbose = 3)
  526. # The merged module colors
  527. mColors = merge$colors;
  528. # Eigengenes of the new merged modules:
  529. mergedMEs = merge$newMEs;
  531. #replot
  532. sizeGrWindow(12, 9)
  533. pdf(file = "Allnet_postmerg.1MEClustering.pdf", wi = 9, he = 6)
  535. plotDendroAndColors(AllnetTree, mColors,
  536. "Module colors",
  537. dendroLabels = FALSE, hang = 0.03,
  538. addGuide = TRUE, guideHang = 0.05,
  539. main = "All Samples Cluster Dendrogram")
  540. dev.off()
  542. #=====================================================================================
  543. # Rename to moduleColors
  544. moduleColors = mColors
  545. # Construct numerical labels corresponding to the colors
  546. colorOrder = c("grey", standardColors(50));
  547. moduleLabels = match(moduleColors, colorOrder)-1;
  548. MEs = mergedMEs;
  549. # Save module colors and labels for use in subsequent parts
  550. save(MEs, moduleLabels, moduleColors, AllnetTree, datTraits, numerictraits,datExpr,Allnet,file = "allsamplesWTMGtimecourse-networkpostmerge.1.RData")
  552. #load("allsamplesbilboMGtimecourse-networkpostmerge.25.RData")
  553. #=====================================================================================
  554. #repeat clustering of new MEs
  555. # Calculate dissimilarity of module eigengenes
  556. MEDiss = 1-cor(MEs);
  557. # Cluster module eigengenes
  558. METree = hclust(as.dist(MEDiss), method = "average");
  560. # Plot the result
  561. sizeGrWindow(7, 6)
  562. pdf(file = "Pretrim_allnet_dendrogram_postmerge.pdf", wi = 9, he = 6)
  563. plot(METree, main = "Clustering of module eigengenes after merging modules",
  564. xlab = "", sub = "")
  566. dev.off()
  568. #=====================================================================================
  569. #=====================================================================================
  570. #correlate eigengenes with external traits and look for the most signi cant associations
  571. # Define numbers of genes and samples
  572. nGenes = ncol(datExpr);
  573. nSamples = nrow(datExpr);
  574. # Recalculate MEs with color labels
  575. MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
  576. MEs = orderMEs(MEs0)
  577. moduleTraitCor = cor(MEs, numerictraits, use = "p");
  578. moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
  580. #FDR correct the pvalues
  581. #sapply(pval,p.adjust,method="fdr") #per column
  583. #FDR correction on entire matrix
  584. FDR <- matrix(p.adjust(as.vector(moduleTraitPvalue), method='fdr'),ncol=ncol(moduleTraitPvalue))
  585. colnames(FDR) <- colnames(moduleTraitPvalue)
  586. rownames(FDR) <- rownames(moduleTraitPvalue)
  589. corout <- cbind(moduleTraitCor,FDR)
  590. write.csv(corout, "PearsonsCorrelations_mod-trait.csv",row.names = T)
  592. #We color code each association by the correlation value:
  594. sizeGrWindow(10,6)
  596. pdf(file = "module-traitrelationshipsFDRcorrected.pdf", wi = 8.5, he = 11)
  598. # Will display correlations and their p-values
  599. textMatrix = paste(signif(moduleTraitCor, 2), "\ (", #pearson correlation coefficient, space, (FDR corrected pvalue)
  600. signif(FDR, 4), ")", sep = "");
  601. dim(textMatrix) = dim(moduleTraitCor)
  602. par(mar = c(6, 8.5, 3, 3));
  603. # Display the correlation values within a heatmap plot
  604. labeledHeatmap(Matrix = moduleTraitCor,
  605. xLabels = names(numerictraits),
  606. yLabels = names(MEs),
  607. ySymbols = names(MEs),
  608. colorLabels = FALSE,
  609. colors = blueWhiteRed(50),
  610. textMatrix = textMatrix,
  611. setStdMargins = FALSE,
  612. cex.text = 0.9,
  613. zlim = c(-1,1),
  614. main = paste("Module-trait relationships"))
  616. dev.off()
  618. #=====================================================================================
  619. #ANOVA for trait associations
  620. #=====================================================================================
  622. library(nlme)
  623. library(lsmeans)
  624. datTraits <- as.data.frame(datTraits)
  625. datTraits$sex <- factor(datTraits$sex)
  626. datTraits$Dam.Treatment.Group <- factor(datTraits$Dam.Treatment.Group)
  627. datTraits$Pup.Treatment.Group <- factor(datTraits$Pup.Treatment.Group)
  628. datTraits$Brain.Region <- factor(datTraits$Brain.Region)
  629. datTraits$subject <- datTraits$Sample.ID
  632. MEs$Sample.ID <- rownames(MEs)
  633. datComb <- merge(datTraits,MEs, by = "Sample.ID")
  634. datComb2 <- datComb %>% gather(module, ME, 25:ncol(datComb))
  636. write.csv(datComb2,"Moduel-ME_dataforANOVA.csv")
  638. mod <- unique(datComb2$module)
  639. #remove grey
  640. mod <- head(mod, -1)
  642. anova_out <- NULL
  643. for (i in mod) {
  644. print(i)
  645. #select data
  646. dattmp <- datComb2 %>% dplyr::filter(module == i)
  647. #model
  648. tmp <- lme(ME ~ Dam.Treatment.Group*Pup.Treatment.Group*sex*Brain.Region, ~1|Sample.ID, data = dattmp)
  649. #model indcludes random effects subject ID
  650. anova <- anova(tmp, type = "marginal") #marginal gives Type 3 SS for ANOVA
  651. #write anova to output file
  652. anova <- as.data.frame(anova)
  653. anova$condition <- rownames(anova)
  654. anova <- anova[-1,]
  655. anova2 <- anova %>% gather(stats,values,1:4) %>%
  656. mutate(col = paste(condition,stats, sep="_")) %>%
  657. dplyr::select(col,values) %>%
  658. spread(col,values)
  659. anova2$module <- paste(i)
  660. #save to loop
  661. anova_out <- rbind(anova_out,anova2)
  662. }
  665. write.csv(anova_out,"MixedModelAnova_moudle-traits.csv")
  669. #include post hocs
  670. anova_out <- NULL
  671. ME_posthoc <- NULL
  672. for (i in mod) {
  673. print(i)
  674. #select data
  675. dattmp <- datComb2 %>% dplyr::filter(module == i)
  676. #model
  677. tmp <- lme(ME ~ Dam.Treatment.Group*Pup.Treatment.Group*sex*Brain.Region, ~1|Sample.ID, data = dattmp)
  678. #model indcludes random effects subject ID
  679. anova <- anova(tmp, type = "marginal") #marginal gives Type 3 SS for ANOVA
  680. #write anova to output file
  681. anova <- as.data.frame(anova)
  682. anova$condition <- rownames(anova)
  683. anova <- anova[-1,]
  684. anova2 <- anova %>% gather(stats,values,1:4) %>%
  685. mutate(col = paste(condition,stats, sep="_")) %>%
  686. dplyr::select(col,values) %>%
  687. spread(col,values)
  688. anova2$module <- paste(i)
  689. #save to loop
  690. anova_out <- rbind(anova_out,anova2)
  692. ## POSTHOC starts here ##
  693. #create a reference grid model object
  694. ME_refgrid <- ref.grid(tmp)
  696. #creates a fitted model using the reference grid, based on treatment and brain region interactions
  697. ME_lsmeans <- lsmeans(ME_refgrid, ~Dam.Treatment.Group*Pup.Treatment.Group|sex*Brain.Region)
  699. #summarize paired comparisons
  700. ME_lsmeans_summary <- summary(pairs(ME_lsmeans, adjust = "none"))
  701. ME_lsmeans_summary <- as.data.frame(ME_lsmeans_summary)
  703. #posthocHSD = contrast(refgrid, method = "pairwise",adjust = "none")
  704. #outsum <- as.data.frame(summary(posthocHSD))
  707. #get only control vs LPS
  708. # outsum <- outsum[grepl("CONTROL", outsum$contrast),]
  710. #correct for multiple comparisons
  711. ME_lsmeans_summary$padjust <- p.adjust(ME_lsmeans_summary$p.value, method = "BH")
  713. #paste module information
  714. ME_lsmeans_summary$module <- paste(i)
  716. #save to output
  717. ME_posthoc <- rbind(ME_posthoc,ME_lsmeans_summary)
  719. }
  721. write.csv(ME_posthoc,"BHposthocs_moudle-traits.csv")
  724. #graph MEs
  725. pdf(file = "boxplot_MEs.pdf", wi = 22, he = 30)
  726. #
  727. p <- ggplot(datComb2, aes(x=Dam.Treatment.Group, y=ME), group = Pup.Treatment.Group) +
  728. facet_wrap(~module*Brain.Region*Sex.x,ncol = 6, scales = "free") +
  729. stat_summary(geom = "boxplot",
  730. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  731. position = "dodge", aes(fill=Pup.Treatment.Group))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  732. geom_point(position = position_dodge(width = 0.90),aes(group=Pup.Treatment.Group, colour=Sex.x)) +
  733. scale_fill_manual(values = cbPalette) +
  734. theme_cowplot(font_size = 15)+
  735. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  736. strip.background = element_rect(fill="#B0C4DE"),
  737. # legend.position="none",
  738. strip.text.x = element_text(size = 20))+
  739. xlab(label = c("Dam Treatment")) +
  740. ylab(label = c("ME"))
  741. p
  742. dev.off()
  743. #
  746. #graph MEs Cerebellum ME blue
  747. pdf(file = "boxplot_MEs_CB_blue.pdf", wi = 6, he = 4)
  748. #
  749. CBblue <- datComb2 %>% filter(module == "MEblue") %>% filter(Brain.Region == "Cerebellum")
  751. p <- ggplot(CBblue, aes(x=Dam.Treatment.Group, y=ME), group = Pup.Treatment.Group) +
  752. facet_wrap(~Sex.x,ncol = 2) +
  753. stat_summary(geom = "boxplot",
  754. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  755. position = "dodge", aes(fill=Pup.Treatment.Group))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  756. geom_point(position = position_dodge(width = 0.90),aes(group=Pup.Treatment.Group)) +
  757. scale_fill_manual(values = cbPalette) +
  758. theme_cowplot(font_size = 15)+
  759. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  760. strip.background = element_rect(fill="#B0C4DE"),
  761. # legend.position="none",
  762. strip.text.x = element_text(size = 20))+
  763. xlab(label = c("Dam Treatment")) +
  764. ylab(label = c("ME"))
  765. p
  766. dev.off()
  767. #
  768. #graph MEs Cerebellum ME turquoise
  769. pdf(file = "boxplot_MEs_CB_turquoise.pdf", wi = 6, he = 4)
  770. #
  771. CBturquoise <- datComb2 %>% filter(module == "MEturquoise") %>% filter(Brain.Region == "Cerebellum")
  773. p <- ggplot(CBturquoise, aes(x=Dam.Treatment.Group, y=ME), group = Pup.Treatment.Group) +
  774. facet_wrap(~Sex.x,ncol = 2) +
  775. stat_summary(geom = "boxplot",
  776. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  777. position = "dodge", aes(fill=Pup.Treatment.Group))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  778. geom_point(position = position_dodge(width = 0.90),aes(group=Pup.Treatment.Group)) +
  779. scale_fill_manual(values = cbPalette) +
  780. theme_cowplot(font_size = 15)+
  781. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  782. strip.background = element_rect(fill="#B0C4DE"),
  783. # legend.position="none",
  784. strip.text.x = element_text(size = 20))+
  785. xlab(label = c("Dam Treatment")) +
  786. ylab(label = c("ME"))
  787. p
  788. dev.off()
  789. #
  791. #graph MEs HC ME turquoise
  792. pdf(file = "boxplot_MEs_HC_turquoise.pdf", wi = 6, he = 4)
  793. #
  794. HCturquoise <- datComb2 %>% filter(module == "MEturquoise") %>% filter(Brain.Region == "Hippocampus")
  796. p <- ggplot(HCturquoise, aes(x=Dam.Treatment.Group, y=ME), group = Pup.Treatment.Group) +
  797. facet_wrap(~Sex.x,ncol = 2) +
  798. stat_summary(geom = "boxplot",
  799. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  800. position = "dodge", aes(fill=Pup.Treatment.Group))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  801. geom_point(position = position_dodge(width = 0.90),aes(group=Pup.Treatment.Group)) +
  802. scale_fill_manual(values = cbPalette) +
  803. theme_cowplot(font_size = 15)+
  804. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  805. strip.background = element_rect(fill="#B0C4DE"),
  806. # legend.position="none",
  807. strip.text.x = element_text(size = 20))+
  808. xlab(label = c("Dam Treatment")) +
  809. ylab(label = c("ME"))
  810. p
  811. dev.off()
  812. #
  814. #graph MEs Cerebellum ME turquoise
  815. pdf(file = "boxplot_MEs_FC_green.pdf",wi = 6, he = 4)
  816. #
  817. FCgreen <- datComb2 %>% filter(module == "MEgreen") %>% filter(Brain.Region == "Frontal Cortex")
  819. p <- ggplot(FCgreen, aes(x=Dam.Treatment.Group, y=ME), group = Pup.Treatment.Group) +
  820. facet_wrap(~Sex.x,ncol = 2) +
  821. stat_summary(geom = "boxplot",
  822. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  823. position = "dodge", aes(fill=Pup.Treatment.Group))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  824. geom_point(position = position_dodge(width = 0.90),aes(group=Pup.Treatment.Group)) +
  825. scale_fill_manual(values = cbPalette) +
  826. theme_cowplot(font_size = 15)+
  827. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  828. strip.background = element_rect(fill="#B0C4DE"),
  829. # legend.position="none",
  830. strip.text.x = element_text(size = 20))+
  831. xlab(label = c("Dam Treatment")) +
  832. ylab(label = c("ME"))
  833. p
  834. dev.off()
  835. #
  837. pdf(file = "boxplot_MEs_FC_yellow.pdf", wi = 6, he = 4)
  838. #
  839. FCyellow <- datComb2 %>% filter(module == "MEyellow") %>% filter(Brain.Region == "Frontal Cortex")
  841. p <- ggplot(FCyellow, aes(x=Dam.Treatment.Group, y=ME), group = Pup.Treatment.Group) +
  842. facet_wrap(~Sex.x,ncol = 2) +
  843. stat_summary(geom = "boxplot",
  844. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  845. position = "dodge", aes(fill=Pup.Treatment.Group))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  846. geom_point(position = position_dodge(width = 0.90),aes(group=Pup.Treatment.Group)) +
  847. scale_fill_manual(values = cbPalette) +
  848. theme_cowplot(font_size = 15)+
  849. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  850. strip.background = element_rect(fill="#B0C4DE"),
  851. # legend.position="none",
  852. strip.text.x = element_text(size = 20))+
  853. xlab(label = c("Dam Treatment")) +
  854. ylab(label = c("ME"))
  855. p
  856. dev.off()
  857. #
  862. datComb2 <- read.csv("/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/WGCNA/Moduel-ME_dataforANOVA.csv")
  864. #graph MEs
  865. datComb2$Pup.Treatment.Group <- gsub("Treg","Tregs",datComb2$Pup.Treatment.Group)
  866. datComb2$Treatment <- paste(datComb2$Dam.Treatment.Group, datComb2$Pup.Treatment.Group, sep="-")
  867. datComb2$Treatment <- factor(datComb2$Treatment,levels=c("Saline-Saline","Saline-Tregs","PolyIC-Saline","PolyIC-Tregs"))
  869. #Saline–Saline (dark blue) — #777787
  871. #Saline–Treg (light blue) — #747788
  873. #Poly I:C–Saline (red) — #b15b53
  875. #Poly I:C–Treg (light red/pink) — #b2817f
  876. cbPalette <- c("#5756f9","#92bffa","#fe0100","#fe8080")
  878. #graph MEs Cerebellum ME blue
  879. pdf(file = "boxplot_MEs_CB_blue_wide.pdf", wi = 8, he = 6)
  880. #
  881. CBblue <- datComb2 %>% filter(module == "MEblue") %>% filter(Brain.Region == "Cerebellum")
  883. p <- ggplot(CBblue, aes(x = Treatment, y = ME, fill = Treatment)) +
  884. facet_wrap(~ Sex.x, ncol = 2) +
  885. stat_summary(
  886. geom = "boxplot",
  887. fun.data = function(x) {
  888. qs <- quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  889. setNames(qs, c("ymin", "lower", "middle", "upper", "ymax"))
  890. },
  891. position = position_dodge2(width = 0.9, preserve = "single")
  892. ) +
  893. geom_point(position = position_dodge(width = 0.9), aes(group = Treatment), alpha = 0.8) +
  894. scale_fill_manual(values = cbPalette) +
  895. # two-line x labels: keep "-" on the top line
  896. scale_x_discrete(labels = function(x) {
  897. parts <- strsplit(x, "\\s*[\\-–—]\\s*", perl = TRUE)
  898. sapply(seq_along(parts), function(i) {
  899. if (length(parts[[i]]) >= 2) paste0(parts[[i]][1], " -\n", parts[[i]][2]) else x[i]
  900. })
  901. }) +
  902. theme_cowplot(font_size = 15) +
  903. theme(
  904. axis.text.x = element_text(angle = 0, hjust = 0.5, lineheight = 0.9),
  905. strip.background = element_rect(fill = "#B0C4DE"),
  906. strip.text.x = element_text(size = 20),
  907. legend.position = "none"
  908. ) +
  909. xlab("Treatment") +
  910. ylab("ME")
  912. p
  914. dev.off()
  915. #
  916. #graph MEs Cerebellum ME turquoise
  917. pdf(file = "boxplot_MEs_CB_turquoise_wide.pdf", wi = 8, he = 6)
  918. #
  919. CBturquoise <- datComb2 %>% filter(module == "MEturquoise") %>% filter(Brain.Region == "Cerebellum")
  921. p <- ggplot(CBturquoise, aes(x = Treatment, y = ME, fill = Treatment)) +
  922. facet_wrap(~ Sex.x, ncol = 2) +
  923. stat_summary(
  924. geom = "boxplot",
  925. fun.data = function(x) {
  926. qs <- quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  927. setNames(qs, c("ymin", "lower", "middle", "upper", "ymax"))
  928. },
  929. position = position_dodge2(width = 0.9, preserve = "single")
  930. ) +
  931. geom_point(position = position_dodge(width = 0.9), aes(group = Treatment), alpha = 0.8) +
  932. scale_fill_manual(values = cbPalette) +
  933. # two-line x labels: keep "-" on the top line
  934. scale_x_discrete(labels = function(x) {
  935. parts <- strsplit(x, "\\s*[\\-–—]\\s*", perl = TRUE)
  936. sapply(seq_along(parts), function(i) {
  937. if (length(parts[[i]]) >= 2) paste0(parts[[i]][1], " -\n", parts[[i]][2]) else x[i]
  938. })
  939. }) +
  940. theme_cowplot(font_size = 15) +
  941. theme(
  942. axis.text.x = element_text(angle = 0, hjust = 0.5, lineheight = 0.9),
  943. strip.background = element_rect(fill = "#B0C4DE"),
  944. strip.text.x = element_text(size = 20),
  945. legend.position = "none"
  946. ) +
  947. xlab("Treatment") +
  948. ylab("ME")
  950. p
  952. dev.off()
  953. #
  955. #graph MEs HC ME turquoise
  956. pdf(file = "boxplot_MEs_HC_turquoise_wide.pdf", wi = 8, he = 6)
  957. #
  958. HCturquoise <- datComb2 %>% filter(module == "MEturquoise") %>% filter(Brain.Region == "Hippocampus")
  960. p <- ggplot(HCturquoise, aes(x = Treatment, y = ME, fill = Treatment)) +
  961. facet_wrap(~ Sex.x, ncol = 2) +
  962. stat_summary(
  963. geom = "boxplot",
  964. fun.data = function(x) {
  965. qs <- quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  966. setNames(qs, c("ymin", "lower", "middle", "upper", "ymax"))
  967. },
  968. position = position_dodge2(width = 0.9, preserve = "single")
  969. ) +
  970. geom_point(position = position_dodge(width = 0.9), aes(group = Treatment), alpha = 0.8) +
  971. scale_fill_manual(values = cbPalette) +
  972. # two-line x labels: keep "-" on the top line
  973. scale_x_discrete(labels = function(x) {
  974. parts <- strsplit(x, "\\s*[\\-–—]\\s*", perl = TRUE)
  975. sapply(seq_along(parts), function(i) {
  976. if (length(parts[[i]]) >= 2) paste0(parts[[i]][1], " -\n", parts[[i]][2]) else x[i]
  977. })
  978. }) +
  979. theme_cowplot(font_size = 15) +
  980. theme(
  981. axis.text.x = element_text(angle = 0, hjust = 0.5, lineheight = 0.9),
  982. strip.background = element_rect(fill = "#B0C4DE"),
  983. strip.text.x = element_text(size = 20),
  984. legend.position = "none"
  985. ) +
  986. xlab("Treatment") +
  987. ylab("ME")
  989. p
  990. dev.off()
  991. #
  993. #graph MEs Cerebellum ME turquoise
  994. pdf(file = "boxplot_MEs_FC_green_wide.pdf",wi = 8, he = 6)
  995. #
  996. FCgreen <- datComb2 %>% filter(module == "MEgreen") %>% filter(Brain.Region == "Frontal Cortex")
  998. p <- ggplot(FCgreen, aes(x = Treatment, y = ME, fill = Treatment)) +
  999. facet_wrap(~ Sex.x, ncol = 2) +
  1000. stat_summary(
  1001. geom = "boxplot",
  1002. fun.data = function(x) {
  1003. qs <- quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  1004. setNames(qs, c("ymin", "lower", "middle", "upper", "ymax"))
  1005. },
  1006. position = position_dodge2(width = 0.9, preserve = "single")
  1007. ) +
  1008. geom_point(position = position_dodge(width = 0.9), aes(group = Treatment), alpha = 0.8) +
  1009. scale_fill_manual(values = cbPalette) +
  1010. # two-line x labels: keep "-" on the top line
  1011. scale_x_discrete(labels = function(x) {
  1012. parts <- strsplit(x, "\\s*[\\-–—]\\s*", perl = TRUE)
  1013. sapply(seq_along(parts), function(i) {
  1014. if (length(parts[[i]]) >= 2) paste0(parts[[i]][1], " -\n", parts[[i]][2]) else x[i]
  1015. })
  1016. }) +
  1017. theme_cowplot(font_size = 15) +
  1018. theme(
  1019. axis.text.x = element_text(angle = 0, hjust = 0.5, lineheight = 0.9),
  1020. strip.background = element_rect(fill = "#B0C4DE"),
  1021. strip.text.x = element_text(size = 20),
  1022. legend.position = "none"
  1023. ) +
  1024. xlab("Treatment") +
  1025. ylab("ME")
  1026. p
  1027. dev.off()
  1028. #
  1030. pdf(file = "boxplot_MEs_FC_yellow_wide.pdf", wi = 8, he = 6)
  1031. #
  1032. FCyellow <- datComb2 %>% filter(module == "MEyellow") %>% filter(Brain.Region == "Frontal Cortex")
  1035. p <- ggplot(FCyellow, aes(x = Treatment, y = ME, fill = Treatment)) +
  1036. facet_wrap(~ Sex.x, ncol = 2) +
  1037. stat_summary(
  1038. geom = "boxplot",
  1039. fun.data = function(x) {
  1040. qs <- quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  1041. setNames(qs, c("ymin", "lower", "middle", "upper", "ymax"))
  1042. },
  1043. position = position_dodge2(width = 0.9, preserve = "single")
  1044. ) +
  1045. geom_point(position = position_dodge(width = 0.9), aes(group = Treatment), alpha = 0.8) +
  1046. scale_fill_manual(values = cbPalette) +
  1047. # two-line x labels: keep "-" on the top line
  1048. scale_x_discrete(labels = function(x) {
  1049. parts <- strsplit(x, "\\s*[\\-–—]\\s*", perl = TRUE)
  1050. sapply(seq_along(parts), function(i) {
  1051. if (length(parts[[i]]) >= 2) paste0(parts[[i]][1], " -\n", parts[[i]][2]) else x[i]
  1052. })
  1053. }) +
  1054. theme_cowplot(font_size = 15) +
  1055. theme(
  1056. axis.text.x = element_text(angle = 0, hjust = 0.5, lineheight = 0.9),
  1057. strip.background = element_rect(fill = "#B0C4DE"),
  1058. strip.text.x = element_text(size = 20),
  1059. legend.position = "none"
  1060. ) +
  1061. xlab("Treatment") +
  1062. ylab("ME")
  1064. p
  1065. dev.off()
  1066. #
  1069. pdf(file = "boxplot_MEs_FC_brown_wide.pdf", wi = 8, he = 6)
  1070. #
  1071. FCbrown <- datComb2 %>% filter(module == "MEbrown") %>% filter(Brain.Region == "Frontal Cortex")
  1074. p <- ggplot(FCbrown, aes(x = Treatment, y = ME, fill = Treatment)) +
  1075. facet_wrap(~ Sex.x, ncol = 2) +
  1076. stat_summary(
  1077. geom = "boxplot",
  1078. fun.data = function(x) {
  1079. qs <- quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  1080. setNames(qs, c("ymin", "lower", "middle", "upper", "ymax"))
  1081. },
  1082. position = position_dodge2(width = 0.9, preserve = "single")
  1083. ) +
  1084. geom_point(position = position_dodge(width = 0.9), aes(group = Treatment), alpha = 0.8) +
  1085. scale_fill_manual(values = cbPalette) +
  1086. # two-line x labels: keep "-" on the top line
  1087. scale_x_discrete(labels = function(x) {
  1088. parts <- strsplit(x, "\\s*[\\-–—]\\s*", perl = TRUE)
  1089. sapply(seq_along(parts), function(i) {
  1090. if (length(parts[[i]]) >= 2) paste0(parts[[i]][1], " -\n", parts[[i]][2]) else x[i]
  1091. })
  1092. }) +
  1093. theme_cowplot(font_size = 15) +
  1094. theme(
  1095. axis.text.x = element_text(angle = 0, hjust = 0.5, lineheight = 0.9),
  1096. strip.background = element_rect(fill = "#B0C4DE"),
  1097. strip.text.x = element_text(size = 20),
  1098. legend.position = "none"
  1099. ) +
  1100. xlab("Treatment") +
  1101. ylab("ME")
  1103. p
  1104. dev.off()
  1105. #
  1112. #save.image(file = "AllFiles.RData")
  1113. #load(file = "/Users/aciernia/Sync/collaborations/Ashwood/MIATcellRNAseq/WGCNA/AllFiles.RData")
  1115. #quantify associations of individual genes with our trait of interest (years) by de ning Gene Signi cance GS as
  1116. #(the absolute value of) the correlation between the gene and the trait
  1118. #For each module, we also de ne a quantitative
  1119. #measure of module membership MM as the correlation of the module eigengene and the gene expression pro le. This
  1120. #allows us to quantify the similarity of all genes on the array to every module.
  1123. ##module membership MM
  1124. #correlation between expression of each gene and MEs (uncorrected p values)
  1125. MEs <- MEs %>% dplyr::select(-Sample.ID)
  1126. # names (colors) of the modules
  1127. modNames = substring(names(MEs), 3)
  1129. geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
  1130. MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
  1132. names(geneModuleMembership) = paste("MM", modNames, sep="");
  1133. names(MMPvalue) = paste("p.MM", modNames, sep="");
  1135. #or
  1136. # calculate the module membership values (aka. module eigengene based
  1137. # connectivity kME):
  1138. #datKME = signedKME(PFCdatExpr, MEs)
  1140. #correlate each gene expression with variable of interest
  1141. # Define variable age containing the years column of datTrait
  1142. #dam
  1143. Dam.Treatment = as.data.frame(numerictraits$Dam.Treatment.Group2);
  1144. names(Dam.Treatment) = "Dam.Treatment"
  1146. Dam_geneTraitSignificance = as.data.frame(cor(datExpr, Dam.Treatment, use = "p"));
  1147. Dam_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Dam_geneTraitSignificance), nSamples));
  1149. names(Dam_geneTraitSignificance) = paste("GS.", names(Dam.Treatment), sep="");
  1150. names(Dam_GSPvalue) = paste("p.GS.", names(Dam.Treatment), sep="")
  1152. #pups
  1153. Pup.Treatment = as.data.frame(numerictraits$Pup.Treatment.Group2);
  1154. names(Pup.Treatment) = "Pup.Treatment"
  1156. Pup_geneTraitSignificance = as.data.frame(cor(datExpr, Pup.Treatment, use = "p"));
  1157. Pup_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Pup_geneTraitSignificance), nSamples));
  1159. names(Pup_geneTraitSignificance) = paste("GS.", names(Pup.Treatment), sep="");
  1160. names(Pup_GSPvalue) = paste("p.GS.", names(Pup.Treatment), sep="")
  1162. #Regions
  1163. Region = as.data.frame(numerictraits$Brain.Region2);
  1164. names(Region) = "Brain Region"
  1166. Region_geneTraitSignificance = as.data.frame(cor(datExpr, Region, use = "p"));
  1167. Region_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Region_geneTraitSignificance), nSamples));
  1169. names(Region_geneTraitSignificance) = paste("GS.", names(Region), sep="");
  1170. names(Region_GSPvalue) = paste("p.GS.", names(Region), sep="")
  1172. #=====================================================================================
  1173. #Intramodular analysis: identifying genes with high GS and MM
  1174. #Using the GS and MM measures, we can identify genes that have a high signi cance for weight as well as high module
  1175. #membership in interesting modules.
  1178. #As an example, we look at the brown module that has the highest association with age.
  1180. pdf(file = "Damtreatment_geneModuleMembershipvsGeneSignificanceAGE.pdf", wi = 10, he = 10)
  1182. par(mfrow = c(2,2));
  1184. module = "blue"
  1185. column = match(module, modNames);
  1186. moduleGenes = moduleColors==module;
  1187. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1188. abs(Dam_geneTraitSignificance[moduleGenes, 1]),
  1189. xlab = paste("Module Membership in", module, "module"),
  1190. ylab = "Dam Treatment Gene significance for age",
  1191. main = paste("Module membership vs. Dam treatment gene significance\n"),
  1192. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1195. module = "brown"
  1196. column = match(module, modNames);
  1197. moduleGenes = moduleColors==module;
  1198. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1199. abs(Dam_geneTraitSignificance[moduleGenes, 1]),
  1200. xlab = paste("Module Membership in", module, "module"),
  1201. ylab = "Dam Treatment Gene significance for age",
  1202. main = paste("Module membership vs. Dam treatment gene significance\n"),
  1203. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1205. module = "turquoise"
  1206. column = match(module, modNames);
  1207. moduleGenes = moduleColors==module;
  1208. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1209. abs(Dam_geneTraitSignificance[moduleGenes, 1]),
  1210. xlab = paste("Module Membership in", module, "module"),
  1211. ylab = "Dam Treatment Gene significance for age",
  1212. main = paste("Module membership vs. Dam treatment gene significance\n"),
  1213. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1215. module = "green"
  1216. column = match(module, modNames);
  1217. moduleGenes = moduleColors==module;
  1218. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1219. abs(Dam_geneTraitSignificance[moduleGenes, 1]),
  1220. xlab = paste("Module Membership in", module, "module"),
  1221. ylab = "Dam Treatment Gene significance for age",
  1222. main = paste("Module membership vs. Dam treatment gene significance\n"),
  1223. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1225. module = "yellow"
  1226. column = match(module, modNames);
  1227. moduleGenes = moduleColors==module;
  1228. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1229. abs(Dam_geneTraitSignificance[moduleGenes, 1]),
  1230. xlab = paste("Module Membership in", module, "module"),
  1231. ylab = "Dam Treatment Gene significance for age",
  1232. main = paste("Module membership vs. Dam treatment gene significance\n"),
  1233. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1235. module = "grey"
  1236. column = match(module, modNames);
  1237. moduleGenes = moduleColors==module;
  1238. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1239. abs(Dam_geneTraitSignificance[moduleGenes, 1]),
  1240. xlab = paste("Module Membership in", module, "module"),
  1241. ylab = "Dam Treatment Gene significance for age",
  1242. main = paste("Module membership vs. Dam treatment gene significance\n"),
  1243. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1245. dev.off()
  1249. pdf(file = "Puptreatment_geneModuleMembershipvsGeneSignificanceAGE.pdf", wi = 10, he = 10)
  1251. par(mfrow = c(2,2));
  1253. module = "blue"
  1254. column = match(module, modNames);
  1255. moduleGenes = moduleColors==module;
  1256. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1257. abs(Pup_geneTraitSignificance[moduleGenes, 1]),
  1258. xlab = paste("Module Membership in", module, "module"),
  1259. ylab = "Pup Treatment Gene significance for age",
  1260. main = paste("Module membership vs. Pup treatment gene significance\n"),
  1261. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1264. module = "brown"
  1265. column = match(module, modNames);
  1266. moduleGenes = moduleColors==module;
  1267. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1268. abs(Pup_geneTraitSignificance[moduleGenes, 1]),
  1269. xlab = paste("Module Membership in", module, "module"),
  1270. ylab = "Pup Treatment Gene significance for age",
  1271. main = paste("Module membership vs. Pup treatment gene significance\n"),
  1272. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1274. module = "turquoise"
  1275. column = match(module, modNames);
  1276. moduleGenes = moduleColors==module;
  1277. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1278. abs(Pup_geneTraitSignificance[moduleGenes, 1]),
  1279. xlab = paste("Module Membership in", module, "module"),
  1280. ylab = "Pup Treatment Gene significance for age",
  1281. main = paste("Module membership vs. Pup treatment gene significance\n"),
  1282. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1284. module = "green"
  1285. column = match(module, modNames);
  1286. moduleGenes = moduleColors==module;
  1287. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1288. abs(Pup_geneTraitSignificance[moduleGenes, 1]),
  1289. xlab = paste("Module Membership in", module, "module"),
  1290. ylab = "Pup Treatment Gene significance for age",
  1291. main = paste("Module membership vs. Pup treatment gene significance\n"),
  1292. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1294. module = "yellow"
  1295. column = match(module, modNames);
  1296. moduleGenes = moduleColors==module;
  1297. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1298. abs(Pup_geneTraitSignificance[moduleGenes, 1]),
  1299. xlab = paste("Module Membership in", module, "module"),
  1300. ylab = "Pup Treatment Gene significance for age",
  1301. main = paste("Module membership vs. Pup treatment gene significance\n"),
  1302. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1304. module = "grey"
  1305. column = match(module, modNames);
  1306. moduleGenes = moduleColors==module;
  1307. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1308. abs(Pup_geneTraitSignificance[moduleGenes, 1]),
  1309. xlab = paste("Module Membership in", module, "module"),
  1310. ylab = "Pup Treatment Gene significance for age",
  1311. main = paste("Module membership vs. Pup treatment gene significance\n"),
  1312. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1314. dev.off()
  1317. pdf(file = "Region_geneModuleMembershipvsGeneSignificance.pdf", wi = 10, he = 10)
  1319. par(mfrow = c(2,2));
  1321. module = "blue"
  1322. column = match(module, modNames);
  1323. moduleGenes = moduleColors==module;
  1324. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1325. abs(Region_geneTraitSignificance[moduleGenes, 1]),
  1326. xlab = paste("Module Membership in", module, "module"),
  1327. ylab = "Region Treatment Gene significance for Region",
  1328. main = paste("Module membership vs. Region treatment gene significance\n"),
  1329. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1332. module = "brown"
  1333. column = match(module, modNames);
  1334. moduleGenes = moduleColors==module;
  1335. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1336. abs(Region_geneTraitSignificance[moduleGenes, 1]),
  1337. xlab = paste("Module Membership in", module, "module"),
  1338. ylab = "Region Treatment Gene significance for Region",
  1339. main = paste("Module membership vs. Region treatment gene significance\n"),
  1340. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1342. module = "turquoise"
  1343. column = match(module, modNames);
  1344. moduleGenes = moduleColors==module;
  1345. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1346. abs(Region_geneTraitSignificance[moduleGenes, 1]),
  1347. xlab = paste("Module Membership in", module, "module"),
  1348. ylab = "Region Treatment Gene significance for Region",
  1349. main = paste("Module membership vs. Region treatment gene significance\n"),
  1350. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1352. module = "green"
  1353. column = match(module, modNames);
  1354. moduleGenes = moduleColors==module;
  1355. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1356. abs(Region_geneTraitSignificance[moduleGenes, 1]),
  1357. xlab = paste("Module Membership in", module, "module"),
  1358. ylab = "Region Treatment Gene significance for Region",
  1359. main = paste("Module membership vs. Region treatment gene significance\n"),
  1360. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1362. module = "yellow"
  1363. column = match(module, modNames);
  1364. moduleGenes = moduleColors==module;
  1365. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1366. abs(Region_geneTraitSignificance[moduleGenes, 1]),
  1367. xlab = paste("Module Membership in", module, "module"),
  1368. ylab = "Region Treatment Gene significance for Region",
  1369. main = paste("Module membership vs. Region treatment gene significance\n"),
  1370. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1372. module = "grey"
  1373. column = match(module, modNames);
  1374. moduleGenes = moduleColors==module;
  1375. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  1376. abs(Region_geneTraitSignificance[moduleGenes, 1]),
  1377. xlab = paste("Module Membership in", module, "module"),
  1378. ylab = "Region Treatment Gene significance for Region",
  1379. main = paste("Module membership vs. Region treatment gene significance\n"),
  1380. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  1382. dev.off()
  1385. #Clearly, GS and MM are highly correlated, illustrating that genes highly significantly associated with a trait are
  1386. #often also the most important (central) elements of modules associated with the trait.
  1388. #=====================================================================================
  1390. #We have found modules with high association with our trait of interest, and have identified their central players by the Module Membership measure.
  1391. #We now merge this statistical information with gene annotation and write out a file that summarizes the most important results and can be inspected in standard spreadsheet software
  1392. # Create the starting data frame
  1393. Ensemble = colnames(datExpr)
  1395. #df <- nc %>% dplyr::select(ensembl_gene_id,mgi_symbol,description,genelength) %>% distinct()
  1397. #get matching data for each ensemble id from the row data from Brainspain input
  1398. ids <- genes[which(genes$ensembl_gene_id %in% Ensemble),]
  1400. #match order
  1402. idRows = match(Ensemble, ids$ensembl_gene_id)
  1403. ids2 = ids[idRows, ]
  1406. geneInfo0 = data.frame(ensembl_gene_id = colnames(datExpr),
  1407. geneSymbol = ids2$mgi_symbol,
  1408. description = ids2$description,
  1409. genelength = ids2$genelength,
  1410. moduleColor = moduleColors,
  1411. Dam_geneTraitSignificance,
  1412. Dam_GSPvalue,
  1413. Pup_geneTraitSignificance,
  1414. Pup_GSPvalue,
  1415. Region_geneTraitSignificance,
  1416. Region_GSPvalue,
  1417. t(datExpr))
  1418. # Order modules by their significance for age
  1419. modOrder = order(-abs(cor(MEs, Dam.Treatment, use = "p")));
  1420. # Add module membership information in the chosen order
  1421. for (mod in 1:ncol(geneModuleMembership))
  1422. {
  1423. oldNames = names(geneInfo0)
  1424. geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],
  1425. MMPvalue[, modOrder[mod]]);
  1426. names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
  1427. paste("p.MM.", modNames[modOrder[mod]], sep=""))
  1428. }
  1432. write.csv(geneInfo0, file = "geneInfo_Master_WGCNA.csv")
  1434. #geneInfo0 <- read.csv("geneInfo_BilboMGtimecourse.csv")
  1436. #=====================================================================================
  1437. #output a character vector of genes, where the genes are the hub gene picked for each module,
  1438. #and the names correspond to the module in which each gene is a hub.
  1440. hubs <- chooseTopHubInEachModule(datExpr,moduleColors, type = "signed",power =10) #https://support.bioconductor.org/p/46342/
  1441. #power of 2 for unsigned, 4 for signed
  1442. write.csv(hubs,"TopHubInEachModule.csv")
  1444. #These genes represent the top 10 genes per module based on kME
  1445. topGenesKME = NULL
  1446. for (i in 1:length(colnames(geneModuleMembership))){
  1447. KMErank = geneModuleMembership[(order(-geneModuleMembership[,i])),] #order by column, - decreasing
  1448. KMErank$Ensemble <- rownames(KMErank)
  1449. GenesKME = KMErank[c(1:10),c(i,7)] #where column 7 is the ensemble id
  1451. #get gene info
  1452. topKMEinfo <- ids2[which(ids2$ensembl_gene_id %in% rownames(GenesKME)),]
  1454. #get kME
  1455. merge <- merge(topKMEinfo,GenesKME,by.x="ensembl_gene_id", by.y = "Ensemble")
  1456. colnames(merge)[9] <- c("kME") #module column
  1458. merge$module <- substr(colnames(geneModuleMembership[i]),3,nchar(colnames(geneModuleMembership[i])))
  1460. topGenesKME = rbind(topGenesKME,merge)
  1461. }
  1463. write.csv(topGenesKME,"Top10HubInEachModule_KME.csv")
  1466. #topGenesKME <- read.csv("TopHubInEachModule.csv")
  1468. #=====================================================================================
  1469. #Extract modules
  1470. module_colors = setdiff(unique(moduleColors), "grey")
  1471. # for (color in module_colors){
  1472. # module=datExpr[,moduleColors==module_colors]
  1473. # write.table(module, paste("module_",color, "log2RPKM+1expression_DLPFC.txt",sep=""), sep="\t", row.names=T, col.names=T,quote=FALSE)
  1474. #
  1475. # }
  1477. #Look at expression patterns of these genes, as they are clustered
  1478. #heatmap colors:
  1479. #install.packages("RColorBrewer")
  1480. library("RColorBrewer")
  1481. library(gplots)
  1483. #gene-module info
  1485. heatDF <- geneInfo0 %>% filter(moduleColor != "grey") %>% arrange(moduleColor)
  1486. heatDF$moduleColor <- factor(heatDF$moduleColor)
  1488. data <- as.matrix(heatDF[,12:188]) #RPKM values
  1489. rownames(data) <- heatDF$ensembl_gene_id
  1492. #data <- t(m)
  1493. library(pheatmap)
  1494. my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 299)
  1496. #define column groups:
  1497. annotation_col <- data.frame(
  1498. #Sample = factor(d$donor_name),
  1499. Dam.Treatment = factor(numerictraits$Dam.Treatment.Group2, levels = unique(numerictraits$Dam.Treatment.Group2)),
  1500. Pup.Treatment = factor(numerictraits$Pup.Treatment.Group2, levels = unique(numerictraits$Pup.Treatment.Group2)),
  1501. Region = factor(numerictraits$Brain.Region2, levels = unique(numerictraits$Brain.Region2)),
  1502. Sex = factor(numerictraits$sex2))
  1504. rownames(annotation_col) = colnames(data)
  1506. head(annotation_col)
  1508. annotation_row <- data.frame(heatDF$moduleColor)
  1510. rownames(annotation_row) = rownames(data)
  1512. colnames(annotation_row) <- c("module")
  1514. head(annotation_row)
  1516. # Specify colors
  1517. #https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
  1518. library(RColorBrewer)
  1519. n <- 27
  1520. qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
  1521. col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
  1522. #pie(rep(1,n), col=sample(col_vector, n))
  1525. # change the color of annotation to what you want: (eg: "navy", "darkgreen")
  1526. Var1 <- sample(col_vector, 2)
  1527. names(Var1) <- unique(annotation_col$Dam.Treatment)
  1529. Var2 <- sample(col_vector, 2)
  1530. names(Var2) <- unique(annotation_col$Pup.Treatment)
  1532. Var3 <- c("cornflowerblue","darksalmon")
  1533. names(Var3) <- unique(annotation_col$Sex)
  1535. Var4 <- sample(col_vector, 3)
  1536. names(Var4) <- unique(annotation_col$Region)
  1538. Var5 <- as.character(unique(heatDF$moduleColor))
  1539. names(Var5) <- unique(heatDF$moduleColor)
  1541. anno_colors <- list(Dam.Treatment = Var1, Pup.Treatment = Var2, Sex = Var3, Brain.Region = Var4, moduleColor = Var5)
  1543. pdf(file = "module-Heatmap_log2RPKM_zscore.pdf", wi = 20, he = 16)
  1544. pheatmap::pheatmap(data,
  1545. cluster_row = F,
  1546. cluster_cols = T,
  1547. annotation_col = annotation_col,
  1548. annotation_row = annotation_row,
  1549. color = my_palette,
  1550. fontsize = 10,
  1551. fontsize_row=6,
  1552. show_rownames = F,
  1553. fontsize_col = 10,
  1554. annotation_colors = anno_colors,
  1555. scale = c("row"),
  1556. main = "Gene Expression Modules")
  1558. dev.off()
  1560. #=====================================================================================
  1561. #average heatmap by condition
  1562. heatDF2 <- heatDF[,c(1,5,12:188)] #ensembl ID, modulecolour, log2PRKM values
  1564. heatDF2 <- heatDF2 %>% gather(Sample.ID, log2RPKM, 3:ncol(heatDF2))
  1565. headDF3 <- merge(heatDF2,datTraits, by = "Sample.ID", all.x=T)
  1568. heatDF4 <- headDF3 %>% group_by(moduleColor,ensembl_gene_id,Group) %>% summarize(meanlog2RPKM= mean(log2RPKM)) %>%
  1569. spread(Group,meanlog2RPKM)%>% arrange(moduleColor)
  1571. matrix <- as.matrix(heatDF4[,3:ncol(heatDF4)])
  1572. rownames(matrix) <- heatDF4$ensembl_gene_id
  1574. #order columns
  1575. colnames(matrix)
  1576. #newcolorder <- c( "7 F" , "7 M" ,"15 F", "15 M", "35 F" ,"35 M" , "84 F", "84 M")
  1578. #matrix <- matrix[,newcolorder]
  1579. #colnames(matrix)
  1581. my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 299)
  1583. #define column groups:
  1584. Dam.Treatment <- substr(colnames(matrix),1,6)
  1585. tmp <- as.data.frame(colnames(matrix))
  1586. colnames(tmp) <- c("group")
  1587. tmp <- tmp %>% separate(group,into=c("Dam","Pup","sex","region"),sep="_")
  1589. sex <- substr(colnames(matrix),(nchar(colnames(matrix))+1)-1,nchar(colnames(matrix)))
  1591. annotation_col <- tmp
  1593. rownames(annotation_col) = colnames(matrix)
  1595. head(annotation_col)
  1597. annotation_row <- data.frame(heatDF4$moduleColor)
  1599. rownames(annotation_row) = rownames(matrix)
  1601. colnames(annotation_row) <- c("moduleColor")
  1603. head(annotation_row)
  1605. # Specify colors
  1606. #https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
  1607. library(RColorBrewer)
  1608. n <- 27
  1609. qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
  1610. col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
  1611. #pie(rep(1,n), col=sample(col_vector, n))
  1614. # change the color of annotation to what you want: (eg: "navy", "darkgreen")
  1615. Var1 <- sample(col_vector, 2)
  1616. names(Var1) <- unique(annotation_col$Dam)
  1618. Var2 <- sample(col_vector, 2)
  1619. names(Var2) <- unique(annotation_col$Pup)
  1621. Var3 <- c("cornflowerblue","darksalmon")
  1622. names(Var3) <- unique(annotation_col$sex)
  1624. Var4 <- sample(col_vector, 3)
  1625. names(Var4) <- unique(annotation_col$region)
  1627. Var5 <- as.character(unique(heatDF4$moduleColor))
  1628. names(Var5) <- unique(heatDF4$moduleColor)
  1630. anno_colors <- list(DamTreatment = Var1, PupTreatment = Var2, Sex = Var3, BrainRegion = Var4, moduleColor = Var5)
  1632. pdf(file = "Heatmap_log2RPKM_zscore_noReplicates.pdf", wi = 11, he = 8.5)
  1633. pheatmap::pheatmap(matrix,
  1634. cluster_row = F,
  1635. cluster_cols = T,
  1636. annotation_col = annotation_col,
  1637. annotation_row = annotation_row,
  1638. color = my_palette,
  1639. show_rownames = F,
  1640. fontsize = 12,
  1641. fontsize_row=6,
  1642. fontsize_col = 12,
  1643. annotation_colors = anno_colors,
  1644. scale = c("row"),
  1645. main = "Gene Expression Modules")
  1647. dev.off()
  1650. #=====================================================================================
  1651. #split by brain region
  1652. # #=====================================================================================
  1654. ##############Cerebellum###################
  1655. CB <- headDF3 %>% filter(Brain.Region == "Cerebellum")
  1657. heatDF4 <- CB %>% group_by(moduleColor,ensembl_gene_id,Group) %>% summarize(meanlog2RPKM= mean(log2RPKM)) %>%
  1658. spread(Group,meanlog2RPKM)%>% arrange(moduleColor)
  1660. matrix <- as.matrix(heatDF4[,3:ncol(heatDF4)])
  1661. rownames(matrix) <- heatDF4$ensembl_gene_id
  1663. #order columns
  1664. colnames(matrix)
  1665. newcolorder <- c("Saline_Saline_M_Cerebellum", "Saline_Saline_F_Cerebellum",
  1666. "Saline_Treg_M_Cerebellum","Saline_Treg_F_Cerebellum",
  1667. "PolyIC_Saline_M_Cerebellum","PolyIC_Saline_F_Cerebellum",
  1668. "PolyIC_Treg_M_Cerebellum", "PolyIC_Treg_F_Cerebellum")
  1670. matrix <- matrix[,newcolorder]
  1671. colnames(matrix)
  1673. my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 299)
  1675. #define column groups:
  1676. tmp <- as.data.frame(colnames(matrix))
  1677. colnames(tmp) <- c("group")
  1678. tmp <- tmp %>% separate(group,into=c("Dam","Pup","sex","region"),sep="_")
  1680. sex <- substr(colnames(matrix),(nchar(colnames(matrix))+1)-1,nchar(colnames(matrix)))
  1682. annotation_col <- tmp
  1684. rownames(annotation_col) = colnames(matrix)
  1686. head(annotation_col)
  1688. annotation_row <- data.frame(heatDF4$moduleColor)
  1690. rownames(annotation_row) = rownames(matrix)
  1692. colnames(annotation_row) <- c("moduleColor")
  1694. head(annotation_row)
  1696. # Specify colors
  1697. #https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
  1698. library(RColorBrewer)
  1699. n <- 27
  1700. qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
  1701. col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
  1702. #pie(rep(1,n), col=sample(col_vector, n))
  1705. # change the color of annotation to what you want: (eg: "navy", "darkgreen")
  1706. Var1 <- sample(col_vector, 2)
  1707. names(Var1) <- unique(annotation_col$Dam)
  1709. Var2 <- sample(col_vector, 2)
  1710. names(Var2) <- unique(annotation_col$Pup)
  1712. Var3 <- c("cornflowerblue","darksalmon")
  1713. names(Var3) <- unique(annotation_col$sex)
  1715. Var4 <- sample(col_vector, 3)
  1716. names(Var4) <- unique(annotation_col$region)
  1718. Var5 <- as.character(unique(heatDF4$moduleColor))
  1719. names(Var5) <- unique(heatDF4$moduleColor)
  1721. anno_colors <- list(DamTreatment = Var1, PupTreatment = Var2, Sex = Var3, BrainRegion = Var4, moduleColor = Var5)
  1723. pdf(file = "CB_Heatmap_log2RPKM_zscore_noReplicates.pdf", wi = 11, he = 8.5)
  1724. pheatmap::pheatmap(matrix,
  1725. cluster_row = F,
  1726. cluster_cols = F,
  1727. annotation_col = annotation_col,
  1728. annotation_row = annotation_row,
  1729. color = my_palette,
  1730. show_rownames = F,
  1731. fontsize = 12,
  1732. fontsize_row=6,
  1733. fontsize_col = 12,
  1734. annotation_colors = anno_colors,
  1735. scale = c("row"),
  1736. main = "Gene Expression Modules")
  1738. dev.off()
  1740. ##############Frontal Cortex###################
  1741. FC <- headDF3 %>% filter(Brain.Region == "Frontal Cortex")
  1743. heatDF4 <- FC %>% group_by(moduleColor,ensembl_gene_id,Group) %>% summarize(meanlog2RPKM= mean(log2RPKM)) %>%
  1744. spread(Group,meanlog2RPKM)%>% arrange(moduleColor)
  1746. matrix <- as.matrix(heatDF4[,3:ncol(heatDF4)])
  1747. rownames(matrix) <- heatDF4$ensembl_gene_id
  1749. #order columns
  1750. #order columns
  1751. colnames(matrix)
  1752. newcolorder <- c("Saline_Saline_M_Frontal Cortex", "Saline_Saline_F_Frontal Cortex",
  1753. "Saline_Treg_M_Frontal Cortex","Saline_Treg_F_Frontal Cortex",
  1754. "PolyIC_Saline_M_Frontal Cortex","PolyIC_Saline_F_Frontal Cortex",
  1755. "PolyIC_Treg_M_Frontal Cortex", "PolyIC_Treg_F_Frontal Cortex")
  1757. matrix <- matrix[,newcolorder]
  1758. colnames(matrix)
  1760. my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 299)
  1762. #define column groups:
  1763. tmp <- as.data.frame(colnames(matrix))
  1764. colnames(tmp) <- c("group")
  1765. tmp <- tmp %>% separate(group,into=c("Dam","Pup","sex","region"),sep="_")
  1767. sex <- substr(colnames(matrix),(nchar(colnames(matrix))+1)-1,nchar(colnames(matrix)))
  1769. annotation_col <- tmp
  1771. rownames(annotation_col) = colnames(matrix)
  1773. head(annotation_col)
  1775. annotation_row <- data.frame(heatDF4$moduleColor)
  1777. rownames(annotation_row) = rownames(matrix)
  1779. colnames(annotation_row) <- c("moduleColor")
  1781. head(annotation_row)
  1783. # Specify colors
  1784. #https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
  1785. library(RColorBrewer)
  1786. n <- 27
  1787. qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
  1788. col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
  1789. #pie(rep(1,n), col=sample(col_vector, n))
  1792. # change the color of annotation to what you want: (eg: "navy", "darkgreen")
  1793. Var1 <- sample(col_vector, 2)
  1794. names(Var1) <- unique(annotation_col$Dam)
  1796. Var2 <- sample(col_vector, 2)
  1797. names(Var2) <- unique(annotation_col$Pup)
  1799. Var3 <- c("cornflowerblue","darksalmon")
  1800. names(Var3) <- unique(annotation_col$sex)
  1802. Var4 <- sample(col_vector, 3)
  1803. names(Var4) <- unique(annotation_col$region)
  1805. Var5 <- as.character(unique(heatDF4$moduleColor))
  1806. names(Var5) <- unique(heatDF4$moduleColor)
  1808. anno_colors <- list(DamTreatment = Var1, PupTreatment = Var2, Sex = Var3, BrainRegion = Var4, moduleColor = Var5)
  1810. pdf(file = "FC_Heatmap_log2RPKM_zscore_noReplicates.pdf", wi = 11, he = 8.5)
  1811. pheatmap::pheatmap(matrix,
  1812. cluster_row = F,
  1813. cluster_cols = F,
  1814. annotation_col = annotation_col,
  1815. annotation_row = annotation_row,
  1816. color = my_palette,
  1817. show_rownames = F,
  1818. fontsize = 12,
  1819. fontsize_row=6,
  1820. fontsize_col = 12,
  1821. annotation_colors = anno_colors,
  1822. scale = c("row"),
  1823. main = "Gene Expression Modules")
  1825. dev.off()
  1827. ##############FHippocampus###################
  1828. HC <- headDF3 %>% filter(Brain.Region == "Hippocampus")
  1830. heatDF4 <- HC %>% group_by(moduleColor,ensembl_gene_id,Group) %>% summarize(meanlog2RPKM= mean(log2RPKM)) %>%
  1831. spread(Group,meanlog2RPKM)%>% arrange(moduleColor)
  1833. matrix <- as.matrix(heatDF4[,3:ncol(heatDF4)])
  1834. rownames(matrix) <- heatDF4$ensembl_gene_id
  1836. #order columns
  1837. colnames(matrix)
  1838. newcolorder <- c("Saline_Saline_M_Hippocampus", "Saline_Saline_F_Hippocampus",
  1839. "Saline_Treg_M_Hippocampus","Saline_Treg_F_Hippocampus",
  1840. "PolyIC_Saline_M_Hippocampus","PolyIC_Saline_F_Hippocampus",
  1841. "PolyIC_Treg_M_Hippocampus", "PolyIC_Treg_F_Hippocampus")
  1843. matrix <- matrix[,newcolorder]
  1844. colnames(matrix)
  1847. my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 299)
  1849. #define column groups:
  1850. tmp <- as.data.frame(colnames(matrix))
  1851. colnames(tmp) <- c("group")
  1852. tmp <- tmp %>% separate(group,into=c("Dam","Pup","sex","region"),sep="_")
  1855. annotation_col <- tmp
  1857. rownames(annotation_col) = colnames(matrix)
  1859. head(annotation_col)
  1861. annotation_row <- data.frame(heatDF4$moduleColor)
  1863. rownames(annotation_row) = rownames(matrix)
  1865. colnames(annotation_row) <- c("moduleColor")
  1867. head(annotation_row)
  1869. # Specify colors
  1870. #https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
  1871. library(RColorBrewer)
  1872. n <- 27
  1873. qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
  1874. col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
  1875. #pie(rep(1,n), col=sample(col_vector, n))
  1878. # change the color of annotation to what you want: (eg: "navy", "darkgreen")
  1879. Var1 <- sample(col_vector, 2)
  1880. names(Var1) <- unique(annotation_col$Dam)
  1882. Var2 <- sample(col_vector, 2)
  1883. names(Var2) <- unique(annotation_col$Pup)
  1885. Var3 <- c("cornflowerblue","darksalmon")
  1886. names(Var3) <- unique(annotation_col$sex)
  1888. Var4 <- sample(col_vector, 3)
  1889. names(Var4) <- unique(annotation_col$region)
  1891. Var5 <- as.character(unique(heatDF4$moduleColor))
  1892. names(Var5) <- unique(heatDF4$moduleColor)
  1894. anno_colors <- list(DamTreatment = Var1, PupTreatment = Var2, Sex = Var3, BrainRegion = Var4, moduleColor = Var5)
  1896. pdf(file = "HC_Heatmap_log2RPKM_zscore_noReplicates.pdf", wi = 11, he = 8.5)
  1897. pheatmap::pheatmap(matrix,
  1898. cluster_row = F,
  1899. cluster_cols = F,
  1900. annotation_col = annotation_col,
  1901. annotation_row = annotation_row,
  1902. color = my_palette,
  1903. show_rownames = F,
  1904. fontsize = 12,
  1905. fontsize_row=6,
  1906. fontsize_col = 12,
  1907. annotation_colors = anno_colors,
  1908. scale = c("row"),
  1909. main = "Gene Expression Modules")
  1911. dev.off()
  1914. #=====================================================================================
  1915. #line plots
  1916. # #=====================================================================================
  1918. # # #scale data same as pheatmap: subtract mean and divide by std. dev
  1919. # # #https://www.biostars.org/p/223532/
  1920. # scale_rows = function(x){
  1921. # m = apply(x, 1, mean, na.rm = T)
  1922. # s = apply(x, 1, sd, na.rm = T)
  1923. # return((x - m) / s)
  1924. # }
  1925. # #
  1926. # scale_mat = function(mat, scale){
  1927. # if(!(scale %in% c("none", "row", "column"))){
  1928. # stop("scale argument shoud take values: 'none', 'row' or 'column'")
  1929. # }
  1930. # mat = switch(scale, none = mat, row = scale_rows(mat), column = t(scale_rows(t(mat))))
  1931. # return(mat)
  1932. # }
  1933. # #
  1934. # # #data:
  1935. # heatDF4 <- headDF3 %>% group_by(moduleColor,ensembl_gene_id,Group) %>% summarize(meanlog2exp = mean(log2RPKM)) %>%
  1936. # spread(Group,meanlog2exp)%>% arrange(moduleColor)
  1937. # #
  1938. # matrix <- as.matrix(heatDF4[,3:ncol(heatDF4)])
  1939. # rownames(matrix) <- heatDF4$ensembl_gene_id
  1940. # #
  1941. # #
  1942. # # #scale by row:
  1943. # merPFCmatrix <- scale_mat(matrix, scale = "row")
  1944. # merPFCmatrix <- as.data.frame(merPFCmatrix)
  1945. # #
  1946. # merPFCmatrix$module <- heatDF4$moduleColor
  1947. # #
  1948. # merPFCDF <- merPFCmatrix %>% gather(group, scaled_expression, 1:(ncol(merPFCmatrix)-1)) %>%
  1949. # group_by(module,group)
  1950. #%>%
  1951. # # summarize(meanscaled_expression = mean(scaled_expression))
  1952. # #
  1953. # merPFCDF <- as.data.frame(merPFCDF)
  1954. # #
  1955. # # #arrange
  1956. # # newcolorder <- c( "7 F" , "7 M" ,"15 F", "15 M", "35 F" ,"35 M" , "84 F", "84 M")
  1957. # #
  1958. # # merPFCDF$age <- factor(merPFCDF$age,levels = newcolorder)
  1959. # #
  1960. # # #set colors:
  1961. # my_palette2 <- unique(merPFCDF$module)
  1962. # #
  1963. # my_palette2 <- col2hex(my_palette2)
  1964. # #
  1965. # library(cowplot)
  1966. # pdf(file = "lineplot_log2RPKM_scaled_noReplicates.pdf", wi = 11, he = 8.5)
  1967. # #
  1968. # p <- ggplot(merPFCDF, aes(x=group, y=meanscaled_expression, group=module,colour=module)) +
  1969. # geom_line(size=2) +
  1970. # scale_color_manual(values = my_palette2)+
  1971. # theme(axis.text.x = element_text(angle=45, vjust=0.5))+
  1972. # xlab(label = c("Treatment")) +
  1973. # ylab(label = c("Scaled log2(RPKM+1)"))
  1974. # #
  1975. # p
  1976. # #
  1977. # dev.off()
  1978. #
  1979. # #facet by module
  1981. merPFCDF2 <- merPFCDF %>% separate(group, into=c("DamTreatment","PupTreatment","sex","region"), sep="_")
  1982. merPFCDF2$DamTreatment <- factor( merPFCDF2$DamTreatment,levels =c("Saline","PolyIC" ))
  1984. theme_set(theme_classic(base_size = 24))
  1985. cbPalette <- c("lightblue","pink")
  1987. pdf(file = "boxplot_log2RPKM+1_scaled_noReplicates_facet.pdf", wi = 18, he = 30)
  1988. #
  1989. p <- ggplot(merPFCDF2, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  1990. facet_wrap(~region*module*sex,ncol = 5) +
  1991. stat_summary(geom = "boxplot",
  1992. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  1993. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  1994. # geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  1995. scale_fill_manual(values = cbPalette) +
  1996. theme_cowplot(font_size = 15)+
  1997. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  1998. strip.background = element_rect(fill="#B0C4DE"),
  1999. # legend.position="none",
  2000. strip.text.x = element_text(size = 30))+
  2001. xlab(label = c("Dam Treatment")) +
  2002. ylab(label = c("Scaled log2(RPKM)"))
  2003. p
  2004. dev.off()
  2005. #
  2007. pdf(file = "boxplot_log2RPKM+1_scaled_noReplicates_facet2.pdf", wi = 22, he = 30)
  2008. #
  2009. p <- ggplot(merPFCDF2, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  2010. facet_wrap(~module*region*sex,ncol = 6) +
  2011. stat_summary(geom = "boxplot",
  2012. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  2013. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  2014. #geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  2015. scale_fill_manual(values = cbPalette) +
  2016. theme_cowplot(font_size = 15)+
  2017. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  2018. strip.background = element_rect(fill="#B0C4DE"),
  2019. # legend.position="none",
  2020. strip.text.x = element_text(size = 20))+
  2021. xlab(label = c("Dam Treatment")) +
  2022. ylab(label = c("Scaled log2(RPKM)"))
  2023. p
  2024. dev.off()
  2025. #
  2027. pdf(file = "boxplot_log2RPKM+1_scaled_CB_bluemodule.pdf", wi = 6, he = 4)
  2028. #
  2029. CBblueRPKM <- merPFCDF2 %>% filter(module == "blue") %>% filter(region == "Cerebellum")
  2030. p <- ggplot( CBblueRPKM, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  2031. facet_wrap(~sex,ncol = 2) +
  2033. stat_summary(geom = "boxplot",
  2034. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  2035. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  2036. #geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  2037. scale_fill_manual(values = cbPalette) +
  2038. theme_cowplot(font_size = 15)+
  2039. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  2040. strip.background = element_rect(fill="#B0C4DE"),
  2041. # legend.position="none",
  2042. strip.text.x = element_text(size = 20))+
  2043. xlab(label = c("Dam Treatment")) +
  2044. ylab(label = c("Scaled log2(RPKM)"))
  2045. p
  2046. dev.off()
  2048. pdf(file = "boxplot_log2RPKM+1_scaled_CB_turquoisemodule.pdf", wi = 6, he = 4)
  2049. #
  2050. CBturquoiseRPKM <- merPFCDF2 %>% filter(module == "turquoise") %>% filter(region == "Cerebellum")
  2051. p <- ggplot( CBturquoiseRPKM, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  2052. facet_wrap(~sex,ncol = 2) +
  2054. stat_summary(geom = "boxplot",
  2055. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  2056. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  2057. #geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  2058. scale_fill_manual(values = cbPalette) +
  2059. theme_cowplot(font_size = 15)+
  2060. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  2061. strip.background = element_rect(fill="#B0C4DE"),
  2062. # legend.position="none",
  2063. strip.text.x = element_text(size = 20))+
  2064. xlab(label = c("Dam Treatment")) +
  2065. ylab(label = c("Scaled log2(RPKM)"))
  2066. p
  2067. dev.off()
  2069. pdf(file = "boxplot_log2RPKM+1_scaled_HC_turquoisemodule.pdf", wi = 6, he = 4)
  2070. #
  2071. HCturquoiseRPKM <- merPFCDF2 %>% filter(module == "turquoise") %>% filter(region == "Hippocampus")
  2072. p <- ggplot( HCturquoiseRPKM, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  2073. facet_wrap(~sex,ncol = 2) +
  2075. stat_summary(geom = "boxplot",
  2076. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  2077. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  2078. #geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  2079. scale_fill_manual(values = cbPalette) +
  2080. theme_cowplot(font_size = 15)+
  2081. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  2082. strip.background = element_rect(fill="#B0C4DE"),
  2083. # legend.position="none",
  2084. strip.text.x = element_text(size = 20))+
  2085. xlab(label = c("Dam Treatment")) +
  2086. ylab(label = c("Scaled log2(RPKM)"))
  2087. p
  2088. dev.off()
  2090. pdf(file = "boxplot_log2RPKM+1_scaled_FC_Greenmodule.pdf", wi = 6, he = 4)
  2091. #
  2092. FCgreenRPKM <- merPFCDF2 %>% filter(module == "green") %>% filter(region == "Frontal Cortex")
  2093. p <- ggplot(FCgreenRPKM, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  2094. facet_wrap(~sex,ncol = 2) +
  2096. stat_summary(geom = "boxplot",
  2097. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  2098. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  2099. #geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  2100. scale_fill_manual(values = cbPalette) +
  2101. theme_cowplot(font_size = 15)+
  2102. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  2103. strip.background = element_rect(fill="#B0C4DE"),
  2104. # legend.position="none",
  2105. strip.text.x = element_text(size = 20))+
  2106. xlab(label = c("Dam Treatment")) +
  2107. ylab(label = c("Scaled log2(RPKM)"))
  2108. p
  2109. dev.off()
  2111. pdf(file = "boxplot_log2RPKM+1_scaled_FC_yellowmodule.pdf", wi = 6, he = 4)
  2112. #
  2113. FCyellowRPKM <- merPFCDF2 %>% filter(module == "yellow") %>% filter(region == "Frontal Cortex")
  2114. p <- ggplot(FCyellowRPKM, aes(x=DamTreatment, y=scaled_expression), group = PupTreatment) +
  2115. facet_wrap(~sex,ncol = 2) +
  2117. stat_summary(geom = "boxplot",
  2118. fun.data = function(x) setNames(quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95)), c("ymin", "lower", "middle", "upper", "ymax")),
  2119. position = "dodge", aes(fill=PupTreatment))+ #gives line as median, top of box and bottom of box as 25th and 75th quartile, line = 5th and 9th percentile
  2120. #geom_point(position = position_dodge(width = 0.90),aes(group=PupTreatment, colour=sex)) +
  2121. scale_fill_manual(values = cbPalette) +
  2122. theme_cowplot(font_size = 15)+
  2123. theme(axis.text.x = element_text(angle=45, vjust=0.5),
  2124. strip.background = element_rect(fill="#B0C4DE"),
  2125. # legend.position="none",
  2126. strip.text.x = element_text(size = 20))+
  2127. xlab(label = c("Dam Treatment")) +
  2128. ylab(label = c("Scaled log2(RPKM)"))
  2129. p
  2130. dev.off()
  2133. #=====================================================================================
  2134. #GO term enrichment:https://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/GeneAnnotation/
  2135. #source("http://bioconductor.org/biocLite.R");
  2136. #biocLite(c("impute", "AnnotationDBI", "GO.db", "org.Hs.eg.db", "org.Mm.eg.db"))
  2138. #install.packages("anRichment_1.01-2.tar.gz")
  2139. #source("https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/GeneAnnotation/installAnRichment.R");
  2140. # installAnRichment();
  2142. library(anRichment)
  2143. options(stringsAsFactors = FALSE)
  2146. ###################################################################
  2147. GOcollection = buildGOcollection(organism = "mouse")
  2149. #get entrez IDs:
  2150. entrez <- getBM(attributes = c("ensembl_gene_id", "entrezgene_id"),
  2151. filter= "ensembl_gene_id",
  2152. values = geneInfo0$ensembl_gene_id,
  2153. mart = mouse)
  2155. geneInfo <- merge(geneInfo0, entrez, by="ensembl_gene_id")
  2157. # evaluates the enrichment of the gene modules in the collection of GO terms
  2159. GOenrichment = enrichmentAnalysis(
  2160. classLabels = geneInfo$moduleColor, identifiers = geneInfo$entrezgene_id,
  2161. refCollection = GOcollection,
  2162. useBackground = "intersection", #intersection of genes in identifiers and the organism database
  2163. threshold = 0.05,
  2164. thresholdType = "FDR",
  2165. getOverlapEntrez = FALSE,
  2166. getOverlapSymbols = TRUE,
  2167. entrySeparator = ",",
  2168. maxReportedOverlapGenes = 1000,
  2169. ignoreLabels = "grey") #ignore genes not assigned to modules (grey)
  2170. #
  2171. collectGarbage()
  2172. #
  2173. names(GOenrichment)
  2174. #
  2175. #enrichment results are summarized in the component enrichmentTable:
  2176. names(GOenrichment$enrichmentTable)
  2177. #
  2179. write.csv(GOenrichment$enrichmentTable, file = "GOenrichment-enrichmentTable.csv")
  2180. #
  2181. #

WGCNA_allconditions.R at commit 432ab07, no license · at the source

Overview

Authors: Destanie Rose1,2, Rachel J Moreno1,2, Hadley Osman1,2, Megan E Rowland3,4, Jill Silverman2,5, Annie Ciernia3,4, Paul Ashwood1,2
  1. Department of Microbiology and Immunology, UC Davis, 50th Street Sacramento, Davis, CA 95817 USA
  2. The M.I.N.D. Institute, University of California at Davis, Davis, CA USA
  3. Department of Biochemistry and Molecular Biology, University of British Columbia, Vancouver, BC Canada
  4. Djavad Mowafaghian Centre for Brain Health, University of British Columbia, Vancouver, BC Canada
  5. Deperatment of Psychiatry and Behavioral Sciences, University of California at Davis, Davis, CA USA
Institutions: University of California, Davis (United States); University of British Columbia (Canada); Djavad Mowafaghian Centre for Brain Health (Canada)
Journal: Journal of neuroinflammation, volume 23, issue 1, article 133
Dates: received 20 December 2025; accepted 12 February 2026; published online 12 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1186/s12974-026-03739-w · PMID 41821031 · PMCID PMC13097898 · OpenAlex W7135079076
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), autism (population), schizophrenia / psychosis (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning
Keywords: Maternal immune activation (MIA), Autism, Schizophrenia, Immune regulation, Regulatory T cells (Tregs), Inflammation, Behavior, TH17 cells, Sex, Gene expression
MeSH: Adoptive Transfer*, Autism Spectrum Disorder*, Brain*, Prenatal Exposure Delayed Effects*, Sex Characteristics*, T-Lymphocytes, Regulatory*, Animals, Cytokines, Disease Models, Animal, Female, Male, Mice, Mice, Inbred C57BL, Pregnancy (* major topic)
Topic: Autism Spectrum Disorder Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: National Institute of Child Health and Human Development
Citations: cited by 1 paper (Europe PMC); 67 references in the paper

Abstract

Autism spectrum disorder (ASD) is characterized by atypical communication, social interactions, and restricted interests. In ASD, there are dysfunctional immune regulatory control mechanisms that can lead to immune activation. Notably lower frequencies of regulatory T cells (Tregs) and reduced immunosuppressive cytokines are reported and associated with more impaired behaviors impacting both individuals with autism and their families. Therefore, therapeutic approaches that enhance immune regulation may offer substantial benefits. Using the maternal immune activation (MIA) model, we investigated whether adoptive transfer of wildtype Tregs into MIA offspring recipients could rescue immune activation, brain transcriptome changes and behaviors exhibited in adult MIA offspring. We also aimed to explore potential sex-differences in responses. In male but not female MIA offspring, Tregs transfer reduced the frequency of T helper (TH)-17 (RORγT⁺ CD4⁺) T cells in both the mesenteric lymph node (MLN) and spleen. Moreover, the frequency of CD25⁺Foxp3⁺ T cells was increased in the MLN and spleen of male but not female MIA offspring following Tregs transfer. Splenocytes from male MIA offspring receiving Tregs showed reduced production of inflammatory cytokines (e.g., IL-6 and TNFα) following PMA/Ionomycin stimulations. In contrast, female MIA offspring that received Tregs exhibited different cytokine profiles characterized by increased production of cytokines, including GM-CSF, IFNγ, and IL-10. In the brain, bulk mRNA sequencing in the cerebellum, frontal cortex, and hippocampus revealed that Tregs-treated male MIA offspring had differentially expressed genes involved in neurodevelopmental disorders, synaptic function, and epigenetic regulation. Minimal gene expression differences were observed in female counterparts. There was significant improvement in self-grooming behaviors in males MIA offspring that received Tregs. In females, social novelty improved after Tregs treatment. In summary, adoptive Tregs transfer reduced systemic inflammation, brain transcription and behavior alterations in a sex dependent manner in the context of MIA. These findings suggest that adoptive Tregs transfer may represent a viable therapeutic avenue for mitigating systemic inflammation and comorbidities associated with MIA and neurodevelopmental disorders such as ASD.

Supplementary Information: The online version contains supplementary material available at 10.1186/s12974-026-03739-w.

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

Repository

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

ciernialab/Rose-et-al.-2026

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 432ab0775a42fb94f6fff3bff53ab2574a85dc4a, 16 November 2025
Languages: R (6), Shell (2)
Size: 10 files, 8 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: cowplot (6 files), ggplot2 (6 files), tidyverse (6 files), clusterProfiler (5 files), edgeR (5 files), limma (5 files), pheatmap (5 files), SAMtools (2 files), FastQC (1 file), nlme (1 file), STAR (1 file), Subread (featureCounts) (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
9 files
At the source: github.com/ciernialab/Rose-et-al.-2026

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

Tracing map

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

What the map holds:

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

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

Data

Datasets cited

Data availability

All processed gene lists and statistics are available in Supplementary Tables. Raw fastq files and counts matrix for the RNA-seq are available on NCBI GEO at X. Analysis code is on our lab github ( [https://github.com/ciernialab/Rose-et-al.-2026](https://github.com/ciernialab/Rose-et-al.-2026).

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 10 keywords, 14 MeSH terms, 1 funder, 65 references.

Cite

This paper

Rose, D., Moreno, R. J., Osman, H., Rowland, M. E., Silverman, J., Ciernia, A., & Ashwood, P. (2026). Sex specific effects of adoptive Tregs transfer on the brain and periphery in maternal immune activation offspring rescuing immune dysregulation. Journal of neuroinflammation, 23(1), 133. https://doi.org/10.1186/s12974-026-03739-w

BibTeX

@article{rose2026sex,
author = {Rose, Destanie and Moreno, Rachel J and Osman, Hadley and Rowland, Megan E and Silverman, Jill and Ciernia, Annie and Ashwood, Paul},
title = {{Sex specific effects of adoptive Tregs transfer on the brain and periphery in maternal immune activation offspring rescuing immune dysregulation}},
journal = {Journal of neuroinflammation},
year = {2026},
month = mar,
volume = {23},
number = {1},
pages = {133},
publisher = {BMC},
issn = {1742-2094},
doi = {10.1186/s12974-026-03739-w},
url = {https://doi.org/10.1186/s12974-026-03739-w},
pmid = {41821031},
pmcid = {PMC13097898}
}

RIS

TY - JOUR
AU - Rose, Destanie
AU - Moreno, Rachel J
AU - Osman, Hadley
AU - Rowland, Megan E
AU - Silverman, Jill
AU - Ciernia, Annie
AU - Ashwood, Paul
TI - Sex specific effects of adoptive Tregs transfer on the brain and periphery in maternal immune activation offspring rescuing immune dysregulation
T2 - Journal of neuroinflammation
J2 - J Neuroinflammation
PY - 2026
DA - 2026/03/12
VL - 23
IS - 1
SP - 133
SN - 1742-2094
PB - BMC
DO - 10.1186/s12974-026-03739-w
UR - https://doi.org/10.1186/s12974-026-03739-w
LA - en
ER -

CSL-JSON

{
"id": "10.1186/s12974-026-03739-w",
"type": "article-journal",
"title": "Sex specific effects of adoptive Tregs transfer on the brain and periphery in maternal immune activation offspring rescuing immune dysregulation",
"container-title": "Journal of neuroinflammation",
"author": [
{
"family": "Rose",
"given": "Destanie"
},
{
"family": "Moreno",
"given": "Rachel J"
},
{
"family": "Osman",
"given": "Hadley"
},
{
"family": "Rowland",
"given": "Megan E"
},
{
"family": "Silverman",
"given": "Jill"
},
{
"family": "Ciernia",
"given": "Annie"
},
{
"family": "Ashwood",
"given": "Paul"
}
],
"container-title-short": "J Neuroinflammation",
"volume": "23",
"issue": "1",
"page": "133",
"DOI": "10.1186/s12974-026-03739-w",
"PMID": "41821031",
"PMCID": "PMC13097898",
"ISSN": "1742-2094",
"publisher": "BMC",
"URL": "https://doi.org/10.1186/s12974-026-03739-w",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
12
]
]
}
}

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.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: FastQC, Subread (featureCounts), STAR, 9 other tools
[2] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Subread (featureCounts), nlme, WGCNA, 7 other tools, mouse, cellular / molecular
[3] doi:10.21203/rs.3.rs-9927928/v1 [code]
Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: FastQC, Subread (featureCounts), STAR, 6 other tools, mouse
[4] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: Subread (featureCounts), STAR, SAMtools, 6 other tools, mouse, cellular / molecular
[5] doi:10.1038/s42003-026-10059-5 [code]
Transcriptomic analysis in autism spectrum disorder suggests three molecular subtypes with distinct phenotypic profiles and functional pathways.
Journal: Communications biology
In common: FastQC, STAR, WGCNA, 4 other tools, autism, cellular / molecular, 1 reference
[6] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: STAR, WGCNA, SAMtools, 6 other tools, autism
[7] doi:10.7554/elife.109104 [code]
Complementary vertebrate &lt;i&gt;Wac&lt;/i&gt; models exhibit phenotypes relevant to DeSanto-Shinawi Syndrome.
Journal: eLife
In common: FastQC, Subread (featureCounts), STAR, 5 other tools, mouse, cellular / molecular
[8] doi:10.1038/s41467-026-71432-w [code]
MeCP2 gene dosage-dependent neurodevelopmentally restricted defects arise by aberrant activation of cell fate-determining bivalent genes.
Journal: Nature communications
In common: FastQC, Subread (featureCounts), STAR, 5 other tools, mouse, cellular / molecular
[9] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Subread (featureCounts), WGCNA, edgeR, 6 other tools, cellular / molecular
[10] doi:10.1038/s41386-026-02406-1 [code]
Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators.
Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology
In common: FastQC, STAR, SAMtools, 5 other tools, schizophrenia / psychosis, cellular / molecular

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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