OSCR

Sex and tissue resolved co-expression networks reveal a female placental-brain axis protective against prenatal PCB exposure.

Code ↔ Paper

2 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 2 matches
  1. [1] § Methods › Bioinformatic Analyses › Weighted gene correlation network analysis ↔ scripts/ConsensusBrainComparisonPCBonly.Rmd, lines 213–269 · score 0.59 · Scale free topology, model fit, connectivity, soft, threshold, power
  2. [2] § Methods › Bioinformatic Analyses › Weighted gene correlation network analysis ↔ scripts/ConsensusPlacentaComparisonPCBonly.Rmd, lines 187–242 · score 0.59 · Scale free topology, model fit, connectivity, soft, threshold, power

Paper

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

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 1,077 lines · 31 KB · MIT · 1 match

  1. ---
  2. title: "ConsensusBrainEnrichRComparisons"
  3. output: html_document
  4. date: "2023-05-22"
  5. ---
  7. ```{r setup, include=FALSE}
  8. knitr::opts_chunk$set(echo = TRUE)
  9. ```
  12. ```{r}
  13. library(WGCNA)
  14. library(readxl)
  15. #took out samples that had Folic Acid
  16. F_Brain <- read_excel("/Users/kchau/Documents/PEBBLES/RNA-seq_ComparisonBrain_PCBonly/F_Brain_voom_PCBonly.xlsx")
  17. F_Brain <- F_Brain[!grepl("PLAC", colnames(F_Brain))]
  19. M_Brain <- read_excel("/Users/kchau/Documents/PEBBLES/RNA-seq_ComparisonBrain_PCBonly/M_Brain_voom_PCBonly.xlsx")
  20. M_Brain <- M_Brain[!grepl("PLAC", colnames(M_Brain))]
  22. # Take a quick look at what is in the data sets (caution, longish output):
  23. dim(F_Brain)
  24. names(F_Brain)
  25. dim(M_Brain)
  26. names(M_Brain)
  27. ```
  30. ```{r}
  31. intersectedgenes <- intersect(rownames(F_Brain), rownames(M_Brain))
  32. F_Brain1 <- F_Brain[intersectedgenes, ]
  33. M_Brain1 <- M_Brain[intersectedgenes, ]
  34. ```
  37. ```{r}
  38. # We work with two sets:
  39. nSets = 2;
  40. # For easier labeling of plots, create a vector holding descriptive names of the two sets.
  41. setLabels = c("Female brain", "Male brain")
  42. shortLabels = c("Female", "Male")
  43. # Form multi-set expression data: columns starting from 9 contain actual expression data.
  44. multiExpr = vector(mode = "list", length = nSets)
  45. multiExpr[[1]] = list(data = as.data.frame(t(F_Brain1[-c(1:1)])));
  46. names(multiExpr[[1]]$data) = F_Brain1$Gene;
  47. rownames(multiExpr[[1]]$data) = names(F_Brain1)[-c(1:1)];
  48. multiExpr[[2]] = list(data = as.data.frame(t(M_Brain1[-c(1:1)])));
  49. names(multiExpr[[2]]$data) = M_Brain1$Gene;
  50. rownames(multiExpr[[2]]$data) = names(M_Brain1)[-c(1:1)];
  51. ```
  54. ```{r}
  55. # Check that the data has the correct format for many functions operating on multiple sets:
  56. exprSize = checkSets(multiExpr)
  57. ```
  59. ```{r}
  60. # Check that all genes and samples have sufficiently low numbers of missing values.
  61. gsg = goodSamplesGenesMS(multiExpr, verbose = 3);
  62. gsg$allOK
  63. ```
  65. ```{r}
  66. if (!gsg$allOK)
  67. {
  68. # Print information about the removed genes:
  69. if (sum(!gsg$goodGenes) > 0)
  70. printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]$data)[!gsg$goodGenes],
  71. collapse = ", ")))
  72. for (set in 1:exprSize$nSets)
  73. {
  75. if (sum(!gsg$goodSamples[[set]]))
  76. printFlush(paste("In set", setLabels[set], "removing samples",
  77. paste(rownames(multiExpr[[set]]$data)[!gsg$goodSamples[[set]]], collapse = ", ")))
  78. # Remove the offending genes and samples
  79. multiExpr[[set]]$data = multiExpr[[set]]$data[gsg$goodSamples[[set]], gsg$goodGenes];
  80. }
  81. # Update exprSize
  82. exprSize = checkSets(multiExpr)
  83. }
  84. ```
  87. ```{r}
  88. sampleTrees = list()
  89. for (set in 1:nSets)
  90. {
  91. sampleTrees[[set]] = hclust(dist(multiExpr[[set]]$data), method = "average")
  92. }
  93. ```
  95. ```{r}
  96. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/SamplesClustering.pdf", width = 12, height = 12);
  97. par(mfrow=c(2,1))
  98. par(mar = c(0, 4, 2, 0))
  99. for (set in 1:nSets)
  100. plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
  101. xlab="", sub="", cex = 0.7);
  102. dev.off()
  103. ```
  106. ```{r}
  107. # Choose the "base" cut height for the female data set
  108. baseHeight = 60
  109. # Adjust the cut height for the male data set for the number of samples
  110. cutHeights = c(60, 60*exprSize$nSamples[2]/exprSize$nSamples[1]);
  111. #cutHeights = c(60, 100);
  112. # Re-plot the dendrograms including the cut lines
  113. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/SamplesClusteringCutoff.pdf", width = 12, height = 12);
  114. par(mfrow=c(2,1))
  115. par(mar = c(0, 4, 2, 0))
  116. for (set in 1:nSets)
  117. {
  118. plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
  119. xlab="", sub="", cex = 0.7);
  120. abline(h=cutHeights[set], col = "red");
  121. }
  122. dev.off()
  123. ```
  125. ```{r}
  126. #outlier removal
  127. for (set in 1:nSets)
  128. {
  129. #Find clusters cut by the line
  130. labels = cutreeStatic(sampleTrees[[set]], cutHeight = cutHeights[set], minSize = 10)
  131. #Keep the largest one (labeled by the number 1)
  132. keep = (labels==1)
  133. multiExpr[[set]]$data = multiExpr[[set]]$data[keep, ]
  134. }
  136. labels1 <- cutreeStatic(sampleTrees[[1]], cutHeight = 60)
  137. labels2 <- cutreeStatic(sampleTrees[[2]], cutHeight = 100, minSize = 10)
  138. ```
  140. ```{r}
  141. collectGarbage();
  142. # Check the size of the leftover data
  143. exprSize = checkSets(multiExpr)
  144. exprSize
  145. ```
  148. ```{r}
  149. PCBtraitData = read_excel("/Users/kchau/Documents/PEBBLES/sampletraits/M_F_Brains copy_combinations_PCB_only.xls")
  150. #FAtraitData = read_excel("/Users/kchau/Documents/PEBBLES/RNA-seq/RawData/FolicAcidM_F_Brainscopy_combinations.xls")
  151. #traitData <- traitData[-c(97:999),]
  152. PCBtraitData <- PCBtraitData[!grepl("PLAC", PCBtraitData$Name),]
  153. dim(PCBtraitData)
  154. names(PCBtraitData)
  155. ```
  157. ```{r}
  158. #allTraits = traitData[, c(1, 5, 7, 8)]
  159. PCBallTraits = PCBtraitData[, -c(2,3,4, 6)]
  160. PCBallTraits = PCBallTraits[-c(59:96),] #remove samples that had Folic Acid added
  161. #FAallTraits = FAtraitData[, -c(6)]
  162. #allTraits = PCBtraitData[, -c(4)]
  163. #allTraits = traitData
  164. # See how big the traits are and what are the trait and sample names
  165. dim(PCBallTraits)
  166. names(PCBallTraits)
  167. PCBallTraits$Name
  168. ```
  170. ```{r}
  171. #allTraits$Folic_Acid <- as.integer(as.factor(allTraits$Folic_Acid))
  172. #library(plyr)
  173. #FAallTraits$Folic_Acid <- revalue(FAallTraits$Folic_Acid, c("YES"=1))
  174. #FAallTraits$Folic_Acid <- revalue(FAallTraits$Folic_Acid, c("NO"=0))
  175. #allTraits$Sex <- revalue(allTraits$Sex, c("F"=0))
  176. #allTraits$Sex <- revalue(allTraits$Sex, c("M"=1))
  177. ```
  180. ```{r}
  181. #remove female brain outlier: 36_01; there are no male brain outliers
  182. #allTraits <- allTraits[!grepl("35_05_BRAIN", allTraits$Name), ]
  183. #allTraits <- allTraits[!grepl("35_06_BRAIN", allTraits$Name), ]
  184. #PCBallTraits <- PCBallTraits[!grepl("36_01_BRAIN", PCBallTraits$Name), ]
  185. #FAallTraits <- FAallTraits[!grepl("36_01_BRAIN", FAallTraits$Name), ]
  186. ```
  188. ```{r}
  189. PCBallTraits <- PCBallTraits[ , -c(2, 4)]
  190. ```
  194. ```{r}
  195. #PCB
  196. # Form a multi-set structure that will hold the sample traits.
  197. PCBTraits = vector(mode="list", length = nSets);
  198. for (set in 1:nSets)
  199. {
  200. PCBsetSamples = rownames(multiExpr[[set]]$data);
  201. PCBtraitRows = match(PCBsetSamples, PCBallTraits$Name);
  202. PCBTraits[[set]] = list(data = PCBallTraits[PCBtraitRows, -1]);
  203. rownames(PCBTraits[[set]]$data) = as.data.frame(PCBallTraits)[PCBtraitRows, 1]
  204. #rownames(PCBTraits[[set]]$data) = as.data.frame(PCBallTraits)[PCBtraitRows, 1];
  205. }
  206. collectGarbage();
  207. # Define data set dimensions
  208. PCBnGenes = exprSize$nGenes;
  209. PCBnSamples = exprSize$nSamples;
  210. ```
  213. ```{r}
  214. # Choose a set of soft-thresholding powers
  215. powers = c(seq(4,10,by=1), seq(12,40, by=2));
  216. # Initialize a list to hold the results of scale-free analysis
  217. powerTables = vector(mode = "list", length = nSets);
  218. # Call the network topology analysis function for each set in turn
  219. for (set in 1:nSets)
  220. powerTables[[set]] = list(data = pickSoftThreshold(multiExpr[[set]]$data, powerVector=powers,
  221. verbose = 2)[[2]]);
  222. collectGarbage();
  223. # Plot the results:
  224. colors = c("black", "red")
  225. # Will plot these columns of the returned scale free analysis tables
  226. plotCols = c(2,5,6,7)
  227. colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "Median connectivity",
  228. "Max connectivity");
  229. # Get the minima and maxima of the plotted points
  230. ylim = matrix(NA, nrow = 2, ncol = 4);
  231. for (set in 1:nSets)
  232. {
  233. for (col in 1:length(plotCols))
  235. {
  236. ylim[1, col] = min(ylim[1, col], powerTables[[set]]$data[, plotCols[col]], na.rm = TRUE);
  237. ylim[2, col] = max(ylim[2, col], powerTables[[set]]$data[, plotCols[col]], na.rm = TRUE);
  238. }
  239. }
  240. # Plot the quantities in the chosen columns vs. the soft thresholding power
  241. sizeGrWindow(8, 6)
  242. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/scaleFreeAnalysis.pdf", wi = 8, he = 6)
  243. par(mfcol = c(2,2));
  244. par(mar = c(4.2, 4.2 , 2.2, 0.5))
  245. cex1 = 0.7;
  246. for (col in 1:length(plotCols)) for (set in 1:nSets)
  247. {
  248. if (set==1)
  249. {
  250. plot(powerTables[[set]]$data[,1], -sign(powerTables[[set]]$data[,3])*powerTables[[set]]$data[,2],
  251. xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],
  252. main = colNames[col]);
  253. addGrid();
  254. }
  255. if (col==1)
  256. {
  257. text(powerTables[[set]]$data[,1], -sign(powerTables[[set]]$data[,3])*powerTables[[set]]$data[,2],
  258. labels=powers,cex=cex1,col=colors[set]);
  259. } else
  260. text(powerTables[[set]]$data[,1], powerTables[[set]]$data[,plotCols[col]],
  261. labels=powers,cex=cex1,col=colors[set]);
  262. if (col==1)
  263. {
  264. legend("bottomright", legend = setLabels, col = colors, pch = 20) ;
  265. } else
  266. legend("topright", legend = setLabels, col = colors, pch = 20) ;
  267. }
  268. dev.off();
  269. ```
  272. #increase minModuleSize to minimize number of modules
  273. ```{r}
  274. #change back to 50 min module size
  275. bnet = blockwiseConsensusModules(
  276. multiExpr, maxBlockSize = PCBnGenes, power = 22, minModuleSize = 50,
  277. networkType = "signed",
  278. TOMType = "signed",
  279. deepSplit = 2,
  280. pamRespectsDendro = FALSE,
  281. mergeCutHeight = 0.25, numericLabels = TRUE,
  282. minKMEtoStay = 0,
  283. saveTOMs = TRUE, verbose = 5)
  284. ```
  288. ```{r}
  289. PCBconsMEs = bnet$multiMEs;
  290. PCBmoduleLabels = bnet$colors;
  291. # Convert the numeric labels to color labels
  292. PCBmoduleColors = labels2colors(PCBmoduleLabels)
  293. PCBconsTree = bnet$dendrograms[[1]]
  294. ```
  297. ```{r}
  298. sizeGrWindow(12,6)
  300. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/PCBBlockwiseGeneDendrosAndColors.pdf", wi = 12, he = 6);
  301. # Use the layout function for more involved screen sectioning
  302. layout(matrix(c(1:4), 2, 2), heights = c(0.8, 0.2), widths = c(1,1))
  303. #layout.show(4);
  304. nBlocks = length(bnet$dendrograms)
  305. # Plot the dendrogram and the module colors underneath for each block
  306. for (block in 1:nBlocks)
  307. plotDendroAndColors(bnet$dendrograms[[block]], PCBmoduleColors[bnet$blockGenes[[block]]],
  308. "Module colors",
  309. main = paste("Consensus Brain Gene dendrogram and module colors in block", block),
  310. dendroLabels = FALSE, hang = 0.03,
  311. addGuide = TRUE, guideHang = 0.05,
  312. setLayout = FALSE)
  313. dev.off()
  314. ```
  318. ```{r}
  319. #find number of modules
  320. table(bnet$colors)
  321. ```
  323. ```{r}
  324. library(dplyr)
  325. PCBTraits[[set]]$data <- PCBTraits[[set]]$data %>%
  326. dplyr::select(-Sex, -LitterCoded)
  327. ```
  329. ```{r}
  330. # Set up variables to contain the module-trait correlations
  331. PCBmoduleTraitCor = list();
  332. PCBmoduleTraitPvalue = list();
  333. # Calculate the correlations
  334. for (set in 1:nSets)
  335. {
  336. PCBmoduleTraitCor[[set]] = cor(PCBconsMEs[[set]]$data, PCBTraits[[set]]$data, use = "p", method = 'spearman');
  337. PCBmoduleTraitPvalue[[set]] = corPvalueFisher(PCBmoduleTraitCor[[set]], exprSize$nSamples[set])
  338. }
  339. ```
  342. ```{r}
  343. # Convert numerical lables to colors for labeling of modules in the plot
  344. PCBMEColors = labels2colors(as.numeric(substring(names(PCBconsMEs[[1]]$data), 3)));
  345. PCBMEColorNames = paste("ME", PCBMEColors, sep="");
  346. # Open a suitably sized window (the user should change the window size if necessary)
  347. sizeGrWindow(10,7)
  348. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/hublitterConsensusModuleTraitRelationships-PCBfemale.pdf", wi = 12, he = 9);
  349. # Plot the module-trait relationship table for set number 1
  350. set = 1
  351. textMatrix = paste(signif(PCBmoduleTraitCor[[set]], 2), "\n(",
  352. signif(PCBmoduleTraitPvalue[[set]], 1), ")", sep = "");
  353. dim(textMatrix) = dim(PCBmoduleTraitCor[[set]])
  354. par(mar = c(6, 8.8, 3, 2.2));
  355. labeledHeatmap(Matrix = PCBmoduleTraitCor[[set]],
  356. xLabels = names(PCBTraits[[set]]$data),
  357. yLabels = PCBMEColorNames,
  358. ySymbols = PCBMEColorNames,
  359. colorLabels = FALSE,
  360. colors = greenWhiteRed(50),
  361. textMatrix = textMatrix,
  362. setStdMargins = FALSE,
  363. cex.text = 0.7,
  364. zlim = c(-1,1),
  365. main = paste("ConsensusModule PCB trait relationships in", setLabels[set]))
  366. dev.off();
  367. # Plot the module-trait relationship table for set number 2
  368. set = 2
  369. textMatrix = paste(signif(PCBmoduleTraitCor[[set]], 2), "\n(",
  370. signif(PCBmoduleTraitPvalue[[set]], 1), ")", sep = "");
  371. dim(textMatrix) = dim(PCBmoduleTraitCor[[set]])
  372. sizeGrWindow(10,7)
  373. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/hublitterConsensusModuleTraitRelationships-PCBmale.pdf", wi = 12, he = 9);
  374. par(mar = c(6, 8.8, 3, 2.2));
  375. labeledHeatmap(Matrix = PCBmoduleTraitCor[[set]],
  376. xLabels = names(PCBTraits[[set]]$data),
  377. yLabels = PCBMEColorNames,
  378. ySymbols = PCBMEColorNames,
  379. colorLabels = FALSE,
  380. colors = greenWhiteRed(50),
  381. textMatrix = textMatrix,
  382. setStdMargins = FALSE,
  383. cex.text = 0.7,
  384. zlim = c(-1,1),
  385. main = paste("ConsensusModule PCB Trait Relationships in", setLabels[set]))
  386. dev.off()
  387. ```
  390. ```{r}
  391. # Initialize matrices to hold the consensus correlation and p-value
  392. PCBconsensusCor = matrix(NA, nrow(PCBmoduleTraitCor[[1]]), ncol(PCBmoduleTraitCor[[1]]));
  393. PCBconsensusPvalue = matrix(NA, nrow(PCBmoduleTraitCor[[1]]), ncol(PCBmoduleTraitCor[[1]]));
  394. # Find consensus negative correlations
  395. PCBnegative = PCBmoduleTraitCor[[1]] < 0 & PCBmoduleTraitCor[[2]] < 0;
  396. PCBconsensusCor[PCBnegative] = pmax(PCBmoduleTraitCor[[1]][PCBnegative], PCBmoduleTraitCor[[2]][PCBnegative]);
  397. PCBconsensusPvalue[PCBnegative] = pmax(PCBmoduleTraitPvalue[[1]][PCBnegative], PCBmoduleTraitPvalue[[2]][PCBnegative]);
  398. # Find consensus positive correlations
  399. PCBpositive = PCBmoduleTraitCor[[1]] > 0 & PCBmoduleTraitCor[[2]] > 0;
  400. PCBconsensusCor[PCBpositive] = pmin(PCBmoduleTraitCor[[1]][PCBpositive], PCBmoduleTraitCor[[2]][PCBpositive]);
  401. PCBconsensusPvalue[PCBpositive] = pmax(PCBmoduleTraitPvalue[[1]][PCBpositive], PCBmoduleTraitPvalue[[2]][PCBpositive])
  402. ```
  404. ```{r}
  405. textMatrix = paste(signif(PCBconsensusCor, 2), "\n(",
  406. signif(PCBconsensusPvalue, 1), ")", sep = "");
  407. dim(textMatrix) = dim(PCBmoduleTraitCor[[set]])
  408. sizeGrWindow(10,7)
  409. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/ModuleTraitRelationships-PCBconsensus.pdf", wi = 12, he = 9);
  410. par(mar = c(6, 8.8, 3, 2.2));
  411. labeledHeatmap(Matrix = PCBconsensusCor,
  412. xLabels = names(PCBTraits[[set]]$data),
  413. yLabels = PCBMEColorNames,
  414. ySymbols = PCBMEColorNames,
  415. colorLabels = FALSE,
  416. colors = greenWhiteRed(50),
  417. textMatrix = textMatrix,
  418. setStdMargins = FALSE,
  419. cex.text = 0.7,
  420. zlim = c(-1,1),
  421. main = paste("Consensus module PCB trait relationships across\n",
  422. paste(setLabels, collapse = " and ")))
  423. ```
  427. ```{r}
  428. # Create a variable weight that will hold just the body weight of mice in both sets
  429. PCBDose = vector(mode = "list", length = nSets);
  430. for (set in 1:nSets)
  431. {
  432. PCBDose[[set]] = list(data = as.data.frame(PCBTraits[[set]]$data$PCB_dose));
  433. names(PCBDose[[set]]$data) = "PCB Dose"
  434. }
  435. # Recalculate consMEs to give them color names
  436. consMEsC = multiSetMEs(multiExpr, universalColors = PCBmoduleColors);
  437. # We add the weight trait to the eigengenes and order them by consesus hierarchical clustering:
  438. MET = consensusOrderMEs(addTraitToMEs(consMEsC, PCBDose))
  439. ```
  441. ```{r}
  442. sizeGrWindow(8,10);
  443. pdf(file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/PCBEigengeneNetworks.pdf", width= 10, height = 12);
  444. par(cex = 0.9)
  445. plotEigengeneNetworks(MET, setLabels, marDendro = c(0,2,2,1), marHeatmap = c(3,3,2,1),
  446. zlimPreservation = c(0.5, 1), xLabelsAngle = 90)
  447. dev.off()
  448. ```
  451. #HUB GENE
  453. ```{r}
  454. hub <- chooseTopHubInEachModule(
  455. datExpr,
  456. bnet$colors,
  457. omitColors = "grey",
  458. power = 2,
  459. type = "signed")
  460. ```
  462. ```{r}
  463. moduleLookup <- data.frame(
  464. ModuleNumber = sort(unique(bnet$colors)),
  465. ModuleColor = labels2colors(sort(unique(bnet$colors)))
  466. )
  468. moduleLookup
  469. ```
  471. ```{r}
  472. hub_df <- data.frame(
  473. ModuleNumber = as.numeric(names(hub)),
  474. HubGene = unname(hub)
  475. )
  477. hub_df <- merge(hub_df, moduleLookup, by = "ModuleNumber")
  478. ```
  480. ```{r}
  481. ?chooseOneHubInEachModule
  482. ```
  485. ```{r}
  486. Onehubgene <- chooseOneHubInEachModule(
  487. datExpr,
  488. bnet$colors,
  489. numGenes = 100,
  490. omitColors = "grey",
  491. power = 2,
  492. type = "signed")
  494. ```
  496. ```{r}
  497. onehubmoduleLookup <- data.frame(
  498. ModuleNumber = sort(unique(bnet$colors)),
  499. ModuleColor = labels2colors(sort(unique(bnet$colors)))
  500. )
  502. onehubmoduleLookup
  503. ```
  505. ```{r}
  506. onehub_df <- data.frame(
  507. ModuleNumber = as.numeric(names(Onehubgene)),
  508. HubGene = unname(Onehubgene)
  509. )
  511. onehub_df <- merge(onehub_df, onehubmoduleLookup, by = "ModuleNumber")
  512. ```
  516. ```{r}
  517. library(biomaRt)
  519. ensembl <- useEnsembl(
  520. biomart = "genes",
  521. dataset = "mmusculus_gene_ensembl",
  522. mirror = "useast"
  523. )
  525. annot <- getBM(
  526. attributes = c("ensembl_gene_id", "mgi_symbol"),
  527. filters = "ensembl_gene_id",
  528. values = hub_df$HubGene,
  529. mart = ensembl
  530. )
  532. hub_df <- merge(
  533. hub_df,
  534. annot,
  535. by.x = "HubGene",
  536. by.y = "ensembl_gene_id",
  537. all.x = TRUE
  538. )
  540. ```
  542. ```{r}
  544. ensembl <- useEnsembl(
  545. biomart = "genes",
  546. dataset = "mmusculus_gene_ensembl",
  547. mirror = "useast"
  548. )
  550. annot <- getBM(
  551. attributes = c("ensembl_gene_id", "mgi_symbol"),
  552. filters = "ensembl_gene_id",
  553. values = onehub_df$HubGene,
  554. mart = ensembl
  555. )
  557. onehub_df <- merge(
  558. onehub_df,
  559. annot,
  560. by.x = "HubGene",
  561. by.y = "ensembl_gene_id",
  562. all.x = TRUE
  563. )
  565. ```
  567. ```{r}
  568. hub_kME <- sapply(colnames(MEs), function(ME) {
  569. module <- sub("^ME", "", ME) # e.g. "blue"
  570. genesInModule <- moduleColors == module
  572. rownames(kME)[genesInModule][
  573. which.max(abs(kME[genesInModule, ME]))
  574. ]
  575. })
  576. ```
  578. ```{r}
  579. MET <- consensusModuleEigengenes(
  580. multiExpr,
  581. colors = moduleColors,
  582. excludeGrey = TRUE
  583. )$eigengenes
  584. ```
  586. ```{r}
  587. consensusColors <- bnet$colors
  588. moduleColors <- consensusColors
  589. moduleColors <- labels2colors(moduleColors)
  591. nSets <- length(multiExpr)
  593. colorsMat <- matrix(
  594. moduleColors,
  595. nrow = length(moduleColors),
  596. ncol = nSets
  597. )
  599. colnames(colorsMat) <- paste0("set", 1:nSets)
  601. colorsList <- list(
  602. moduleColors,
  603. moduleColors
  604. )
  607. MET <- multiSetMEs(
  608. multiExpr,
  609. colors = moduleColors,
  610. excludeGrey = TRUE
  611. )
  612. ```
  615. ```{r}
  616. datExpr <- multiExpr[[1]]$data
  617. kMEtable <- signedKME(multiExpr, MET)
  619. ```
  622. ```{r}
  623. kME1 <- kME[[1]]$cor
  624. kME2 <- kME[[2]]$cor
  626. moduleLabels <- as.numeric(factor(moduleColors))
  629. hubGenes <- sapply(colnames(kME1), function(ME) {
  631. moduleLabel <- as.numeric(sub("^ME", "", ME)) # 1, 2, ...
  633. genesInModule <- moduleLabels == moduleLabel
  634. geneNames <- rownames(kME1)[genesInModule]
  636. meanKME <- (abs(kME1[genesInModule, ME]) +
  637. abs(kME2[genesInModule, ME])) / 2
  639. geneNames[which.max(meanKME)]
  640. })
  642. hubGenes
  643. ```
  649. ```{r}
  650. library("AnnotationDbi")
  651. library("org.Mm.eg.db")
  652. #columns(org.Mm.eg.db) # returns list of available keytypes
  653. F_Brain1$entrez = mapIds(org.Mm.eg.db,
  654. keys=F_Brain1$Gene, #Column containing Ensembl gene ids
  655. column="ENTREZID",
  656. keytype="ENSEMBL",
  657. multiVals="first")
  658. ```
  662. ```{r}
  663. #columns(org.Mm.eg.db) # returns list of available keytypes
  664. M_Brain1$entrez = mapIds(org.Mm.eg.db,
  665. keys=M_Brain1$Gene, #Column containing Ensembl gene ids
  666. column="ENTREZID",
  667. keytype="ENSEMBL",
  668. multiVals="first")
  669. ```
  672. ```{r}
  673. #GO Terms Female Brain
  674. GOenrFemaleBrain = GOenrichmentAnalysis(PCBmoduleColors, F_Brain1$entrez, organism = "mouse", evidence = "all", nBestP = 10)
  676. #GO Terms Male Brain
  677. GOenrMaleBrain = GOenrichmentAnalysis(PCBmoduleColors, M_Brain1$entrez, organism = "mouse", evidence = "all", nBestP = 10)
  678. ```
  682. ```{r}
  683. tabFemaleBrain = GOenrFemaleBrain$bestPTerms[[4]]$enrichment
  684. tabMaleBrain = GOenrMaleBrain$bestPTerms[[4]]$enrichment
  685. ```
  687. ```{r}
  688. write.table(tabFemaleBrain, file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/ConsensusFemaleBrainGOEnrichmentTable.csv", sep = ",", quote = TRUE, row.names = FALSE)
  690. write.table(tabMaleBrain, file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/ConsensusMaleBrainGOEnrichmentTable.csv", sep = ",", quote = TRUE, row.names = FALSE)
  691. ```
  697. #For Consensus GO terms
  698. ```{r}
  699. consMEs.unord = multiSetMEs(multiExpr, universalColors = PCBmoduleLabels, excludeGrey = TRUE)
  700. GS = list();
  701. kME = list();
  702. for (set in 1:nSets)
  703. {
  704. GS[[set]] = corAndPvalue(multiExpr[[set]]$data, PCBTraits[[set]]$data);
  705. kME[[set]] = corAndPvalue(multiExpr[[set]]$data, consMEs.unord[[set]]$data);
  706. }
  707. ```
  709. ```{r}
  710. GS.metaZ = (GS[[1]]$Z + GS[[2]]$Z)/sqrt(2);
  711. kME.metaZ = (kME[[1]]$Z + kME[[2]]$Z)/sqrt(2);
  712. GS.metaP = 2*pnorm(abs(GS.metaZ), lower.tail = FALSE);
  713. kME.metaP = 2*pnorm(abs(kME.metaZ), lower.tail = FALSE)
  714. ```
  716. ```{r}
  717. GSmat = rbind(GS[[1]]$cor, GS[[2]]$cor, GS[[1]]$p, GS[[2]]$p, GS.metaZ, GS.metaP);
  718. nTraits = checkSets(PCBTraits)$nGenes
  719. traitNames = colnames(PCBTraits[[1]]$data)
  720. dim(GSmat) = c(PCBnGenes, 6*nTraits)
  721. #rownames(GSmat) = probes;
  722. rownames(GSmat) = M_Brain1$Gene;
  723. colnames(GSmat) = spaste(
  724. c("GS.set1.", "GS.set2.", "p.GS.set1.", "p.GS.set2.", "Z.GS.meta.", "p.GS.meta"),
  725. rep(traitNames, rep(6, nTraits)))
  726. # Same code for kME:
  727. kMEmat = rbind(kME[[1]]$cor, kME[[2]]$cor, kME[[1]]$p, kME[[2]]$p, kME.metaZ, kME.metaP);
  728. MEnames = colnames(consMEs.unord[[1]]$data);
  729. nMEs = checkSets(consMEs.unord)$nGenes
  730. dim(kMEmat) = c(PCBnGenes, 6*nMEs)
  731. #rownames(kMEmat) = probes;
  732. rownames(kMEmat) = M_Brain1$Gene;
  733. colnames(kMEmat) = spaste(
  734. c("kME.set1.", "kME.set2.", "p.kME.set1.", "p.kME.set2.", "Z.kME.meta.", "p.kME.meta"),
  735. rep(MEnames, rep(6, nMEs)))
  736. ```
  738. ```{r}
  739. info = data.frame(Probe = M_Brain1$Gene, EntrezID = M_Brain1$entrez,
  740. ModuleLabel = PCBmoduleLabels,
  741. ModuleColor = labels2colors(PCBmoduleLabels), GSmat, kMEmat);
  742. write.csv(info, file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/testingbrain_consensusAnalysis-CombinedNetworkResults.csv",
  743. row.names = FALSE, quote = FALSE)
  744. ```
  746. ```{r}
  747. info <- read.csv("/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/testingbrain_consensusAnalysis-CombinedNetworkResults.csv")
  748. ```
  751. ```{r}
  752. hubTable <- do.call(rbind, lapply(sort(unique(info$ModuleLabel)), function(m) {
  754. colname <- paste0("Z.kME.meta.ME", m)
  756. # skip modules without eigengenes (e.g. grey)
  757. if (!colname %in% colnames(info)) return(NULL)
  759. sub <- info[info$ModuleLabel == m, ]
  761. # force numeric (guards against factor/character conversion)
  762. zvals <- as.numeric(sub[[colname]])
  764. sub[which.max(abs(zvals)), ]
  765. }))
  767. ```
  770. ```{r}
  771. nrow(hubTable)
  772. length(unique(hubTable$ModuleLabel))
  773. any(hubTable$ModuleColor == "grey")
  774. ```
  779. ```{r}
  780. # kME1 is the correlation matrix: genes × MEs
  781. ME2color <- sapply(colnames(kME1), function(ME) {
  783. # genes whose strongest membership is this ME
  784. genesInModule <- apply(abs(kME1), 1, function(x)
  785. colnames(kME1)[which.max(x)] == ME
  786. )
  788. # most common color among those genes
  789. names(sort(table(moduleColors[genesInModule]), decreasing = TRUE))[1]
  790. })
  792. ME2color
  793. ```
  797. ```{r}
  798. info <- info[, c(1:4)]
  799. info$GeneSymbol <- 'NA'
  801. #Convert entrezID to gene symbols using biomart
  802. library(biomaRt)
  803. library(dplyr)
  804. library(xlsx)
  807. mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
  808. #got the probe and gene symbol that matches
  809. geneSymbols <- getBM(filters = "ensembl_gene_id", attributes= c("ensembl_gene_id", "external_gene_name"), values=info$Probe, mart= mart)
  811. #join the dataframe with moduleinfo by Probe column
  812. joinGeneSymbol <- left_join(info, geneSymbols, by = c("Probe" = "ensembl_gene_id"))
  813. ```
  816. ```{r}
  817. #split function
  818. modulecolors_split <- split(joinGeneSymbol, joinGeneSymbol$ModuleColor)
  821. # Assuming you have a list of module colors and their corresponding data frames
  822. #module_colors <- c('red', 'blue', 'green', ...) # Add other module colors as needed
  823. #module_data_frames <- list(red = redmodule, blue = bluemodule, green = greenmodule, ...) # Replace with actual data frames for other modules
  825. # Create a function to perform enrichment analysis and save the results
  826. perform_enrichment <- function(module_color, module_data) {
  827. enriched_data <- enrichR::enrichr(
  828. module_data$external_gene_name,
  829. c('GO_Biological_Process_2023', 'GO_Cellular_Component_2023', 'GO_Molecular_Function_2023')
  830. )
  831. file_path <- sprintf("/Users/kchau/Documents/WGCNA/ConsensusBrainComparisonPCBonly/GO_terms/%s_excel.xlsx", module_color)
  832. openxlsx::write.xlsx(enriched_data, file_path)
  833. }
  835. # Loop through each module color and perform enrichment analysis
  836. for (color in names(modulecolors_split)) {
  837. module_data <- modulecolors_split[[color]]
  838. perform_enrichment(color, module_data)
  839. }
  840. ```
  843. ```{r}
  844. # Choose interesting modules
  845. intModules = c("MEbrown")
  846. for (module in intModules)
  847. {
  848. #Select module probes
  849. modGenes = (moduleColors==module)
  850. #Get their entrez ID codes
  851. modLLIDs = F_Brain1$entrez[modGenes];
  852. #Write them into a file
  853. fileName = paste("LocusLinkIDs-", module, ".txt", sep="");
  854. write.table(as.data.frame(modLLIDs), file = "/Users/kchau/Documents/WGCNA/ConsensusBrainComparison/outputgenelist",
  855. row.names = FALSE, col.names = FALSE)
  856. }
  857. ```
  859. ```{r}
  860. install.packages("onewaytests")
  861. library(onewaytests)
  862. bf.test()
  863. ```
  865. ```{r}
  866. MEtrait <- PCBtraitData[c(1:6)]
  868. femaleMEtrait <- MEtrait %>% filter(PCBtraitData$Sex == "F")
  869. femaleMEtrait <- femaleMEtrait[-c(32:51),] #remove samples with Folate
  871. femaleME_lightcyan_traits <- cbind(femaleMEtrait, MET[[1]]$data$MElightcyan)
  873. femaleME_lightcyan_traits$PCB_dose <- as.factor(femaleME_lightcyan_traits$PCB_dose)
  876. ggplotPCBEffect_lightcyan <- ggplot(data = femaleME_lightcyan_traits, aes(x = PCB_dose, y = MET[[1]]$data$MElightcyan, color = PCB_dose)) + geom_boxplot() + geom_point() + ggtitle("PCB Dosage Effect in Female Brain: Eigengene Values in Light Cyan Module") + xlab("PCB Dose (mg/kg/day)") + ylab("Module Eigengene Value")
  878. ggplotPCBEffect_lightcyan + theme(text = element_text(size = 10)) + geom_point(size=2) + ylim(-0.4, 0.4)
  879. ```
  881. ```{r}
  882. MEtrait <- PCBtraitData[c(1:6)]
  884. maleMEtrait <- MEtrait %>% filter(PCBtraitData$Sex == "M")
  885. maleMEtrait <- maleMEtrait[-c(28:45),] #remove samples with Folate
  887. maleME_lightcyan_traits <- cbind(maleMEtrait, MET[[2]]$data$MElightcyan)
  889. maleME_lightcyan_traits$PCB_dose <- as.factor(maleME_lightcyan_traits$PCB_dose)
  892. ggplotPCBEffect_lightcyan_male <- ggplot(data = maleME_lightcyan_traits, aes(x = PCB_dose, y = MET[[2]]$data$MElightcyan, color = PCB_dose)) + geom_boxplot() + geom_point() + ggtitle("PCB Dosage Effect in Male Brain: Eigengene Values in Light Cyan Module") + xlab("PCB Dose (mg/kg/day)") + ylab("Module Eigengene Value")
  894. ggplotPCBEffect_lightcyan_male + theme(text = element_text(size = 10)) + geom_point(size=2) + ylim(-0.4, 0.4)
  895. ```
  905. ```{r}
  906. ggplotTurquoise_male_new <- ggplot(data = MEturquoise_male, aes(x = FA_PCB, y = MET[[2]]$data$MEturquoise, color = FA_PCB)) + geom_boxplot() + geom_point() + ggtitle("Male Brain: Eigengene Values in Turquoise Module") + xlab("Treatment") + ylab("Module Eigengene Values")
  907. ```
  912. ```
  920. ```{r}
  921. #make the female Module Eigengene Dataframe
  922. MEtrait <- traitData[c(1:4)]
  925. MEtrait <- MEtrait[c(1:6, 14:23, 29:37, 44:49, 59:68, 77:86), ] #filter out the males
  926. MEtrait <- subset(MEtrait, MEtrait$Name != '36_01_BRAIN') #remove the outlier in female samples
  927. ```
  929. ```{r}
  930. MEturquoiseandtraits <- cbind(MEtrait, MET[[1]]$data$MEturquoise)
  931. ```
  933. ```{r}
  934. summary(MEturquoiseandtraits$`MET[[1]]$data$MEturquoise`)
  935. ```
  937. ```{r}
  938. MEturquoiseandtraits$highorlowdose <- c(1:50)
  939. MEturquoiseandtraits$highorlowdose[1:16] <- "low_dose"
  940. MEturquoiseandtraits$highorlowdose[17:31] <- "high_dose"
  941. MEturquoiseandtraits$highorlowdose[32:50] <- "low_dose"
  942. ```
  946. #want a boxplot separated by Yes Folic Acid (2 bars referring to 0 and 0.1 PCB) and No Folic Acid (2 bars reffering to 0 and 0.1 PCB dose)
  949. ```{r}
  950. MEturquoiseandtraits$FA_PCB <- paste(MEturquoiseandtraits$Folic_Acid, MEturquoiseandtraits$PCB_dose, sep = "-")
  951. ```
  953. ```{r}
  954. MEturquoiseandtraits$FA_PCB <- as.factor(MEturquoiseandtraits$FA_PCB)
  955. #levels(MEturquoiseandtraits$FA_PCB) <- c("No FA, No PCB", "No FA, Yes PCB", "Yes FA, No PBC", "Yes FA, Yes PCB")
  956. ```
  959. ```{r}
  960. ggplotTurquoise_female_new <- ggplot(data = MEturquoiseandtraits, aes(x = FA_PCB, y = MET[[1]]$data$MEturquoise, color = FA_PCB)) + geom_boxplot() + geom_point() + ggtitle("Female Brain: Eigengene Values in Turquoise Module") + xlab("Treatment") + ylab("Module Eigengene Values")
  961. ```
  963. ```{r}
  964. ggplotTurquoise_female_new
  965. ```
  967. 1) red and blue paired together, yellow and pink together (Folate Effect)
  968. subset
  971. ```{r}
  972. MEturquoise_female_folateeffect <- MEturquoiseandtraits[c(1:16, 32:50), ]
  974. MEturquoise_female_folateeffect$FA_PCB <- factor(MEturquoise_female_folateeffect$FA_PCB, levels=c("NO-0", "YES-0", "NO-0.1", "YES-0.1"))
  976. ggplotTurquoise_female_folateeffect <- ggplot(data = MEturquoise_female_folateeffect, aes(x = FA_PCB, y = MEturquoise_female_folateeffect$`MET[[1]]$data$MEturquoise`, color = FA_PCB)) + geom_boxplot() + geom_point() + ggtitle("Folate Effect in Female Brain: Eigengene Values in Turquoise Module") + xlab("Folic Acid - PCB Dose") + ylab("Module Eigengene Value")
  978. ggplotTurquoise_female_folateeffect
  983. ```
  990. 2) without folate (PCB dosage effect) #add color
  993. ```{r}
  994. MEfolateeffect <- MEturquoiseandtraits
  995. MEfolateeffect <- MEturquoiseandtraits[1:31, ] #only include samples with no folate
  996. MEfolateeffect$PCB_dose <- as.character(MEfolateeffect$PCB_dose)
  998. ggplotFolateEffect <- ggplot(data = MEfolateeffect, aes(x = MEfolateeffect$FA, y = MEfolateeffect$`MET[[1]]$data$MEturquoise`, color = MEfolateeffect$PCB_dose)) + geom_boxplot() + geom_point() + ggtitle("PCB Dosage Effect in Female Brain: Eigengene Values in Turquoise Module") + xlab("PCB Dose (mg/kg/day)") + ylab("Module Eigengene Value")
  1001. ggplotFolateEffect + theme(legend.position = "none")
  1002. ```
  1005. #Plot Male
  1007. ```{r}
  1008. MEtrait_male <- traitData[c(1:4)]
  1009. MEtrait_male <- MEtrait_male[c(7:13, 24:28, 38:43, 50:58, 69:76, 87:96), ] #filter out the females
  1010. #filter out male outliers
  1011. MEtrait_male <- subset(MEtrait_male, MEtrait_male$Name != '35_05_BRAIN') #remove the outlier in male samples
  1012. MEtrait_male <- subset(MEtrait_male, MEtrait_male$Name != '35_06_BRAIN') #remove the outlier in male samples
  1013. ```
  1016. ```{r}
  1017. MEturquoise_male <- cbind(MEtrait_male, MET[[2]]$data$MEturquoise)
  1018. ```
  1021. ```{r}
  1022. MEturquoise_male$FA_PCB <- paste(MEturquoise_male$Folic_Acid, MEturquoise_male$PCB_dose, sep = "-")
  1023. ```
  1025. ```{r}
  1026. MEturquoise_male$FA_PCB <- as.factor(MEturquoise_male$FA_PCB)
  1027. #levels(MEturquoiseandtraits$FA_PCB) <- c("No FA, No PCB", "No FA, Yes PCB", "Yes FA, No PBC", "Yes FA, Yes PCB")
  1028. ```
  1031. ```{r}
  1032. ggplotTurquoise_male_new <- ggplot(data = MEturquoise_male, aes(x = FA_PCB, y = MET[[2]]$data$MEturquoise, color = FA_PCB)) + geom_boxplot() + geom_point() + ggtitle("Male Brain: Eigengene Values in Turquoise Module") + xlab("Treatment") + ylab("Module Eigengene Values")
  1033. ```
  1035. ```{r}
  1036. ggplotTurquoise_male_new
  1037. ```
  1038. ```{r}
  1039. MEturquoise_male_folateeffect <- MEturquoise_male[c(1:12, 28:43), ]
  1042. MEturquoise_male_folateeffect$FA_PCB <- factor(MEturquoise_male_folateeffect$FA_PCB, levels=c("NO-0", "YES-0", "NO-0.1", "YES-0.1"))
  1044. ggplotTurquoise_male_folateeffect <- ggplot(data = MEturquoise_male_folateeffect, aes(x = FA_PCB, y = MEturquoise_male_folateeffect$`MET[[2]]$data$MEturquoise`, color = FA_PCB)) + geom_boxplot() + geom_point() + ggtitle("Folate Effect in Male Brain: Eigengene Values in Turquoise Module") + xlab("Folate - PCB Dose (mg/kg/day)") + ylab("Module Eigengene Value")
  1046. ggplotTurquoise_male_folateeffect
  1049. ```
  1059. 1) red and blue paired together, yellow and pink together (Folate Effect)
  1061. 2) without folate (PCB dosage effect) #add color
  1063. ```{r}
  1066. MEfolateeffect_male <- MEturquoise_male
  1067. MEfolateeffect_male <- MEturquoise_male[1:27, ] #only include samples with no folate
  1068. MEfolateeffect_male$FA_PCB <- MEfolateeffect_male$PCB_dose
  1070. MEfolateeffect_male$PCB_dose <- as.character(MEfolateeffect_male$PCB_dose)
  1072. ggplotFolateEffect_male <- ggplot(data = MEfolateeffect_male, aes(x = MEfolateeffect_male$PCB_dose, y = MEfolateeffect_male$`MET[[2]]$data$MEturquoise`, color = MEfolateeffect_male$PCB_dose)) + geom_boxplot() + geom_point() + ggtitle("PCB Dosage Effect in Male Brain: Eigengene Values in Turquoise Module") + xlab("PCB Dose (mg/kg/day)") + ylab("Module Eigengene Value")
  1075. ggplotFolateEffect_male + theme(legend.position = "none")
  1076. ```

ConsensusBrainComparisonPCBonly.Rmd at commit de8c052, under MIT · at the source

Overview

Authors: Kelly H. Chau1,2,3, Kari Neier1,2,3, Anthony E. Valenzuela4, Rebecca J. Schmidt3,5,6, Blythe Durbin-Johnson2,6, Pamela J. Lein3,4, Ian Korf2,7, Janine M. LaSalle1,2,3,5
ORCID iDs: Anthony E. Valenzuela, Blythe Durbin-Johnson
  1. Department of Medical Microbiology and Immunology, School of Medicine, University of California,Davis, Davis, CA 95616 USA
  2. UC Davis Genome Center, University of California,Davis, Davis, CA 95616 USA
  3. MIND Institute, School of Medicine, University of California,Davis, Sacramento, CA 95817 USA
  4. Department of Molecular Biosciences, School of Veterinary Medicine, University of California,Davis, Davis, CA 95616 USA
  5. Perinatal Origins of Disparities Center, University of California,Davis, Davis, CA 95616 USA
  6. Department of Public Health Sciences, School of Medicine, University of California,Davis, Davis, CA 95616 USA
  7. Department of Molecular and Cellular Biology, University of California,Davis, Davis, CA 95616 USA
Institutions: University of California, Davis (United States)
Journal: Genome biology, volume 27, issue 1, article 171
Dates: received 14 January 2025; accepted 20 March 2026; published online 7 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1186/s13059-026-04052-8 · PMID 41943130 · PMCID PMC13188613 · OpenAlex W7150709832
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), mouse (organism), other condition (population)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Connectivity, Spectral & time-frequency
Keywords: Placental-brain axis, Neurodevelopmental disorders, Polychlorinated biphenyls, Folic acid interactions, Weighted gene correlation network, Differential gene expression, X chromosome, Xist
MeSH: Brain*, Gene Regulatory Networks*, Placenta*, Polychlorinated Biphenyls*, Prenatal Exposure Delayed Effects*, Animals, Developmental Origins of Health and Disease, Female, Humans, Male, Mice, Pregnancy, Sex Characteristics (* major topic)
Topic: Toxic Organic Pollutants Impact (Health, Toxicology and Mutagenesis, Environmental Science), according to OpenAlex
Funding: U.S. Department of Health & Human Services | NIH | National Institute of Environmental Health Sciences (NIEHS) (R01ES029213); National Institute of Child Health and Human Development (P50HD103526)
Citations: cited by 1 paper (Europe PMC); 55 references in the paper

Abstract

Background: Neurodevelopmental disorders have a strong male bias that is poorly understood. The placenta provides molecular information about environmental interactions with genetics (including biological sex) that shape developmental processes in the brain. We investigate placental-brain transcriptional responses in an established mouse model of prenatal exposure to a human-relevant mixture of polychlorinated biphenyls (PCBs).

Results: To understand sex, tissue, and dosage effects in embryonic (E18) brain and placenta RNAseq data, we use weighted gene correlation network analysis (WGCNA) to create gene networks that could be compared across sex or tissue. WGCNA reveals that expression within most correlated gene networks is significantly and strongly associated with PCB exposure, but frequently in opposite directions between male–female and placenta-brain comparisons. In WGCNA and differentially expressed gene analyses, more transcriptional changes are observed in male brain than placenta, but the reverse is seen in females. Furthermore, female X-inactive specific transcript (Xist) levels correlate with sex-specific and non-monotonic PCB dose response, suggesting an X-linked protective epigenetic mechanism. The transcriptomic effects of low-dose PCB exposure are significantly opposed by dietary folic acid supplementation across both sexes but are strongest in female placentas. PCB and folic acid interacting gene networks are enriched in metabolic pathways involved in energy usage and translation, with female-specific protective effects enriched in PPAR, thermogenesis, glycerolipid, and O-glycan biosynthesis, as opposed to toxicant responses in male brain.

Conclusions: A female protective effect in response to prenatal PCB exposure appears to be mediated by dose-dependent sex differences in transcriptional modulation of placental metabolic pathways.

Supplementary Information: The online version contains supplementary material available at 10.1186/s13059-026-04052-8.

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

Repositories

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

kellychau/RNAseqWeightedGeneCorrelationNetworksPCBBrainPlacenta

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: de8c0521690de94a12039188cedf9739c1c91d6f, 13 March 2026
Languages: R (8)
Size: 14 files, 8 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, 8 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (8 files), WGCNA (8 files), data.table (7 files), ggplot2 (5 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
10 files
At the source: github.com/kellychau/RNAseqWeightedGeneCorrelationNetworksPCBBrainPlacenta

Zenodo 19011412

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (8 files), WGCNA (8 files), data.table (7 files), ggplot2 (5 files)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
10 files
At the source:

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 16 scripts, each with its path and the digest of its content;
  • 2 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 RNAseq data have been deposited in GEO accession GSE315769 [53].The code used to generate the analyses in this study is available on GitHub: (https://github.com/kellychau/RNAseqWeightedGeneCorrelationNetworksPCBBrainPlacenta) [54] and is released under the MIT license. A DOI‑minted archival copy is available at (https://doi.org/10.5281/zenodo.19011412) [55].

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

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

Cite

This paper

Chau, K. H., Neier, K., Valenzuela, A. E., Schmidt, R. J., Durbin-Johnson, B., Lein, P. J., Korf, I., & LaSalle, J. M. (2026). Sex and tissue resolved co-expression networks reveal a female placental-brain axis protective against prenatal PCB exposure. Genome biology, 27(1), 171. https://doi.org/10.1186/s13059-026-04052-8

BibTeX

@article{chau2026sex,
author = {Chau, Kelly H. and Neier, Kari and Valenzuela, Anthony E. and Schmidt, Rebecca J. and Durbin-Johnson, Blythe and Lein, Pamela J. and Korf, Ian and LaSalle, Janine M.},
title = {{Sex and tissue resolved co-expression networks reveal a female placental-brain axis protective against prenatal PCB exposure}},
journal = {Genome biology},
year = {2026},
month = apr,
volume = {27},
number = {1},
pages = {171},
publisher = {BMC},
issn = {1474-7596},
doi = {10.1186/s13059-026-04052-8},
url = {https://doi.org/10.1186/s13059-026-04052-8},
pmid = {41943130},
pmcid = {PMC13188613}
}

RIS

TY - JOUR
AU - Chau, Kelly H.
AU - Neier, Kari
AU - Valenzuela, Anthony E.
AU - Schmidt, Rebecca J.
AU - Durbin-Johnson, Blythe
AU - Lein, Pamela J.
AU - Korf, Ian
AU - LaSalle, Janine M.
TI - Sex and tissue resolved co-expression networks reveal a female placental-brain axis protective against prenatal PCB exposure
T2 - Genome biology
J2 - Genome Biol
PY - 2026
DA - 2026/04/07
VL - 27
IS - 1
SP - 171
SN - 1474-7596
PB - BMC
DO - 10.1186/s13059-026-04052-8
UR - https://doi.org/10.1186/s13059-026-04052-8
LA - en
ER -

CSL-JSON

{
"id": "10.1186/s13059-026-04052-8",
"type": "article-journal",
"title": "Sex and tissue resolved co-expression networks reveal a female placental-brain axis protective against prenatal PCB exposure",
"container-title": "Genome biology",
"author": [
{
"family": "Chau",
"given": "Kelly H."
},
{
"family": "Neier",
"given": "Kari"
},
{
"family": "Valenzuela",
"given": "Anthony E."
},
{
"family": "Schmidt",
"given": "Rebecca J."
},
{
"family": "Durbin-Johnson",
"given": "Blythe"
},
{
"family": "Lein",
"given": "Pamela J."
},
{
"family": "Korf",
"given": "Ian"
},
{
"family": "LaSalle",
"given": "Janine M."
}
],
"container-title-short": "Genome Biol",
"volume": "27",
"issue": "1",
"page": "171",
"DOI": "10.1186/s13059-026-04052-8",
"PMID": "41943130",
"PMCID": "PMC13188613",
"ISSN": "1474-7596",
"publisher": "BMC",
"URL": "https://doi.org/10.1186/s13059-026-04052-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
7
]
]
}
}

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.1016/j.isci.2026.115573 [code]
Female cortical cellular mosaicism underlies shared MeCP2 and PCB impacted gene pathways.
Journal: iScience
In common: WGCNA, data.table, ggplot2, 1 other tool, mouse, 4 references, author Anthony Valenzuela
[2] doi:10.3390/ijms27177801
Expression of the Human R163C-&lt;i&gt;RYR1&lt;/i&gt; Gain-of-Function Mutation Modified 2,2',3,5',6-Pentachlorobiphenyl (PCB 95) Developmental Neurotoxicity in Weanling Mice.
Journal: International journal of molecular sciences
In common: other condition, mouse, 4 references
[3] doi:10.1038/s41467-026-75723-0 [code]
Spatial transcriptomics reveals distinct cell type dynamics following opioid dependence in female mice with the common human μ-opioid receptor variant Oprm1 A118G.
Journal: Nature communications
In common: WGCNA, data.table, ggplot2, 1 other tool, other condition, mouse, 1 reference
[4] doi:10.1016/j.xhgg.2026.100652 [code]
CRISPR-engineered deletion of POGZ alters transcription factor binding at promoters of genes involved in synaptic signaling.
Journal: HGG advances
In common: WGCNA, data.table, ggplot2, 1 other tool, 2 references
[5] doi:10.1126/scitranslmed.adq4529 [code]
Modulating alternative splicing of &lt;i&gt;MECP2&lt;/i&gt; is a potential therapeutic strategy for Rett syndrome.
Journal: Science translational medicine
In common: WGCNA, data.table, ggplot2, 1 other tool, other condition, mouse, 1 reference
[6] doi:10.1038/s41467-026-72598-z [code]
Functional impact of genetic background on variable expressivity in neurodevelopmental disorders.
Journal: Nature communications
In common: WGCNA, data.table, ggplot2, 1 other tool, other condition, 1 reference
[7] doi:10.1038/s41380-026-03686-1 [code]
Early oligodendrocyte dysfunction signature in Alzheimer's disease: Insights from DNA methylomics and transcriptomics.
Journal: Molecular psychiatry
In common: WGCNA, data.table, ggplot2, 1 other tool, mouse, 1 reference
[8] doi:10.1038/s41467-026-74038-4 [code]
Semaglutide attenuates neuroinflammation in male mice.
Journal: Nature communications
In common: WGCNA, data.table, ggplot2, 1 other tool, mouse, 1 reference
[9] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: WGCNA, data.table, ggplot2, 1 other tool, mouse, 1 reference
[10] doi:10.1038/s41380-026-03585-5 [code]
Multiomics analysis identifies VPA-induced changes in neural progenitor cells, ventricular-like regions, and cellular microenvironment in dorsal forebrain organoids.
Journal: Molecular psychiatry
In common: WGCNA, data.table, ggplot2, 1 other tool, mouse, 1 reference

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.