OSCR

Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice identifies multiple high confidence candidate modifiers of Enteric Nervous System development.

Code ↔ Paper

20 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 20 matches
  1. [1] § Methods › Generation and processing of fetal ENS single nucleus assay for transposase-accessible chromatin-sequencing (snATAC-seq) data › Analysis of snATAC-seq data. ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 487–530 · score 0.93 · position frequency matrices, TFBSTools, getMatrixSet, TF binding motifs, accessible loci, JASPAR
  2. [2] § Results › Evolutionarily conserved SOX10 binding sites within Sox10Dom modifier intervals overlap or are near genes differentially expressed in the migrating wavefront ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 788–828 · score 0.93 · Cyp7b1, Fam204a, Ppargc1a, conserved SOX10, nearest gene, binding motif
  3. [3] § Methods › Modifier intervals for downstream candidate gene analysis ↔ 2024-10-04_ModifierInterval_CandidateGenePipeline_Part1_V3.R, lines 472–518 · score 0.89 · Browser Tool, BothSexNo15X5, BinBothSexAllChrom, BinFemale, BinMale, gene symbol
  4. [4] § Results › Evolutionarily conserved SOX10 binding sites within Sox10Dom modifier intervals overlap or are near genes differentially expressed in the migrating wavefront ↔ 2024-12-05_CandidateGenesComparisonsCrossModality.R, lines 175–205 · score 0.89 · Cyp7b1, Fam204a, Ppargc1a, nearest gene, Adgrb3, Bcas3
  5. [5] § Methods › Modifier intervals for downstream candidate gene analysis ↔ 2024-10-02_ModifierInterval_CandidateGenePipeline_Part2+3_V1.R, lines 247–292 · score 0.88 · Browser Tool, BothSexNo15X5, BinBothSexAllChrom, BinFemale, BinMale, gene symbol
  6. [6] § Results › Differential expression of modifier interval genes in the migrating wavefront of enteric neural crest-derived cells further prioritizes candidate genes ↔ 2024-10-03_ModifierInterval_CandidateGenePipeline_Part5_V1.R, lines 361–426 · score 0.87 · Hs3st3b1, Sh2b3, Slc10a4, Col1a1, Alkbh5, Hoxb5
  7. [7] § Methods › Generation and processing of fetal ENS single nucleus assay for transposase-accessible chromatin-sequencing (snATAC-seq) data › Analysis of snATAC-seq data. ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 245–287 · score 0.86 · log2 Fold Change, FindAllMarkers, Bonferroni adjusted, accessible locus, snATAC, identities
  8. [8] § Results › Prioritization of top aganglionosis modifier candidate genes across data modalities ↔ 2024-12-18_CandidateGenes_VariantsDifferentThanB6_FunctionalPrediction_V2.R, lines 295–357 · score 0.79 · UCSC Genome Browser, variant predicted, Col1a1, insertion, splice, Phox2b
  9. [9] § Results › Prioritization of top aganglionosis modifier candidate genes across data modalities ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 788–828 · score 0.78 · CellChat, TF binding motifs, conserved SOX10, Col1a1, expressed genes, Mecom
  10. [10] § Results › Chromatin accessibility within enteric neuronal progenitors highlights putative regulatory regions within aganglionosis modifier intervals ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 917–960 · score 0.78 · 33400347–33400362, 97891443–97891462, SOX10 binding sites, SOX10 binding motifs, chr11, Gopinath
  11. [11] § Results › Candidate modifiers of aganglionosis related to human GI phenotypes ↔ 2024-12-18_CandidateGenes_OverlapWithHumanGWASHit.R, lines 363–394 · score 0.78 · stool frequency GWAS, Mb away, HSCR GWAS, Col1a1, Rbpj, Sox10Dom
  12. [12] § Methods › Generation and processing of fetal ENS single nucleus assay for transposase-accessible chromatin-sequencing (snATAC-seq) data › Analysis of snATAC-seq data. ↔ 2024-09-16_snATAC-seq_Phox2bH2B-CFP_16.6dpc_ImportProcessing.R, lines 46–86 · score 0.72 · LSI dimensionality reduction, CellRanger, CFP, ATAC, Harmony, imported
  13. [13] § Methods › Overlap of Sox10Dom aganglionosis modifier interval candidate genes with human HSCR and stool frequency GWAS summary statistics ↔ 2024-12-18_CandidateGenes_OverlapWithHumanGWASHit.R, lines 363–394 · score 0.66 · stool frequency, HSCR GWAS, hg19, FDR, candidate genes, overlapped
  14. [14] § Results › Differential expression of modifier interval genes in the migrating wavefront of enteric neural crest-derived cells further prioritizes candidate genes ↔ 2024-10-03_ModifierInterval_CandidateGenePipeline_Part5_V1.R, lines 361–426 · score 0.62 · Phox2a, Col1a1, Ascl1, Prph, Tgfbi, Uchl1
  15. [15] § Methods › Generation and processing of fetal ENS single nucleus assay for transposase-accessible chromatin-sequencing (snATAC-seq) data › Isolation of nuclei from fetal mouse tissue. ↔ 2024-09-16_snATAC-seq_Phox2bH2B-CFP_16.6dpc_ImportProcessing.R, lines 1–44 · score 0.61 · ATAC seq, neuronal cells, snATAC, Libraries, CFP, Nuclei
  16. [16] § Results › Chromatin accessibility within enteric neuronal progenitors highlights putative regulatory regions within aganglionosis modifier intervals ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 54–95 · score 0.60 · Phox2b H2B CFP, scRNA, snATAC, activity, WT, nucleus
  17. [17] § Results › Prioritization of candidate genes based on expression in the fetal mouse intestine ↔ 2024-10-02_ModifierInterval_CandidateGenePipeline_Part2+3_V1.R, lines 344–382 · score 0.59 · 9.5–15.5, fetal gut, known genes, modifier intervals, dpc, candidate genes
  18. [18] § Results › Differential expression of modifier interval genes in the migrating wavefront of enteric neural crest-derived cells further prioritizes candidate genes ↔ 2024-10-03_ModifierInterval_CandidateGenePipeline_Part5_V1.R, lines 229–273 · score 0.57 · Phox2a, enteric neural crest, Volcano, gene expression, wavefront, modifier interval
  19. [19] § Results › Differential expression of modifier interval genes in the migrating wavefront of enteric neural crest-derived cells further prioritizes candidate genes ↔ 2024-12-26_CandidateGene_ExonPheWAS_ManhattanPlots.R, lines 114–171 · score 0.53 · COL1A1, Hoxb5, Myh10, Prph, Tgfbi, Uchl1
  20. [20] § Results › Chromatin accessibility within enteric neuronal progenitors highlights putative regulatory regions within aganglionosis modifier intervals ↔ 2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, lines 575–627 · score 0.53 · percent.diff, TF binding motifs, background, nuclei, chromatin, neuroblast

Paper

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

The paper is loaded when this pane is shown.

The authors' code

R · 1,092 lines · 55 KB · CC-BY-4.0 · 7 matches

  1. ################################################################################
  2. ### Modifier Interval Candidate Gene Pipeline Part 7:
  3. ### How many differentially accessible loci in WT 16.5dpc ENS nuclei (by
  4. ### cell group: Neuronal Branches, Neuroblast, Progenitor) are within modifier
  5. ### intervals? Which TF binding motifs are enriched within modifier interval-
  6. ### contained differentially accessible loci?
  7. ################################################################################
  9. ### Import packages:
  11. library(Seurat)
  12. library(Signac)
  13. library(cowplot)
  14. library(GenomicRanges)
  15. library(TFBSTools)
  16. library(JASPAR2024)
  17. library(BSgenome.Mmusculus.UCSC.mm10)
  18. library(tidyverse)
  19. library(TxDb.Mmusculus.UCSC.mm10.knownGene)
  20. library(ggVennDiagram)
  22. ### Get filepaths where data are stored:
  24. savepath='path/to/files/'
  25. savepath='C:/Users/jtben/Documents/ms2/MS2_Files/Manuscripts/Sox10DomPedigree_Dach1/'
  27. path10='D:/Projects/Gopinath_2016_SOX10_BindingConsensus_HMC/'
  29. path5='G:/Documents/ms2/MS2_Files/Projects/snATAC-seq/Dim1Included/'
  31. pathtoobjects<-'D:/JustinPaperCodeReview/JustinObjects/'
  33. path4='G:/Documents/ms2/MS2_Files/Projects/snATAC-seq/'
  35. path8='D:/Projects/Zhao_et_al_2022_scRNA-seq_GITractDev9.5-15.5dpc/'
  37. ### Import GWAS Summary Statistics and Modifier Interval Set Tables:
  39. GWASList<-readRDS(file = paste0(savepath,
  40. '2024-08-21_AllGWASCompiled.RDS'))
  41. IntList<-readRDS(file = paste0(savepath,
  42. '2024-08-21_AllGWAS-DerivedModifierIntervalsSetsCompiled.RDS'))
  43. IntGeneExpData<-readRDS(file = paste0(savepath,
  44. '2024-08-21_AllGWAS-DerivedModifierIntervalsSetsCompiled_GenesWithinData.RDS'))
  45. IntGeneExpList<-readRDS(file = paste0(savepath,
  46. '2024-08-21_AllGWAS-DerivedModifierIntervalsSetsCompiled_GenesWithinUnique.RDS'))
  47. TopSNPList<-readRDS(file = paste0(savepath,
  48. '2024-08-21_AllGWAS-DerivedModifierIntervalsSets_TopSNPs.RDS'))
  49. PeaksPerInterval<-readRDS(file = paste0(savepath,
  50. '2024-08-29_AllGWAS-DerivedModifierIntervalsSets_PeakSNPsPerModInt.RDS'))
  51. GroupsSigList<-readRDS(file = paste0(savepath,
  52. '2024-10-04_AllGWASCompiled_WithModInts.RDS'))
  54. ### Now I will load in the E16.5 WT snATAC-seq data
  55. ### generated from Phox2b H2B-CFP:
  56. reprocessed<-readRDS(file = paste0(path5,'2023-08-07_reprocessed_No1718+GeneralLabels+FragUpdate+Links.RDS'))
  57. DimPlot(reprocessed,label=TRUE)
  58. DefaultAssay(reprocessed)<-'MACS2peaks'
  59. Idents(reprocessed)<-'seurat_clusters'
  61. ### Load in data from Avila & Benthal et al., 2024:
  63. JAAE15.5 <- readRDS(paste0(pathtoobjects,
  64. "2024-08-23_ENS_MainClusts_JAA.RDS"))
  66. ### Make sure the grouped identities are present in the object:
  68. DimPlot(JAAE15.5,label=TRUE,group.by = 'Grouped_Ident')
  70. ### Make sure the Sub_cluster identities are present in the object:
  72. JAAE15.5<-SetIdent(JAAE15.5,value = 'sub_cluster')
  73. DimPlot(JAAE15.5,label=TRUE)
  75. ### Isolate only the WT cells from the object:
  77. WT<-subset(JAAE15.5,subset = Condition == c("Wt"))
  79. ### Make the default assay SCT:
  81. DefaultAssay(WT)<-'SCT'
  83. ### Below I use the estimated RNA expression in the snATAC-seq data to match
  84. ### cell identities to the scRNA-seq data:
  86. ### Note that Neuroblast Path 2 cells may not exist in the snATAC-seq data
  87. ### because these are nuclei, and Path 2 cells are cycling. M phase cells may
  88. ### have condensed chromosomes and therefore no nuclear membrane, making nuclei
  89. ### isolation from these cells impossible.
  91. ### quantify gene activity
  92. ### note that we restrict the imputation to variable genes from scRNA-seq,
  93. ### but could impute the full transcriptome if we wanted to
  94. genes.use <- VariableFeatures(WT)
  95. refdata <- GetAssayData(WT, assay = "RNA", slot = "data")[genes.use, ]
  97. ### Identify anchors
  98. transfer.anchors <- FindTransferAnchors(reference = WT,
  99. query = reprocessed,
  100. features = VariableFeatures(object = WT),
  101. reference.assay = "RNA",
  102. query.assay = "RNA", reduction = "cca")
  104. ### Transfer data to estimate which cell types correspond to snATAC-seq nuclei:
  106. celltype.predictions <- TransferData(anchorset = transfer.anchors,
  107. refdata = WT$Grouped_Ident,
  108. weight.reduction = reprocessed[["lsi"]],
  109. dims = 2:30)
  111. ### Add the metadata for the cell type predictions to the snATAC-seq data:
  113. reprocessed <- AddMetaData(reprocessed, metadata = celltype.predictions)
  115. ### Add the predicted IDs in the correct order to the snATAC-seq metadata:
  117. reprocessed$predicted.id <- factor(reprocessed$predicted.id,
  118. levels = c('Progenitors','Path1',
  119. 'Path2',"Branch.A",
  120. "Branch.B","Branch.C"))
  122. ### Make the predicted cell IDs the main cluster identities in the snATAC-seq
  123. ### data and make the grouped cell IDs the main cluster identities in the
  124. ### scRNA-seq reference data:
  126. Idents(reprocessed)<-'predicted.id'
  127. Idents(WT)<-'Grouped_Ident'
  129. ### Plot the UMAPs beside each other:
  131. p1 <- DimPlot(reprocessed)
  132. p2 <- DimPlot(WT)
  133. p1 | p2
  135. ### Set up dataframe to calculate confidence of cell type predictions:
  137. predictions <- table(reprocessed$seurat_clusters, reprocessed$predicted.id)
  138. predictions <- predictions/rowSums(predictions) # normalize for number of cells in each cell type
  139. predictions <- as.data.frame(predictions)
  140. predictions<-predictions[predictions$Var1 %in% c('0','1','2','3','4','5','6',
  141. '7','8','9','10','11','12',
  142. '13','14','15','16'),]
  144. ### Plot a heatmap of the cell type prediction confidence:
  146. p3 <- ggplot(predictions, aes(Var1, Var2, fill = Freq)) +
  147. geom_tile() + scale_fill_gradient(name = "Fraction \nof cells",
  148. low = "#ffffc8", high = "#7d0025") +
  149. xlab("Seurat Clusters") +
  150. ylab("Predicted cell type label from RNA") +
  151. theme_cowplot() +
  152. theme(axis.text.x = element_text(angle = 90,
  153. vjust = 0.5,
  154. hjust = 1))
  155. p1 | p2 | p3
  157. ### Save these plots:
  159. ggsave(plot = p1+
  160. theme(axis.line = element_line(colour = 'black', size = 1.1),
  161. axis.text.y = element_text(size = 10,face='bold'),
  162. legend.text=element_text(size=15,face='bold'),
  163. plot.title = element_text(size=20,face='bold'),
  164. legend.title = element_text(size=15,face='bold'),
  165. axis.text.x = element_text(size = 10,face='bold')),
  166. filename = paste0(savepath,"Images/2024-09-09_Reprocessed_snATAC-seqPhox2bBright2+Dim1+2_predicted.id.jpg"),
  167. device = 'jpg',dpi = 600,units = 'in',width = 5,height = 3)
  168. ggsave(plot = p2+
  169. theme(axis.line = element_line(colour = 'black', size = 1.1),
  170. axis.text.y = element_text(size = 10,face='bold'),
  171. legend.text=element_text(size=15,face='bold'),
  172. plot.title = element_text(size=20,face='bold'),
  173. legend.title = element_text(size=15,face='bold'),
  174. axis.text.x = element_text(size = 10,face='bold')),
  175. filename = paste0(savepath,"Images/2024-09-09_JAA15.5scRNA-seq_Grouped_Ident.jpg"),
  176. device = 'jpg',dpi = 600,units = 'in',width = 5,height = 3)
  177. ggsave(plot = p3+
  178. theme(axis.line = element_line(colour = 'black', size = 1.1),
  179. axis.text.y = element_text(size = 12,face='bold'),
  180. legend.text=element_text(size=15,face='bold'),
  181. plot.title = element_text(size=20,face='bold'),
  182. legend.title = element_text(size=15,face='bold'),
  183. axis.text.x = element_text(size = 12,face='bold'))+RotatedAxis(),
  184. filename = paste0(savepath,"Images/2024-09-09_Predicted.ID_Heatmap.jpg"),
  185. device = 'jpg',dpi = 600,units = 'in',width = 7,height = 5)
  187. ### Read in the Bright and Dim Phox2b-CFP snATAC-seq datasets from which the
  188. ### main snATAC-seq datasets originated for plotting purposes:
  190. Bright2<-readRDS(file = paste0(path4,'2021-08-26_snATAC-seq_Phox2bCFP_Bright2.RDS'))
  191. hm.Dim.integrate<-readRDS(file = paste0(path4,'2021-08-20_snATAC-seq_Phox2bCFP_Dims_IntegratedBatchCorrectedHarmony.RDS'))
  192. DimPlot(hm.Dim.integrate,label=TRUE) | DimPlot(Bright2,label=TRUE)
  193. gc()
  195. ### Plot UMAPs and save the plots:
  197. ggsave(plot = DimPlot(reprocessed)+
  198. theme(axis.line = element_line(colour = 'black', size = 1.1),
  199. axis.text.y = element_text(size = 10,face='bold'),
  200. legend.text=element_text(size=15,face='bold'),
  201. plot.title = element_text(size=20,face='bold'),
  202. legend.title = element_text(size=15,face='bold'),
  203. axis.text.x = element_text(size = 10,face='bold')),
  204. filename = paste0(savepath,"Images/2024-09-08_Reprocessed_snATAC-seqPhox2bBright2+Dim1+2.jpg"),
  205. device = 'jpg',dpi = 600,units = 'in',width = 5,height = 3)
  206. ggsave(plot = DimPlot(reprocessed)+
  207. theme(axis.line = element_line(colour = 'black', size = 1.1),
  208. axis.text.y = element_text(size = 10,face='bold'),
  209. legend.text=element_text(size=15,face='bold'),
  210. plot.title = element_text(size=20,face='bold'),
  211. legend.title = element_text(size=15,face='bold'),
  212. axis.text.x = element_text(size = 10,face='bold')),
  213. filename = paste0(savepath,"Images/2024-09-08_Reprocessed_snATAC-seqPhox2bBright2+Dim1+2_forLegend.jpg"),
  214. device = 'jpg',dpi = 600,units = 'in',width = 7,height = 5)
  216. ggsave(plot = DimPlot(hm.Dim.integrate)+
  217. theme(axis.line = element_line(colour = 'black', size = 1.1),
  218. axis.text.y = element_text(size = 10,face='bold'),
  219. legend.text=element_text(size=15,face='bold'),
  220. plot.title = element_text(size=20,face='bold'),
  221. legend.title = element_text(size=15,face='bold'),
  222. axis.text.x = element_text(size = 10,face='bold')),
  223. filename = paste0(savepath,"Images/2024-09-08_Reprocessed_snATAC-seqPhox2bDim1+2_HBC.jpg"),
  224. device = 'jpg',dpi = 600,units = 'in',width = 5,height = 3)
  226. ggsave(plot = DimPlot(Bright2)+
  227. theme(axis.line = element_line(colour = 'black', size = 1.1),
  228. axis.text.y = element_text(size = 10,face='bold'),
  229. legend.text=element_text(size=15,face='bold'),
  230. plot.title = element_text(size=20,face='bold'),
  231. legend.title = element_text(size=15,face='bold'),
  232. axis.text.x = element_text(size = 10,face='bold')),
  233. filename = paste0(savepath,"Images/2024-09-08_Reprocessed_snATAC-seqPhox2bBright2.jpg"),
  234. device = 'jpg',dpi = 600,units = 'in',width = 5,height = 3)
  235. ggsave(plot = DimPlot(Bright2)+
  236. theme(axis.line = element_line(colour = 'black', size = 1.1),
  237. axis.text.y = element_text(size = 10,face='bold'),
  238. legend.text=element_text(size=15,face='bold'),
  239. plot.title = element_text(size=20,face='bold'),
  240. legend.title = element_text(size=15,face='bold'),
  241. axis.text.x = element_text(size = 10,face='bold')),
  242. filename = paste0(savepath,"Images/2024-09-08_Reprocessed_snATAC-seqPhox2bBright2_forLegend.jpg"),
  243. device = 'jpg',dpi = 600,units = 'in',width = 7,height = 5)
  245. ### Make certain the main cluster identities in the snATAC-seq data are the
  246. ### predicted cell IDs:
  248. Idents(reprocessed)<-'predicted.id'
  250. ### Perform FindAllMarkers to find the differentially accessible peaks per cell
  251. ### ID in the snATAC-seq UMAP:
  253. FAM<-FindAllMarkers(object = reprocessed,
  254. assay = 'MACS2peaks',
  255. return.thresh = 0.05)
  257. ### Filter the differentially accessible peaks to those greater than 1.5
  258. ### average log2FoldChange and less than -1.5 average log2FoldChange:
  260. FAMpos<-FAM[FAM$avg_log2FC>1.5,]
  261. FAMneg<-FAM[FAM$avg_log2FC<((-1)*1.5),]
  262. FAM.sig<-rbind(FAMpos,FAMneg)
  264. ### Filter for Bonferroni-adjusted p-value less than 0.05:
  266. FAM.sig<-FAM.sig[FAM.sig$p_val_adj<0.05,]
  268. ### Change the name of a column to "loc" for loci:
  270. names(FAM.sig)[7]<-'loc'
  272. ### Split up the loci column to seqnames, start, and end columns per
  273. ### differentially accessible locus:
  275. FAM.sig[c('seqnames', 'start','end')] <- stringr::str_split_fixed(FAM.sig$loc, '-', 3)
  277. ### Flank each differentially accessible locus +/-25 base pairs:
  279. FAM.sig$start<-as.numeric(FAM.sig$start)-25
  280. FAM.sig$end<-as.numeric(FAM.sig$end)+25
  282. ### Save the differentially accessible loci dataset to a csv file:
  284. write.csv(x = FAM.sig,
  285. file = paste0(savepath,
  286. '2024-09-09_Phox2bCFPsnATAC-seqWT_DAPeaks_predicted.id.csv'),
  287. quote = FALSE)
  290. ### Find all Differentially Accessible Peaks conserved binding motifs that are
  291. ### within modifier intervals:
  293. ### Set up empty named list to write to:
  295. DAModIntList<-list(Male=NULL,
  296. Female=NULL,
  297. BothSexNo15X5=NULL,
  298. BothSexAllChrom=NULL,
  299. BinMale=NULL,
  300. BinFemale=NULL,
  301. BinBothSexAllChrom=NULL)
  303. ### Set up for loop to loop through the gene lists for each interval set to
  304. ### find which overlap with Differentially Accessible Peaks:
  306. for(i in seq_along(DAModIntList)){
  308. ### Initiate an empty list for writing Differentially Accessible Peaks to:
  310. DAModIntListmodlocilist<-list()
  311. for(j in seq_along(paste0('chr',IntList[[i]]$chrom))){
  312. DAModIntListmodlocilist[[j]]<-FAM.sig[FAM.sig$seqnames==paste0('chr',IntList[[i]]$chrom)[j] &
  313. FAM.sig$start>=IntList[[i]]$start0.5Mb[j] &
  314. FAM.sig$end<=IntList[[i]]$end0.5Mb[j],]
  317. try(DAModIntListmodlocilist[[j]]$ModInt<-IntList[[i]]$IntervalID[j])
  318. try(DAModIntListmodlocilist[[j]]$source<-names(IntList)[[i]])
  319. }
  321. DAModIntList[[i]]<-do.call('rbind',DAModIntListmodlocilist)
  322. }
  324. DAModIntList
  326. ### Add a name column to the data:
  328. for(i in seq_along(DAModIntList)){
  329. DAModIntList[[i]]$assay<-names(DAModIntList)[[i]]
  330. }
  332. ### Combine to one data table:
  334. FAM.sig.all<-do.call(rbind,DAModIntList)
  336. write.csv(x = FAM.sig.all,
  337. file = paste0(savepath,'2024-10-09_DifferentiallyAccessiblePeaksPredictedCellIDs_InModIntervals.csv'),
  338. quote = FALSE,row.names = FALSE)
  340. FAM.sig.all$IntervalID2<-paste0(FAM.sig.all$assay,'_',FAM.sig.all$ModInt)
  341. FAM.sig.all$IntervalID2Clust<-paste0(FAM.sig.all$cluster,'_',
  342. FAM.sig.all$IntervalID2)
  344. table(FAM.sig.all$IntervalID2,FAM.sig.all$cluster)
  346. ### Find overlap of the differentially accessible peaks across the modifier
  347. ### interval definitions:
  349. geneSets <- list(BothSexAllChrom = unique(FAM.sig.all[FAM.sig.all$assay=='BothSexAllChrom',]$loc),
  350. BothSexNo15X5 = unique(FAM.sig.all[FAM.sig.all$assay=='BothSexNo15X5',]$loc),
  351. Female = unique(FAM.sig.all[FAM.sig.all$assay=='Female',]$loc),
  352. Male = unique(FAM.sig.all[FAM.sig.all$assay=='Male',]$loc),
  353. BinBothSexAllChrom = unique(FAM.sig.all[FAM.sig.all$assay=='BinBothSexAllChrom',]$loc),
  354. BinFemale = unique(FAM.sig.all[FAM.sig.all$assay=='BinFemale',]$loc),
  355. BinMale = unique(FAM.sig.all[FAM.sig.all$assay=='BinMale',]$loc))
  357. vnnplt<-ggVennDiagram(geneSets,
  358. category.names = c("All Chrs",
  359. "No Chr5,15,X",
  360. "Female",
  361. "Male",
  362. "Bin All Chrs",
  363. 'BinFemale',
  364. 'BinMale'),
  365. # label_size = 5,
  366. # set_size = 5,
  367. force_upset = TRUE,
  368. label_size = 20,top.bar.numbers.size = 6,
  369. set_size = 20
  370. )
  371. # +scale_fill_distiller(palette = "Purples",direction=1)
  372. vnnplt
  374. vnnplt$plotlist[[1]][["theme"]][["text"]]$face<-'bold'
  375. vnnplt$plotlist[[1]][["theme"]][["text"]]$size<-17
  376. vnnplt$plotlist[[2]][["theme"]][["text"]]$face<-'bold'
  377. vnnplt$plotlist[[2]][["theme"]][["text"]]$size<-17
  378. vnnplt$plotlist[[3]][["theme"]][["text"]]$face<-'bold'
  379. vnnplt$plotlist[[3]][["theme"]][["text"]]$size<-17
  381. vnnplt
  383. ggsave(device = 'jpg',dpi = 600,
  384. plot = vnnplt,width = 13,height = 5,
  385. filename = paste0(savepath,"Images/2024-10-09_DiffAccChrom_snATAC-seq_Phox2b_predicted.id_UpsetPlot.jpg"))
  387. ### Convert differentially accessible loci to GRanges objects so that
  388. ### plotting ChIPseeker annotations can be performed per group:
  390. ### DA Peaks overall:
  392. FAM.sig_GR<-GenomicRanges::makeGRangesFromDataFrame(df = FAM.sig,
  393. keep.extra.columns = TRUE,
  394. seqnames.field = 'seqnames',
  395. start.field = 'start',
  396. end.field = 'end',
  397. ignore.strand = TRUE)
  399. ### Annotate the loci:
  401. bed.annot=ChIPseeker::annotatePeak(FAM.sig_GR,
  402. tssRegion=c(-3000, 3000),
  403. TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
  404. annoDb="org.Mm.eg.db")
  406. annot_peaks=as.data.frame(bed.annot)
  408. ### Annotate the peaks per predicted cell ID:
  410. FAM.sig_GR.list<- split(FAM.sig_GR , f = FAM.sig_GR$cluster)
  411. peakAnnoList <- lapply(FAM.sig_GR.list, ChIPseeker::annotatePeak,
  412. TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
  413. tssRegion=c(-3000, 3000), verbose=FALSE,
  414. annoDb="org.Mm.eg.db")
  416. ### Plot the annotations in stacked bar plots:
  418. ChIPseeker::plotAnnoBar(peakAnnoList)
  419. ggsave(plot = ChIPseeker::plotAnnoBar(peakAnnoList)+
  420. theme(axis.line = element_line(colour = 'black', size = 1.1),
  421. axis.text.y = element_text(size = 10,face='bold'),
  422. legend.text=element_text(size=15,face='bold'),
  423. plot.title = element_text(size=20,face='bold'),
  424. legend.title = element_text(size=15,face='bold'),
  425. axis.text.x = element_text(size = 10,face='bold')),
  426. device = 'jpg',dpi = 600,width = 8,height = 2.5,
  427. filename = paste0(savepath,"Images/2024-09-09_ALLFAMPeaks_DiffAccChrom_snATAC-seq_Phox2b_predicted.id_PeakClass_splitbyCluster.jpg"))
  429. ### Now for the DA peaks in modifier intervals:
  431. FAM.sig.all_GR<-GenomicRanges::makeGRangesFromDataFrame(df = FAM.sig.all,
  432. keep.extra.columns = TRUE,
  433. seqnames.field = 'seqnames',
  434. start.field = 'start',
  435. end.field = 'end',
  436. ignore.strand = TRUE)
  438. ### Annotate the loci:
  440. bed.annot=ChIPseeker::annotatePeak(FAM.sig.all_GR,
  441. tssRegion=c(-3000, 3000),
  442. TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
  443. annoDb="org.Mm.eg.db")
  444. annot_peaks=as.data.frame(bed.annot)
  446. # ChIPseeker::upsetplot(bed.annot, vennpie=TRUE)
  447. # ChIPseeker::vennpie(x = bed.annot)
  448. # ChIPseeker::plotAnnoBar(bed.annot)
  450. ### Then Split by modifier definitions, or "assay" and plot:
  452. FAM.sig.all_GR.list<- split(FAM.sig.all_GR , f = FAM.sig.all_GR$assay)
  453. peakAnnoList <- lapply(FAM.sig.all_GR.list, ChIPseeker::annotatePeak,
  454. TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
  455. tssRegion=c(-3000, 3000), verbose=FALSE,
  456. annoDb="org.Mm.eg.db")
  457. ChIPseeker::plotAnnoBar(peakAnnoList)
  458. ggsave(plot = ChIPseeker::plotAnnoBar(peakAnnoList)+
  459. theme(axis.line = element_line(colour = 'black', size = 1.1),
  460. axis.text.y = element_text(size = 10,face='bold'),
  461. legend.text=element_text(size=15,face='bold'),
  462. plot.title = element_text(size=20,face='bold'),
  463. legend.title = element_text(size=15,face='bold'),
  464. axis.text.x = element_text(size = 10,face='bold')),
  465. device = 'jpg',dpi = 600,width = 8,height = 2,
  466. filename = paste0(savepath,"Images/2024-09-09_DiffAccChrom_snATAC-seq_Phox2b_predicted.id_PeakClass_splitbyAssay.jpg"))
  468. ### Then split by cluster definitions and plot:
  470. FAM.sig.all_GR.list<- split(FAM.sig.all_GR , f = FAM.sig.all_GR$cluster)
  471. peakAnnoList <- lapply(FAM.sig.all_GR.list, ChIPseeker::annotatePeak,
  472. TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
  473. tssRegion=c(-3000, 3000), verbose=FALSE,
  474. annoDb="org.Mm.eg.db")
  475. ChIPseeker::plotAnnoBar(peakAnnoList)
  476. ggsave(plot = ChIPseeker::plotAnnoBar(peakAnnoList)+
  477. theme(axis.line = element_line(colour = 'black', size = 1.1),
  478. axis.text.y = element_text(size = 10,face='bold'),
  479. legend.text=element_text(size=15,face='bold'),
  480. plot.title = element_text(size=20,face='bold'),
  481. legend.title = element_text(size=15,face='bold'),
  482. axis.text.x = element_text(size = 10,face='bold')),
  483. device = 'jpg',dpi = 600,width = 8,height = 2.5,
  484. filename = paste0(savepath,"Images/2024-09-09_DiffAccChrom_snATAC-seq_Phox2b_predicted.id_PeakClass_splitbyCluster.jpg"))
  487. ### Then a combination of both to see how these might differ across cluster and
  488. ### modifier interval definitions:
  490. ### Make combined metadata of cluster and GEMMA run:
  492. FAM.sig.all_GR$assay_cluster<-paste0(FAM.sig.all_GR$assay,'_',
  493. FAM.sig.all_GR$cluster)
  495. ### Add annotations per assay_cluster:
  497. FAM.sig.all_GR.list2<- split(FAM.sig.all_GR , f = FAM.sig.all_GR$assay_cluster)
  498. peakAnnoList <- lapply(FAM.sig.all_GR.list2, ChIPseeker::annotatePeak,
  499. TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
  500. tssRegion=c(-3000, 3000), verbose=FALSE,
  501. annoDb="org.Mm.eg.db")
  503. ### Plot these:
  505. ChIPseeker::plotAnnoBar(peakAnnoList)
  506. ggsave(plot = ChIPseeker::plotAnnoBar(peakAnnoList)+
  507. theme(axis.line = element_line(colour = 'black', size = 1.1),
  508. axis.text.y = element_text(size = 10,face='bold'),
  509. legend.text=element_text(size=15,face='bold'),
  510. plot.title = element_text(size=20,face='bold'),
  511. legend.title = element_text(size=15,face='bold'),
  512. axis.text.x = element_text(size = 10,face='bold')),
  513. device = 'jpg',dpi = 600,width = 8,height = 7,
  514. filename = paste0(savepath,"Images/2024-09-09_DiffAccChrom_snATAC-seq_Phox2b_predicted.id_PeakClass_splitbyAssayandCluster.jpg"))
  516. ### Perform TF binding motif enrichment on the differentially accessible loci
  517. ### within/per modifier intervals via Signac's base function and JASPAR 2024:
  519. ### Get the position frequency matrices from JASPAR 2024 for vertebrates:
  521. pfm <- TFBSTools::getMatrixSet(
  522. x = db(JASPAR2024()),
  523. opts = list(collection = "CORE",
  524. tax_group = 'vertebrates',
  525. all_versions = FALSE))
  526. # pfm.un <- TFBSTools::getMatrixSet(
  527. # x = db(JASPAR2024()),
  528. # opts = list(collection = "UNVALIDATED",
  529. # tax_group = 'vertebrates',
  530. # all_versions = TRUE))
  532. ### Add the UNVALIDATED DACH1 TF binding motif information from JASPAR2020 to
  533. ### the pfm object:
  534. ### https://jaspar2020.genereg.net/matrix/UN0306.1/
  536. pfm@listData[['UN0306.1']]<-pfm@listData[['MA0442.3']]
  538. as.matrix(read.table(file = paste0(savepath,'UN0306.1.jaspar.tsv'),
  539. header = FALSE,sep = '\t',row.names = c('A','C','G','T')))
  541. pfm@listData[['UN0306.1']]@profileMatrix<-as.matrix(read.table(file = paste0(savepath,'UN0306.1.jaspar.tsv'),
  542. header = FALSE,
  543. sep = '\t',
  544. row.names = c('A','C','G','T')))
  545. pfm@listData[['UN0306.1']]@profileMatrix
  547. pfm@listData[['UN0306.1']]@ID<-'UN0306.1'
  548. pfm@listData[['UN0306.1']]@name<-'DACH1'
  549. pfm@listData[['UN0306.1']]@matrixClass<-'SKI / DACH domain factors'
  550. pfm@listData[['UN0306.1']]@tags$centrality_logp<-(-633.186)
  551. pfm@listData[['UN0306.1']]@tags$family<-'DACH / dachshund-related factors'
  552. pfm@listData[['UN0306.1']]@tags$medline<-NULL
  553. pfm@listData[['UN0306.1']]@tags$remap_tf_name<-'DACH1'
  554. pfm@listData[['UN0306.1']]@tags$source<-'http://remap.univ-amu.fr/'
  555. pfm@listData[['UN0306.1']]@tags$tax_group<-'vertebrates'
  556. pfm@listData[['UN0306.1']]@tags$tfbs_shape_id<-NULL
  557. pfm@listData[['UN0306.1']]@tags$unibind<-NULL
  558. pfm@listData[['UN0306.1']]@tags$collection<-'UNVALIDATED'
  559. pfm@listData[['UN0306.1']]@tags$acc<-'Q9UI36'
  561. ### Add motif information to the snATAC-seq object:
  563. reprocessed <- AddMotifs(
  564. object = reprocessed,
  565. genome = BSgenome.Mmusculus.UCSC.mm10,
  566. pfm = pfm)
  568. ### Find background open chromatin to assess differentially accessible
  569. ### loci within modifier loci against for motif analysis:
  571. ### Get all "open" loci from the snATAC-seq object:
  573. open.peaks <- AccessiblePeaks(reprocessed)
  575. ### Match the overall GC content in the peak set:
  577. meta.feature <- GetAssayData(reprocessed,
  578. assay = "MACS2peaks",
  579. layer = "meta.features")
  580. peaks.matched <- MatchRegionStats(
  581. meta.feature = meta.feature[open.peaks, ],
  582. query.feature = meta.feature[unique(FAM.sig$loc), ],
  583. n = 50000)
  585. ### Perform motif enrichment analysis on just the progenitor and neuroblast
  586. ### nuclei by a combined metadata of GEMMA GWAS run + Modifier Interval:
  588. ### Make a GRanges object for these differentially accessible peaks:
  590. FAM.sig.all_GR_Prog<-GenomicRanges::makeGRangesFromDataFrame(df = FAM.sig.all[FAM.sig.all$cluster %in%
  591. c('Progenitors',
  592. 'Neuroblast1',
  593. 'Neuroblast2'),],
  594. keep.extra.columns = TRUE,
  595. seqnames.field = 'seqnames',
  596. start.field = 'start',
  597. end.field = 'end',
  598. ignore.strand = TRUE)
  600. ### Split the GRanges object by modifier interval definition:
  602. FAM.sig.all_GR_Prog.list<- split(FAM.sig.all_GR_Prog,
  603. f = FAM.sig.all_GR_Prog$IntervalID2)
  605. ### Perform TF binding motif enrichment analysis on these differentially
  606. ### accessible peaks per modifier interval:
  608. EnrichedMotifsList<-list()
  609. EnrichedMotifsList<-lapply(X = FAM.sig.all_GR_Prog.list,
  610. FUN = function(x) {FindMotifs(object = reprocessed,
  611. features = unique(x$loc))})
  613. ### Make into a combined dataframe:
  615. EnrichedMotifs_Prog<-EnrichedMotifsList %>%
  616. bind_rows(.id = "IntervalID2")
  618. ### Filter for those with an adjusted p-value of less than 0.05:
  620. EnrichedMotifs_Prog<-EnrichedMotifs_Prog[EnrichedMotifs_Prog$p.adjust<0.05,]
  622. ### Add the percent motif in the test chromatin vs the percent motif in the
  623. ### background chromatin:
  625. EnrichedMotifs_Prog$percent.diff<-EnrichedMotifs_Prog$percent.observed-EnrichedMotifs_Prog$percent.background
  627. ### Convert the motif names into mouse gene names:
  629. convert_to_mouse_gene <- function(gene_name) {
  630. # Split the gene name by "::"
  631. parts <- unlist(strsplit(gene_name, "::"))
  632. # Convert each part to mouse gene format
  633. mouse_genes <- sapply(parts, function(part) {
  634. part <- tolower(part)
  635. part <- paste0(toupper(substring(part, 1, 1)), substring(part, 2))
  636. return(part)
  637. })
  638. ### Combine the parts back if there were multiple
  639. return(mouse_genes)
  640. }
  642. ### Apply the function to the GeneNames column
  643. MouseGeneNames <- gsub(pattern = '\\.',
  644. replacement = '',
  645. x = unique(unlist(sapply(EnrichedMotifs_Prog$motif.name,
  646. convert_to_mouse_gene))))
  648. ### There is a specific name that needs to be corrected here: EWSR1-FLI1
  649. MouseGeneNames[717]<-strsplit(MouseGeneNames[142], "-")[[1]][2]
  650. MouseGeneNames[717]<-paste0(toupper(substring(MouseGeneNames[717], 1, 1)),
  651. tolower(substring(MouseGeneNames[717], 2)))
  652. MouseGeneNames[142]<-strsplit(MouseGeneNames[142], "-")[[1]][1]
  654. ### Load in the "Neural Crest" cells from Zhao et al., 2022:
  656. nc1<-readRDS(file=paste0(path8,
  657. '2024-08-21_NeuralCrestOnly_Zhaoetal2022Dev_SCTv5IntegrateSeuratV5.RDS'))
  658. DimPlot(nc1,label=TRUE)
  659. nc1$time<-factor(x=nc1$time,
  660. levels = c("E9.5",
  661. "E10.5",
  662. "E11.5",
  663. "E13.5",
  664. "E15.5"))
  666. ### Isolate all motif gene names that are expressed in the "Neural Crest" cells
  667. ### above a pseudobulk-RNA expression of 0.1 by timepoint:
  669. exptable<-as.data.frame(AverageExpression(object = nc1,
  670. assays = 'RNA',
  671. features = MouseGeneNames,
  672. group.by = 'time'))
  673. exptablefilt_0.1<-filter_all(exptable, any_vars(. > 0.1))
  674. # exptablefilt_0.2<-filter_all(exptable, any_vars(. > 0.2))
  675. # exptablefilt_0.12<-filter_all(exptable, any_vars(. > 0.12))
  677. ### Filter the enriched motif dataset to those genes that are expressed:
  679. EnrichedMotifs_Prog_Exp <- EnrichedMotifs_Prog[sapply(EnrichedMotifs_Prog$motif.name,
  680. function(name) any(sapply(rownames(exptablefilt_0.1),
  681. grepl, name,ignore.case=TRUE))), ]
  683. ### Load in all result tables from each pipeline part:
  685. Sox10TFMotsModInt<-read.csv(file=paste0(savepath,'2024-10-08_Sox10TFMotifsConservedInModIntervals_ClosestFeature.csv'))
  687. Sox10TFMotsModInt10gene<-read.csv(file=paste0(savepath,'2024-10-08_Sox10TFMotifsConservedInModIntervals_Closest10Genes.csv'))
  689. Top2peaks10gene<-read.csv(file=paste0(savepath,'2024-09-04_Top2PeaksPerIntervalSet_Closest10Genes.csv'))
  691. WavefrontModInt<-read.csv(file=paste0(savepath,'2024-10-08_StavelyWavefrontDEG_InModifiers_All_Assays.csv'))
  693. dfcompleteUnfiltModInt_filt<-read.csv(file = paste0(savepath,'2024-10-02_FiltEnrichedCellChatGenesInModInts_netP.csv'))
  695. # FAMfilt<-read.csv(file = paste0(savepath,
  696. # '2024-09-30_FAM_Zhao_2022_ByCellTypeTimeCombined_pvalLT0.05.csv'))[2:8]
  698. DiffExpModInt<-read.csv(file=paste0(savepath,'2024-10-07_Zhao_FAM_DiffExpGenesGenesPerModInt.csv'))
  699. DiffExpModInt[c('CellType', 'Age')]<-stringr::str_split_fixed(DiffExpModInt$cluster,'_',2)
  700. # DiffExpModInt_Adj<-DiffExpModInt[DiffExpModInt$Age!='E15.5',]
  701. # DiffExpModInt_Adj<-subset(DiffExpModInt_Adj,subset=CellType %in%
  702. # unique(DiffExpModInt_Adj$CellType)[c(1,3:9,11:23)])
  703. DiffExpModInt_Adj<-subset(DiffExpModInt,subset=CellType %in%
  704. unique(DiffExpModInt$CellType)[c(1,3:9,11:24)])
  706. ### Convert the motif names into mouse gene names:
  708. convert_to_mouse_gene <- function(gene_name) {
  709. # Split the gene name by "::"
  710. parts <- unlist(strsplit(gene_name, "::"))
  711. # Convert each part to mouse gene format
  712. mouse_genes <- sapply(parts, function(part) {
  713. part <- tolower(part)
  714. part <- paste0(toupper(substring(part, 1, 1)), substring(part, 2))
  715. return(part)
  716. })
  717. ### Combine the parts back if there were multiple
  718. return(mouse_genes)
  719. }
  721. ### Apply the function to the GeneNames column
  722. MouseGeneNames <- gsub(pattern = '\\.',
  723. replacement = '',
  724. x = unique(unlist(sapply(EnrichedMotifs_Prog_Exp$motif.name,
  725. convert_to_mouse_gene))))
  727. ### Find n closest genes per DA peak:
  729. DAPeaksModInt<-FAM.sig.all[FAM.sig.all$cluster %in%
  730. c('Progenitors',
  731. 'Neuroblast1',
  732. 'Neuroblast2'),]
  734. # loci<-GenomicRanges::makeGRangesFromDataFrame(df = DAPeaksModInt,
  735. # keep.extra.columns = TRUE,
  736. # seqnames.field = 'seqnames',
  737. # start.field = 'start',
  738. # end.field = 'end',
  739. # ignore.strand = TRUE)
  740. DAPeaksModInt[c('seqnames', 'start','end')] <- stringr::str_split_fixed(DAPeaksModInt$loc, '-', 3)
  741. loci<-DAPeaksModInt[,8:10]
  742. loci$start<-as.numeric(loci$start)
  743. loci$end<-as.numeric(loci$end)
  744. genesobj<-as.data.frame(genes(TxDb.Mmusculus.UCSC.mm10.knownGene))
  746. # Retrieve gene symbols using AnnotationDbi
  747. gene_symbols <- select(org.Mm.eg.db,
  748. keys = genesobj$gene_id,
  749. columns = "SYMBOL",
  750. keytype = "ENTREZID")
  751. genesobj <- merge(genesobj,
  752. gene_symbols,
  753. by.x = "gene_id",
  754. by.y = "ENTREZID")
  756. # Function to find nearest n genes
  757. find_nearest_genes <- function(loci, genesobj, n) {
  758. result <- lapply(1:nrow(loci), function(i) {
  760. genesobjchrom<-genesobj[genesobj$seqnames==loci$seqnames[i],]
  762. # Calculate distances to each gene
  763. distances <- pmin(abs(genesobjchrom$start - loci$start[i]),
  764. abs(genesobjchrom$end - loci$start[i]),
  765. abs(genesobjchrom$start - loci$end[i]),
  766. abs(genesobjchrom$end - loci$end[i]))
  768. # Check if locus falls within gene
  769. within_gene <- (loci$start[i] >= genesobjchrom$start & loci$end[i] <= genesobjchrom$end) |
  770. (loci$start[i] <= genesobjchrom$start & loci$end[i] >= genesobjchrom$start) |
  771. (loci$start[i] <= genesobjchrom$end & loci$end[i] >= genesobjchrom$end)
  773. # Set distance to 0 if locus falls within gene
  774. distances[within_gene] <- 0
  776. # Get the indices of the nearest genes
  777. nearest_indices <- order(distances)[1:n]
  778. nearest_genes <- genesobjchrom[nearest_indices, ]
  779. data.frame(
  780. locus = paste0(loci$seqnames[i],'-',loci$start[i],'-',loci$end[i]),
  781. gene_ids = nearest_genes$SYMBOL,
  782. distances = distances[nearest_indices]
  783. )
  784. })
  785. do.call(rbind, result)
  786. }
  788. # Set the number of nearest genes to retrieve
  789. n <- 3
  791. # Find the nearest genes
  792. nearest_genes_df <- find_nearest_genes(loci, genesobj, n)
  794. # Display the results
  795. print(nearest_genes_df)
  797. nearest_genes_df_Sox10TFBM <- find_nearest_genes(Sox10TFMotsModInt[,3:5], genesobj, 10)
  799. ### Using the nearest 10 genes to both the conserved Sox10 TF binding motifs,
  800. ### the nearest 3 genes to the differentially accessible peaks, the
  801. ### differentially expressed genes from the gut cells across celltype/time, the
  802. ### differentially expressed genes from the ENCDC migrating wavefront, and the
  803. ### genes from the CellChat results (all of these modalities within modifier
  804. ### intervals), I will label the gene that overlap with the enriched TF motifs
  805. ### for figure-making purposes.
  807. genes<-unique(c(unique(c('St18','Pkhd1',
  808. 'Il17f','Mcm3',
  809. 'Adgrb3','Cyp7b1',
  810. 'Nlgn1','Mecom',
  811. 'Sox2','Slit2',
  812. 'Ppargc1a','Chic2',
  813. 'Pdgfra','Prdm8',
  814. 'Antxr2','Ranbp17',
  815. 'Tlx3','Tenm2',
  816. 'Bcas3','Tbx2',
  817. 'Cadps','Ptprg',
  818. 'Dach1','Fam204a')),
  819. unique(WavefrontModInt$symbol),
  820. unique(nearest_genes_df$gene_ids),
  821. c('Ppia','Col1a1','Lgals9','Ednrb','Nrg1'),
  822. unique(DiffExpModInt_Adj$mm10.kgXref.geneSymbol)))
  824. MouseGeneNamesSub<-intersect(MouseGeneNames,genes)
  826. EnrichedMotifs_Prog_Exp_sub <- EnrichedMotifs_Prog_Exp[sapply(EnrichedMotifs_Prog_Exp$motif.name,
  827. function(name) any(sapply(MouseGeneNamesSub,
  828. grepl, name,ignore.case=TRUE))), ]
  830. ### The above code includes genes Ebf1 and Tbx2. Unfortunately, the way the code
  831. ### is written, this fails to filter out Srebf1, Sox21, Tbx20, and Tbx21, because Ebf1
  832. ### and Tbx2 are substrings of those genes. I manually checked to see if each of
  833. ### these genes were in at least 1 other data modality, and only Srebf2 was in
  834. ### the differentially accessible peaks and in the differentially expressed
  835. ### genes. For the rest which were not in any other of the data modalities, I am
  836. ### manually removing them here:
  837. EnrichedMotifs_Prog_Exp_sub<-EnrichedMotifs_Prog_Exp_sub[EnrichedMotifs_Prog_Exp_sub$motif.name %in%
  838. setdiff(EnrichedMotifs_Prog_Exp_sub$motif.name,
  839. c('Tbx20','Tbx21',
  840. 'TBX21','TBX20',
  841. 'SOX21')),]
  844. EnrichedMotifs_Prog_Exp_subtoplot <- EnrichedMotifs_Prog_Exp_sub %>%
  845. separate(IntervalID2, into = c("IntervalID", "specificinterval"), sep = "_",
  846. extra = "merge")
  848. p4<-ggplot(data=EnrichedMotifs_Prog_Exp_subtoplot,
  849. aes(x=fold.enrichment,
  850. y=-log10(p.adjust),
  851. col=IntervalID,
  852. label=motif.name,
  853. size=percent.diff
  854. )) +
  855. geom_point() +
  856. theme_minimal() +
  857. theme(
  858. # legend.position = c(.93, .93),
  859. legend.title=element_text(size = 20,face='bold'),
  860. legend.text = element_text(size = 20,face='bold'),
  861. # legend.title=element_blank(),
  862. axis.line = element_line(colour = 'black', size = 1.1),
  863. axis.text.y = element_text(size = 10,face='bold'),
  864. axis.text.x = element_text(size = 10,face='bold'),
  865. axis.title = element_text(size = 15,face='bold')) +
  866. guides(colour=guide_legend(override.aes=list(alpha=1, size=3.25))) +
  867. ggrepel::geom_text_repel(show.legend =FALSE,
  868. fontface='bold',
  869. color='black')
  870. p4
  871. ggsave(device = 'jpg',dpi = 600,
  872. plot = p4,units = 'in',width = 11,height = 7,
  873. filename = paste0(savepath,
  874. "Images/2025-02-11_DiffAccChromTFsEnrichProgNeurblst12_",
  875. 'OverlapWithOtherModalities',
  876. "_snATAC-seq_Phox2b_predicted.id_1SidedVolc.jpg"))
  878. p4<-ggplot(data=EnrichedMotifs_Prog_Exp,
  879. aes(x=fold.enrichment,
  880. y=-log10(p.adjust),
  881. col=assay,
  882. label=motif.name,
  883. size=percent.diff
  884. )) +
  885. geom_point() +
  886. theme_minimal() +
  887. theme(
  888. # legend.position = c(.93, .93),
  889. legend.title=element_text(size = 20,face='bold'),
  890. legend.text = element_text(size = 20,face='bold'),
  891. # legend.title=element_blank(),
  892. axis.line = element_line(colour = 'black', size = 1.1),
  893. axis.text.y = element_text(size = 10,face='bold'),
  894. axis.text.x = element_text(size = 10,face='bold'),
  895. axis.title = element_text(size = 15,face='bold')) +
  896. guides(colour=guide_legend(override.aes=list(alpha=1, size=3.25))) +
  897. ggrepel::geom_text_repel(show.legend =FALSE,
  898. fontface='bold',
  899. color='black')
  900. p4
  901. ggsave(device = 'jpg',dpi = 600,
  902. plot = p4,units = 'in',width = 11,height = 7,
  903. filename = paste0(savepath,
  904. "Images/2025-07-08_DiffAccChromTFsEnrichProgNeurblst12_",
  905. 'ExpressedEnENCDCs',
  906. "_snATAC-seq_Phox2b_predicted.id_1SidedVolc.jpg"))
  909. write.csv(x = EnrichedMotifs_Prog,
  910. file = paste0(savepath,'/2024-10-14_EnrichedMotifs_Progen+NB12_DAC-derivedMotifEnrichment_InModInts_byIntervalID2.csv'),
  911. quote = FALSE,row.names = FALSE)
  912. write.csv(x = EnrichedMotifs_Prog_Exp,
  913. file = paste0(savepath,'/2024-10-14_EnrichedMotifs_ExprsdNeurCrst_Progen+NB12_DAC-derivedMotifEnrichment_InModInts_byIntervalID2.csv'),
  914. quote = FALSE,row.names = FALSE)
  917. ### Find differentially accessible peaks that overlap with Sox10 binding sites
  918. ### within modifier intervals, if any:
  920. Sox10TFMotifsInModInts.all.cf<-read.csv(file = paste0(savepath,
  921. '2024-10-08_Sox10TFMotifsConservedInModIntervals_ClosestFeature.csv'))
  923. ### Binding sites are smaller than the peaks, so I am using the below order to
  924. ### find the overlap:
  926. locilist<-list()
  927. for(i in seq_along(FAM.sig.all$seqnames)){
  928. locilist[[i]]<-Sox10TFMotifsInModInts.all.cf[Sox10TFMotifsInModInts.all.cf$seqnames==FAM.sig.all$seqnames[i] &
  929. Sox10TFMotifsInModInts.all.cf$start>=FAM.sig.all$start[i] &
  930. Sox10TFMotifsInModInts.all.cf$end<=FAM.sig.all$end[i],]
  931. }
  932. FAM_inSox10Gopinath<-do.call('rbind',locilist)
  934. unique(FAM_inSox10Gopinath$V2)
  936. ### Two differentially accessible loci overlap with Sox10 TF binding motifs
  937. ### c('chr11-33400347-33400362','chr14-97891443-97891462'). These overlap with
  938. ### Ranbp17 and Dach1 introns. However, these do not have many fragments on
  939. ### These peaks (normalized frags 0-20 in the CoveragePlot) are very small.
  941. Sox10TFMotifsInModInts.all.cf[,4]<-as.numeric(Sox10TFMotifsInModInts.all.cf[,4])
  942. Sox10TFMotifsInModInts.all.cf[,5]<-as.numeric(Sox10TFMotifsInModInts.all.cf[,5])
  944. ### Plot binding motif loci and Dach1 as a whole:
  946. CoveragePlot(
  947. object = reprocessed,links = FALSE,
  948. region = 'chr14-97891443-97891462',extend.downstream = 10000,extend.upstream = 10000) /
  949. Signac::PeakPlot(object = reprocessed,region = 'chr14-97891443-97891462',
  950. peaks = GRanges(Sox10TFMotifsInModInts.all.cf[,3:5]),
  951. extend.downstream = 10000,extend.upstream = 10000)+
  952. patchwork::plot_layout(heights = c(5, 0.5))
  954. CoveragePlot(
  955. object = reprocessed,links = FALSE,
  956. region = 'chr14-97786853-98169543',extend.downstream = 10000,extend.upstream = 10000) /
  957. Signac::PeakPlot(object = reprocessed,region = 'chr14-97786853-98169543',
  958. peaks = GRanges(Sox10TFMotifsInModInts.all.cf[,3:5]),
  959. extend.downstream = 10000,extend.upstream = 10000)+
  960. patchwork::plot_layout(heights = c(5, 0.5))
  962. ### Plot binding motif loci and Ranbp17 as a whole:
  964. CoveragePlot(
  965. object = reprocessed,links = FALSE,
  966. region = 'chr11-33400347-33400362',extend.downstream = 10000,extend.upstream = 10000) /
  967. Signac::PeakPlot(object = reprocessed,region = 'chr11-33400347-33400362',
  968. peaks = GRanges(Sox10TFMotifsInModInts.all.cf[,3:5]),
  969. extend.downstream = 10000,extend.upstream = 10000)+
  970. patchwork::plot_layout(heights = c(5, 0.5))
  972. CoveragePlot(
  973. object = reprocessed,links = FALSE,
  974. region = 'chr11-33211795-33513746',extend.downstream = 10000,extend.upstream = 10000) /
  975. Signac::PeakPlot(object = reprocessed,region = 'chr11-33211795-33513746',
  976. peaks = GRanges(Sox10TFMotifsInModInts.all.cf[,3:5]),
  977. extend.downstream = 10000,extend.upstream = 10000)+
  978. patchwork::plot_layout(heights = c(5, 0.5))
  980. ### These peaks are very small and may not be important. Will investigate these
  981. ### genes further but this is not evidence that there are conserved Sox10
  982. ### binding sites within these differentially accessible peaks.
  984. write.csv(x = FAM_inSox10Gopinath,
  985. file = paste0(savepath,'2024-10-14_Sox10GopinathBS_Within_Phox2bCFPsnATAC-seqWT_DAPeaks_predicted.id_InModInts+CF.csv'),
  986. quote = FALSE,row.names = FALSE)
  988. saveRDS(reprocessed,paste0('D:/ModifierIntervalCodeSeuratV5/DataFiles/2024-10-15_snATAC-seq_16.5dpcPhox2bH2B-CFP+WholeGutPooled_Reprocessed.RDS'))
  990. sessionInfo()
  992. # R version 4.4.1 (2024-06-14 ucrt)
  993. # Platform: x86_64-w64-mingw32/x64
  994. # Running under: Windows 11 x64 (build 22631)
  995. #
  996. # Matrix products: default
  997. #
  998. #
  999. # locale:
  1000. # [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8
  1001. # [4] LC_NUMERIC=C LC_TIME=English_United States.utf8
  1002. #
  1003. # time zone: America/Chicago
  1004. # tzcode source: internal
  1005. #
  1006. # attached base packages:
  1007. # [1] stats4 stats graphics grDevices utils datasets methods base
  1008. #
  1009. # other attached packages:
  1010. # [1] org.Mm.eg.db_3.19.1 ggVennDiagram_1.5.2
  1011. # [3] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0 GenomicFeatures_1.56.0
  1012. # [5] AnnotationDbi_1.66.0 Biobase_2.64.0
  1013. # [7] lubridate_1.9.3 forcats_1.0.0
  1014. # [9] stringr_1.5.1 dplyr_1.1.4
  1015. # [11] purrr_1.0.2 readr_2.1.5
  1016. # [13] tidyr_1.3.1 tibble_3.2.1
  1017. # [15] ggplot2_3.5.1 tidyverse_2.0.0
  1018. # [17] BSgenome.Mmusculus.UCSC.mm10_1.4.3 BSgenome_1.72.0
  1019. # [19] rtracklayer_1.64.0 BiocIO_1.14.0
  1020. # [21] Biostrings_2.72.1 XVector_0.44.0
  1021. # [23] JASPAR2024_0.99.6 BiocFileCache_2.12.0
  1022. # [25] dbplyr_2.5.0 TFBSTools_1.42.0
  1023. # [27] GenomicRanges_1.56.1 GenomeInfoDb_1.40.1
  1024. # [29] IRanges_2.38.0 S4Vectors_0.42.0
  1025. # [31] BiocGenerics_0.50.0 cowplot_1.1.3
  1026. # [33] Signac_1.13.0 Seurat_5.1.0
  1027. # [35] SeuratObject_5.0.2 sp_2.1-4
  1028. #
  1029. # loaded via a namespace (and not attached):
  1030. # [1] fs_1.6.4 matrixStats_1.3.0 spatstat.sparse_3.0-3
  1031. # [4] bitops_1.0-7 enrichplot_1.24.0 DirichletMultinomial_1.46.0
  1032. # [7] HDO.db_0.99.1 httr_1.4.7 RColorBrewer_1.1-3
  1033. # [10] tools_4.4.1 sctransform_0.4.1 utf8_1.2.4
  1034. # [13] R6_2.5.1 lazyeval_0.2.2 uwot_0.2.2
  1035. # [16] withr_3.0.0 gridExtra_2.3 progressr_0.14.0
  1036. # [19] cli_3.6.2 textshaping_0.4.0 spatstat.explore_3.2-7
  1037. # [22] fastDummies_1.7.3 scatterpie_0.2.3 labeling_0.4.3
  1038. # [25] spatstat.data_3.0-4 ggridges_0.5.6 pbapply_1.7-2
  1039. # [28] yulab.utils_0.1.4 Rsamtools_2.20.0 systemfonts_1.1.0
  1040. # [31] DOSE_3.30.1 R.utils_2.12.3 parallelly_1.37.1
  1041. # [34] plotrix_3.8-4 limma_3.60.2 rstudioapi_0.16.0
  1042. # [37] RSQLite_2.3.7 TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 generics_0.1.3
  1043. # [40] gridGraphics_0.5-1 gtools_3.9.5 ica_1.0-3
  1044. # [43] spatstat.random_3.2-3 GO.db_3.19.1 Matrix_1.7-0
  1045. # [46] fansi_1.0.6 abind_1.4-5 R.methodsS3_1.8.2
  1046. # [49] lifecycle_1.0.4 yaml_2.3.8 SummarizedExperiment_1.34.0
  1047. # [52] gplots_3.1.3.1 qvalue_2.36.0 SparseArray_1.4.8
  1048. # [55] Rtsne_0.17 grid_4.4.1 blob_1.2.4
  1049. # [58] promises_1.3.0 crayon_1.5.2 pwalign_1.0.0
  1050. # [61] miniUI_0.1.1.1 lattice_0.22-6 annotate_1.82.0
  1051. # [64] KEGGREST_1.44.0 pillar_1.9.0 fgsea_1.30.0
  1052. # [67] boot_1.3-30 rjson_0.2.21 future.apply_1.11.2
  1053. # [70] codetools_0.2-20 fastmatch_1.1-4 leiden_0.4.3.1
  1054. # [73] glue_1.7.0 ggfun_0.1.5 data.table_1.15.4
  1055. # [76] treeio_1.28.0 vctrs_0.6.5 png_0.1-8
  1056. # [79] spam_2.10-0 gtable_0.3.5 poweRlaw_0.80.0
  1057. # [82] cachem_1.1.0 S4Arrays_1.4.1 mime_0.12
  1058. # [85] tidygraph_1.3.1 pracma_2.4.4 survival_3.6-4
  1059. # [88] RcppRoll_0.3.0 statmod_1.5.0 fitdistrplus_1.1-11
  1060. # [91] ROCR_1.0-11 nlme_3.1-164 ggtree_3.12.0
  1061. # [94] bit64_4.0.5 filelock_1.0.3 RcppAnnoy_0.0.22
  1062. # [97] irlba_2.3.5.1 KernSmooth_2.23-24 colorspace_2.1-0
  1063. # [100] seqLogo_1.70.0 DBI_1.2.3 tidyselect_1.2.1
  1064. # [103] bit_4.0.5 compiler_4.4.1 curl_5.2.1
  1065. # [106] DelayedArray_0.30.1 plotly_4.10.4 shadowtext_0.1.3
  1066. # [109] scales_1.3.0 caTools_1.18.2 ChIPseeker_1.40.0
  1067. # [112] lmtest_0.9-40 digest_0.6.35 goftest_1.2-3
  1068. # [115] spatstat.utils_3.0-4 presto_1.0.0 motifmatchr_1.26.0
  1069. # [118] htmltools_0.5.8.1 pkgconfig_2.0.3 MatrixGenerics_1.16.0
  1070. # [121] fastmap_1.2.0 rlang_1.1.4 htmlwidgets_1.6.4
  1071. # [124] UCSC.utils_1.0.0 shiny_1.8.1.1 farver_2.1.2
  1072. # [127] zoo_1.8-12 jsonlite_1.8.8 BiocParallel_1.38.0
  1073. # [130] GOSemSim_2.30.0 R.oo_1.26.0 RCurl_1.98-1.14
  1074. # [133] magrittr_2.0.3 GenomeInfoDbData_1.2.12 ggplotify_0.1.2
  1075. # [136] dotCall64_1.1-1 patchwork_1.2.0 munsell_0.5.1
  1076. # [139] Rcpp_1.0.12 ape_5.8 viridis_0.6.5
  1077. # [142] reticulate_1.37.0 stringi_1.8.4 ggraph_2.2.1
  1078. # [145] zlibbioc_1.50.0 MASS_7.3-60.2 plyr_1.8.9
  1079. # [148] parallel_4.4.1 listenv_0.9.1 ggrepel_0.9.5
  1080. # [151] deldir_2.0-4 CNEr_1.40.0 graphlayouts_1.1.1
  1081. # [154] splines_4.4.1 tensor_1.5 hms_1.1.3
  1082. # [157] igraph_2.0.3 spatstat.geom_3.2-9 RcppHNSW_0.6.0
  1083. # [160] pkgload_1.3.4 reshape2_1.4.4 TFMPvalue_0.0.9
  1084. # [163] XML_3.99-0.16.1 tweenr_2.0.3 tzdb_0.4.0
  1085. # [166] httpuv_1.6.15 RANN_2.6.1 polyclip_1.10-6
  1086. # [169] future_1.33.2 scattermore_1.2 ggforce_0.4.2
  1087. # [172] xtable_1.8-4 restfulr_0.0.15 tidytree_0.4.6
  1088. # [175] RSpectra_0.16-1 later_1.3.2 viridisLite_0.4.2
  1089. # [178] ragg_1.3.2 aplot_0.2.3 memoise_2.0.1
  1090. # [181] GenomicAlignments_1.40.0 cluster_2.1.6 timechange_0.3.0
  1091. # [184] globals_0.16.3

2025-02-12_ModifierInterval_CandidateGenePipeline_Part7_V2.R, under CC-BY-4.0 · at the source

Overview

Authors: Joseph T Benthal1, Justin A Avila2,3, Jeffrey R Smith4, E Michelle Southard-Smith1,2,4
ORCID iDs: E Michelle Southard-Smith
  1. Program in Human Genetics, Vanderbilt University, Nashville, Tennessee, United States of America
  2. Vanderbilt Brain Institute, Vanderbilt University, Nashville, Tennessee, United States of America
  3. Stowers Institute for Medical Research, Kansas City, Missouri, United States of America
  4. Genetic Medicine, Vanderbilt University School of Medicine, Nashville, Tennessee, United States of America
Institutions: Vanderbilt University (United States); Stowers Institute for Medical Research (United States)
Journal: PLoS computational biology, volume 22, issue 7, article e1014424
Dates: received 3 November 2025; accepted 9 June 2026; published online 6 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014424 · PMID 42406857 · PMCID PMC13372245 · OpenAlex W7167545869
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), mouse (organism)
Methods: Statistics
MeSH: Enteric Nervous System*, Hirschsprung Disease*, SOXE Transcription Factors*, Animals, Computational Biology, Disease Models, Animal, Female, Genetic Predisposition to Disease, Genome-Wide Association Study, Humans, Male, Mice, Multiomics, Neural Crest, Pedigree (* major topic)
Topic: Congenital gastrointestinal and neural anomalies (Surgery, Medicine), according to OpenAlex
Funding: National Institute of Diabetes and Digestive and Kidney Diseases (1F31DK137637-01, R01 DK127178, R01 DK60047, R01-DK127178-03S1); Howard Hughes Medical Institute (GT11505); National Institute of Child Health and Human Development (T32-HD007502); NIDDK NIH HHS (R01 DK127178, P30 DK058404, F31 DK137637, R01 DK060047); National Institute of Diabetes and Kidney Disease (P30-DK058404); NICHD NIH HHS (P50 HD103537, T32 HD007502); Eunice Kennedy Shriver National Institute of Child Health and Human Development (T32-HD007502)
Citations: not cited yet (Europe PMC); 96 references in the paper

Abstract

Hirschsprung disease (HSCR) is characterized by absence of enteric ganglia (aganglionosis) along variable lengths of the distal intestine. This disorder results from deficient colonization of fetal intestine by enteric neural crest-derived cells (ENCDCs). HSCR exhibits complex, multifactorial inheritance with penetrance and severity varying widely even within families. SOX10 is among causal genes that predispose to aganglionosis. Yet, how gene interactions influence severity of HSCR aganglionosis is not understood. Prior mapping of aganglionosis modifiers was achieved in a standard F1-intercross utilizing the Sox10Dom HSCR mouse model. Here we deploy a novel strategy of genotyping an extended pedigree pedigree of Sox10Dom mice on a mixed genetic background. GWAS in this pedigree points to novel aganglionosis modifier intervals with replication and refinement of prior modifier regions. Complementary omics analysis of the developing Enteric Nervous System (ENS) enabled identification of multiple high-priority candidate genes within these modifier intervals based on gene expression, chromatin accessibility, and presence of conserved SOX10 binding motifs. We implemented a prioritization pipeline for ranking potential modifiers that generated candidate lists including several well-known for effects on ENS development as well as multiple novel genes. Among the novel genes, Dach1 ranked as a top priority candidate gene for modifying migration of ENCDCs and thus influencing aganglionosis severity. The results identify genome intervals with intrinsic genes that are logical candidates for modifying Sox10Dom aganglionosis severity. We also note that several human orthologs to aganglionosis modifier candidate genes are within linkage disequilibrium blocks containing genetic variants associated with human gut motility disorders, which offers opportunity for gaining biological insight into human HSCR severity.

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 20 matches between paragraphs and lines of code.

Zenodo 20671852

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: R (21)
Size: 31 files, 21 scripts
Software Heritage: not checked
Found in: “Data Availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (18 files), ggplot2 (16 files), Seurat (14 files), patchwork (6 files), cowplot (4 files), ggpubr (3 files), data.table (2 files), Harmony (2 files), clusterProfiler (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
21 files
At the source: doi.org/10.5281/zenodo.20671852

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;
  • 21 scripts, each with its path and the digest of its content;
  • 20 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 Code, fragment files, and snATAC-seq Seurat Object are shared at Zenodo: https://doi.org/10.5281/zenodo.20671852 SnATAC-seq data is accessible at GEO under accession GSE298406 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE298406).

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

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 15 MeSH terms, 7 funders, 94 references.

Cite

This paper

Benthal, J. T., Avila, J. A., Smith, J. R., & Southard-Smith, E. M. (2026). Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice identifies multiple high confidence candidate modifiers of Enteric Nervous System development. PLoS computational biology, 22(7), e1014424. https://doi.org/10.1371/journal.pcbi.1014424

BibTeX

@article{benthal2026combinatorial,
author = {Benthal, Joseph T and Avila, Justin A and Smith, Jeffrey R and Southard-Smith, E Michelle},
title = {{Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice identifies multiple high confidence candidate modifiers of Enteric Nervous System development}},
journal = {PLoS computational biology},
year = {2026},
month = jul,
volume = {22},
number = {7},
pages = {e1014424},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014424},
url = {https://doi.org/10.1371/journal.pcbi.1014424},
pmid = {42406857},
pmcid = {PMC13372245}
}

RIS

TY - JOUR
AU - Benthal, Joseph T
AU - Avila, Justin A
AU - Smith, Jeffrey R
AU - Southard-Smith, E Michelle
TI - Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice identifies multiple high confidence candidate modifiers of Enteric Nervous System development
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/07/06
VL - 22
IS - 7
SP - e1014424
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014424
UR - https://doi.org/10.1371/journal.pcbi.1014424
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014424",
"type": "article-journal",
"title": "Combinatorial multiomic analysis from a pedigree of Sox10Dom Hirschsprung mice identifies multiple high confidence candidate modifiers of Enteric Nervous System development",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Benthal",
"given": "Joseph T"
},
{
"family": "Avila",
"given": "Justin A"
},
{
"family": "Smith",
"given": "Jeffrey R"
},
{
"family": "Southard-Smith",
"given": "E Michelle"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "7",
"page": "e1014424",
"DOI": "10.1371/journal.pcbi.1014424",
"PMID": "42406857",
"PMCID": "PMC13372245",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014424",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
6
]
]
}
}

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.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, clusterProfiler, Seurat, 7 other tools, genetics / omics, 4 references
[2] doi:10.1038/s41597-026-07185-4 [code]
A multi-center cross-platform single-cell multimodal atlas of the mouse cerebral cortex.
Journal: Scientific data
In common: Harmony, clusterProfiler, Seurat, 6 other tools, genetics / omics, mouse, 4 references
[3] doi:10.1038/s41467-026-76232-w [code]
Th17 effector cytokines induce shared and distinct microglial and endothelial cell responses in a mouse model for post-streptococcal encephalitis.
Journal: Nature communications
In common: Harmony, clusterProfiler, Seurat, 6 other tools, genetics / omics, mouse, 3 references
[4] 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: Harmony, clusterProfiler, Seurat, 7 other tools, 2 references
[5] doi:10.1038/s41398-026-04200-5 [code]
Postmortem brain single-nucleus and bulk gene expression analyses identify shared and distinct abnormalities in bipolar disorder and major depressive disorder.
Journal: Translational psychiatry
In common: Harmony, clusterProfiler, Seurat, 6 other tools, genetics / omics, 3 references
[6] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, clusterProfiler, Seurat, 7 other tools, genetics / omics, mouse, 1 reference
[7] doi:10.1038/s41514-026-00397-3 [code]
Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8&lt;sup&gt;+&lt;/sup&gt; T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.
Journal: npj aging
In common: Harmony, Seurat, cowplot, 6 other tools, 3 references
[8] doi:10.1101/gr.281113.125 [code]
Single-nucleus multiomic profiling of the aging mouse substantia nigra reveals conserved gene alterations linked to Parkinson's disease.
Journal: Genome research
In common: clusterProfiler, Seurat, cowplot, 6 other tools, genetics / omics, mouse, 2 references
[9] doi:10.1186/s12967-026-08266-z [code]
Single-cell multi-omic integration analysis prioritizes druggable genes and reveals cell-type-specific causal effects in glioblastomagenesis.
Journal: Journal of translational medicine
In common: clusterProfiler, Seurat, cowplot, 6 other tools, genetics / omics, 2 references
[10] doi:10.1016/j.cell.2026.05.026 [code]
The critical role of the endogenous immune compartment after CAR T cell therapy in recurrent GBM.
Journal: Cell
In common: Harmony, clusterProfiler, Seurat, 6 other tools, genetics / omics, 2 references

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.