OSCR

Identification and validation of central amygdala FGFR1 as a therapeutic target for alcohol use disorder using single-nucleus sequencing in rats.

Code ↔ Paper

4 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 4 matches
  1. [1] § Methods › snRNA-seq data processing › Co-expression network analysis ↔ master.R, lines 783–862 · score 0.94 · soft power, ConstructNetwork, ModuleConnectivity, ModuleEigengenes, SampleID, WGCNA
  2. [2] § Methods › snRNA-seq data processing › Quality control and clustering ↔ master.R, lines 89–145 · score 0.88 · SCTransform, percent.mt, percent.ribo, sample layers, M22, M33
  3. [3] § Methods › snRNA-seq data processing › Differential expression analysis ↔ master.R, lines 1615–1684 · score 0.77 · distinct_test, min_non_zero_cells, name_cluster, name_sample, cpm, matrix
  4. [4] § Methods › snRNA-seq data processing › Subclustering of GABAergic cells ↔ master.R, lines 432–491 · score 0.69 · percent.mt, percent.ribo, sample layers, split, subset, Subclustering

Paper

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

The paper is loaded when this pane is shown.

The authors' code

R · 2,362 lines · 59 KB · MIT · 4 matches

  1. library(Seurat)
  2. library(SeuratDisk)
  4. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  7. library(BPCells)
  8. adata.data = open_matrix_anndata_hdf5("D:/Leon/snRNAseq/high_depth_processing/Heilig_aggr_noF1_nonormalize_highdepth/outs/count/high_processing.h5ad")
  10. write_matrix_dir(
  11. mat = adata.data,
  12. dir = "D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth/counts"
  13. )
  17. adata.mat = open_matrix_dir(dir = "D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth/counts")
  20. #library(biomaRt)
  21. #ensembl = useMart("ensembl", dataset = "rnorvegicus_gene_ensembl")
  22. #bm = getBM(attributes = c("ensembl_gene_id", "external_gene_name"), values = rownames(adata.mat), mart = ensembl)
  24. #symbols = bm$external_gene_name[match(rownames(adata.mat), bm$ensembl_gene_id)]
  25. #rownames(adata.mat) = symbols
  28. adata = CreateSeuratObject(counts = adata.mat, min.cells = 3, min.features = 200)
  30. adata
  33. SampleID_barcode <- read.csv("D:/Leon/snRNAseq/high_depth_processing/Heilig_aggr_noF1_nonormalize_highdepth/outs/SampleID.csv")
  34. meta_file<-read.csv("D:/Leon/snRNAseq/high_depth_processing/Heilig_aggr_noF1_nonormalize_highdepth/outs/aggregation.csv")
  36. barcode_sampleid<-merge(SampleID_barcode,meta_file,by.x="SampleID",by.y="sample_id")
  37. rownames(barcode_sampleid)<-barcode_sampleid$Barcode
  38. adata<-AddMetaData(adata,barcode_sampleid)
  40. adata = subset(adata, subset = Sex == "Male")
  42. adata[["RNA"]]$counts <- as(object = adata[["RNA"]]$counts, Class = "dgCMatrix")
  43. Features(adata)
  44. saveRDS(adata, "males_unprocessed_high.Rds")
  46. #-----------------------------------------------QC---------------------------------------
  47. library(Seurat)
  48. library(BPCells)
  49. library(scuttle)
  51. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  52. adata = readRDS("males_unprocessed_high.Rds")
  53. adata
  56. library(BiocParallel)
  57. library(scDblFinder)
  59. #adata[["RNA"]]$counts <- as(object = adata[["RNA"]]$counts, Class = "dgCMatrix")
  61. adata = as.SingleCellExperiment(adata)
  63. adata = scDblFinder(adata, samples="SampleID")
  65. table(adata$scDblFinder.class)
  67. adata = as.Seurat(adata, data = NULL)
  68. adata
  69. gc()
  70. Idents(adata) <- "scDblFinder.class"
  71. adata = subset(adata, idents = "singlet") #remove doublets
  72. adata
  75. adata[['percent.mt']] = PercentageFeatureSet(adata, pattern = "^Mt-")
  76. gc()
  77. adata[['percent.ribo']] = PercentageFeatureSet(adata, pattern = "(^Rpl|^Rps|^Mrp)")
  78. gc()
  80. #adata.sce = as.SingleCellExperiment(adata)
  81. #low.total = isOutlier(adata.sce$nCount_RNA, type="lower", log=TRUE, batch = adata.sce$SampleID)
  82. #summary(low.total)
  84. #discard.mito = isOutlier(adata.sce$percent.mt, type="higher", batch=adata.sce$SampleID)
  85. #discard.ribo = isOutlier(adata.sce$percent.ribo, type="higher", batch=adata.sce$SampleID)
  86. #discard.sum = isOutlier(adata.sce$nCount_RNA, type="lower", log=TRUE, batch=adata.sce$SampleID)
  87. #discard.feature = isOutlier(adata.sce$nFeature_RNA, type="lower", log=TRUE, batch=adata.sce$SampleID)
  89. #qc.stats = perCellQCMetrics(adata.sce)
  91. VlnPlot(adata, 'nFeature_RNA', group.by = 'SampleID', pt.size = 0)
  92. VlnPlot(adata, 'nCount_RNA', group.by = 'SampleID', pt.size = 0)
  93. VlnPlot(adata, 'percent.mt', group.by = 'SampleID', pt.size = 0)
  94. VlnPlot(adata, 'percent.ribo', group.by = 'SampleID', pt.size = 0)
  96. FeatureScatter(adata, 'nCount_RNA', 'percent.mt')
  97. FeatureScatter(adata, 'nCount_RNA', 'nFeature_RNA')
  99. adata = adata[, !adata$SampleID=='M22'] #removing bad samples
  100. adata = adata[, !adata$SampleID=='M33']
  101. adata = adata[, !adata$SampleID=='M8']
  102. gc()
  104. adata = subset(adata, subset = nFeature_RNA < 7000)
  105. adata = subset(adata, subset = percent.mt < 5)
  106. adata = subset(adata, subset = percent.ribo < 2.5)
  107. gc()
  108. VlnPlot(adata, 'nFeature_RNA', group.by = 'SampleID', pt.size = 0)
  109. VlnPlot(adata, 'percent.mt', group.by = 'SampleID', pt.size = 0)
  110. VlnPlot(adata, 'percent.ribo', group.by = 'SampleID', pt.size = 0)
  114. adata
  116. FeatureScatter(adata, 'nCount_RNA', 'percent.mt')
  117. FeatureScatter(adata, 'nCount_RNA', 'nFeature_RNA')
  120. saveRDS(adata, "males_filtered_high.rds")
  123. #---------------------------normalization----------------------------------
  125. adata = readRDS("males_filtered_high.rds")
  127. table(adata$Groups)
  130. counts = GetAssayData(adata, assay = "RNA")
  131. #remove mito and ribo genes
  132. discard = which(grepl("^Rpl|^Rps|^Mt-|^Mrp", rownames(counts)))
  133. counts = counts[-c(discard),]
  134. adata = subset(adata, features = rownames(counts))
  135. adata
  136. gc()
  138. adata[["RNA"]] <- split(adata[["RNA"]], f = adata$SampleID) #split the counts into sample-layers
  139. gc()
  142. adata <- SCTransform(adata, vars.to.regress = c("percent.mt", "percent.ribo"))
  144. gc()
  145. saveRDS(adata, "filtered_SCT_high.rds")
  148. #---------------------------------------------- dim reduction and clustering -----------------------------------
  149. library(Seurat)
  150. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  153. adata = readRDS("filtered_SCT_high.rds")
  155. gc()
  156. adata
  158. adata = RunPCA(adata)
  159. gc()
  160. ElbowPlot(adata,ndims = 50)
  165. gc()
  166. adata = RunUMAP(adata, dims = 1:50, metric = 'euclidean', reduction.name = 'euclidean_umap')
  167. gc()
  168. #adata = RunUMAP(adata, dims = 1:50, metric = 'correlation', reduction.name = 'correlation_umap')
  170. adata
  171. gc()
  174. DimPlot(adata, reduction = "euclidean_umap", group.by = "SampleID")
  175. DimPlot(adata, reduction = "euclidean_umap", group.by = "Groups")
  177. gc()
  179. adata <- FindNeighbors(adata, dims = 1:50, reduction = "pca")
  180. gc()
  181. adata <- FindClusters(adata, resolution = 0.1, cluster.name = "SCT_clusters")
  182. gc()
  184. DimPlot(adata, reduction = "euclidean_umap", group.by = "SCT_clusters")
  185. gc()
  187. adata = PrepSCTFindMarkers(adata)
  188. gc()
  190. adata.markers = FindAllMarkers(adata, only.pos = TRUE)
  194. saveRDS(adata, "high_clustered.rds")
  195. write.csv(adata.markers, "high_markers.csv")
  198. #---------------------------anotation------------------------------------
  200. library(Seurat)
  201. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  204. adata = readRDS("high_clustered.rds")
  205. markers = read.csv('high_markers.csv')
  206. gc()
  208. DimPlot(adata, reduction = "euclidean_umap", group.by = "SCT_clusters", label = TRUE)
  210. markers[markers$gene == "Slc17a7",]
  212. gc()
  213. VlnPlot(adata, features = c(""))
  214. FeaturePlot(adata, features = c("Ndst4"), reduction = "euclidean_umap")
  216. VlnPlot(adata, features = c("Slc17a7"))
  221. #adata = RenameIdents(adata,
  222. "0" = "Oligendrocytes", #Mog, Mobp, Mag
  223. "1"= "Gabaergic", #Rbfox3, Gad1, Gad2,
  224. "2" = "Gabaergic", #Rbfox3, Gad1, Gad2,
  225. "3" = "Hbb-b1_Hbb-a2_Blood?", #,
  226. "4" = "OPC", #Pdgfra, Epn2, Neu4, Pcdh15
  227. "5" = "Microglia", # Itgam, Cx3cr1, P2ry12, Tmem119
  228. "6" = "Gabaergic", #Rbfox3, Gad2, Foxp2,
  229. "7" = "Astrocytes", #Slc1a2, Gfap, Agt, Aqp4,
  230. "8" = "Gabaergic", #Rbfox3, Gad1, Foxp2,
  231. "9" = "Gabaergic", #Rbfox3, Gad1, Gad2,
  232. "10" = "Gabaergic", #Gad1, Gad2,
  233. "11" = "Glutaminergic_Slc17a6", #Rbfox3, Slc17a6,
  234. "12" = "Endothelial", #Cd93, Vwf, Emcn, Flt1
  235. "13" = "Cholinergic", #Slc5a7, Acly,
  236. "14" = "Glutaminergic_Slc17a7", #Slc17a7, Rbfox3,
  237. "15" = "Gabaergic", #Rbfox3, Gad1, Gad2,
  238. "16" = "Microglia", #Itgam, Cx3cr1, P2ry12, Tmem119
  239. "17" = "Microglia ", #Itgam, P2ry12,
  241. )
  242. current.cluster.ids = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17)
  243. new.cluster.ids = c("Oligendrocytes","Gabaergic","Gabaergic", "Hbb-b1_Hbb-a2_Blood?",
  244. "OPC","Microglia","Gabaergic","Astrocytes","Gabaergic",
  245. "Gabaergic","Gabaergic","Glutaminergic_Slc17a6","Endothelial",
  246. "Cholinergic","Glutaminergic_Slc17a7","Gabaergic","Microglia","Microglia")
  248. library(plyr)
  249. adata$cell_type = plyr::mapvalues(x = adata$SCT_clusters, from = current.cluster.ids, to = new.cluster.ids)
  253. DimPlot(adata, group.by='cell_type', reduction = "euclidean_umap")
  256. #---------------------- cluster exploration------------------------
  258. Idents(adata) = adata$cell_type
  260. RidgePlot(adata, features = "percent.ribo")
  261. RidgePlot(adata, features = "percent.mt")
  262. RidgePlot(adata, features = "nCount_RNA")
  263. RidgePlot(adata, features = "nFeature_RNA")
  265. table(adata$cell_type, adata$SampleID)
  266. DimPlot(adata, group.by='SampleID', reduction = "euclidean_umap")
  274. features = c("Mog", "Mobp", "Mag", "Gad1", "Gad2", "Hbb-b1", "Hba-a2",
  275. "Pdgfra", "Epn2", "Pcdh15", "Csf1r", "Cx3cr1", "P2ry12", "Tmem119",
  276. "Satb1", "Slc17a6","Slc17a7", "Cd93", "Vwf", "Emcn", "Flt1", "Slc5a7",
  277. "Acly", "Slc1a2", "Gfap", "Agt", "Aqp4", "Plp1"
  278. , "Rbfox3")
  280. DotPlot(adata, features = features, group.by = "cell_type") + RotatedAxis()
  282. VlnPlot(adata, "Mog", layer = "scale.data", pt.size = 0)
  284. RidgePlot(adata, features = c("Gad1", "Gad2"), layer = "scale.data", ncol = 2)
  286. RidgePlot(adata, features = c("Rbfox3"), layer = "scale.data", ncol = 1)
  288. RidgePlot(adata, features = c("Pdgfra", "Pcdh15"), layer = "scale.data", ncol = 1)
  290. RidgePlot(adata, features = c("Csf1r", "Cx3cr1", "P2ry12"), layer = "scale.data", ncol = 1)
  292. RidgePlot(adata, features = c("Slc17a7", "Slc17a6"), ncol = 1)
  294. RidgePlot(adata, features = c("Emcn"), layer = "scale.data", ncol = 1)
  296. RidgePlot(adata, features = c("Slc5a7"), layer = "scale.data", ncol = 1)
  298. RidgePlot(adata, features = c("Slc1a2"), layer = "scale.data", ncol = 1)
  300. saveRDS(adata, "high_annotated.rds")
  303. ############################ decontX ###############################
  305. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  307. library(Seurat)
  308. library(celda)
  309. library(qs)
  311. adata = readRDS('high_annotated.rds')
  313. gc()
  315. DefaultAssay(adata) = 'RNA'
  316. adata
  317. gc()
  319. adata[['RNA']] = JoinLayers(adata[['RNA']])
  320. gc()
  322. counts = GetAssayData(adata, slot = 'counts')
  324. gc()
  326. sce = SingleCellExperiment(list(counts = counts))
  329. gc()
  330. sce = decontX(sce)
  332. metadata(sce)$decontX$estimates$all_cells$delta
  334. #öka delta här
  336. sce.delta = decontX(sce, delta = c(3.55, 20), estimateDelta = F)
  338. library(scater)
  340. plot(sce$decontX_contamination, sce.delta$decontX_contamination,
  341. xlab = "DecontX estimated priors",
  342. ylab = "Setting priors to estimate higher contamination")
  343. abline(0, 1, col = "red", lwd = 2)
  345. umap = reducedDim(sce.delta, "decontX_UMAP")
  347. plotDimReduceCluster(sce.delta$decontX_clusters,
  348. dim1 = umap[, 1], dim2 = umap[, 2])
  350. plotDecontXContamination(sce)
  351. plotDecontXContamination(sce.delta)
  353. differences_sce = rowMeans(counts) - rowMeans(decontXcounts(sce))
  356. gc()
  358. top_genes_sce = sort(abs(differences_sce), decreasing = T)[1:20]
  359. top_genes_sce
  361. differences_sce_delta = rowMeans(counts) - rowMeans(decontXcounts(sce.delta))
  364. gc()
  366. top_genes_sce_delta = sort(abs(differences_sce_delta), decreasing = T)[1:20]
  367. top_genes_sce_delta
  369. gc()
  371. sum(counts(sce) < decontXcounts(sce))
  373. celltypemappings = list(Oligodendrocytes = 2
  374. )
  376. plotDecontXMarkerExpression(sce,
  377. markers = 'Kcnip4',
  378. groupClusters = celltypemappings
  379. )
  381. adata[["decontXcounts"]] = CreateAssayObject(counts = decontXcounts(sce.delta))
  384. gc()
  387. VlnPlot(adata, 'Gad1', group.by = 'cell_type', pt.size = 0, assay = 'decontXcounts') + NoLegend()
  388. VlnPlot(adata, 'Gad1', group.by = 'cell_type', pt.size = 0, assay = 'RNA') + NoLegend()
  390. rm(counts)
  392. rna_matrix = GetAssayData(adata, assay = "RNA", layer = 'counts')
  393. gc()
  394. decontx_matrix = GetAssayData(adata, assay = "decontXcounts", layer = "counts")
  396. common_genes = intersect(rownames(rna_matrix), rownames(decontx_matrix))
  397. rna_matrix = rna_matrix[common_genes,]
  398. decontx_matrix = decontx_matrix[common_genes,]
  399. gc()
  400. differences = rowMeans(rna_matrix) - rowMeans(decontx_matrix)
  402. gc()
  404. top_genes = names(sort(abs(differences), decreasing = T))[1:20]
  405. top_genes
  407. VlnPlot(adata, group.by = 'cell_type', 'Nkain2', pt.size = 0, assay = 'decontXcounts') + NoLegend()
  408. VlnPlot(adata, group.by = 'cell_type', 'Nkain2', pt.size = 0, assay = 'RNA') + NoLegend()
  410. genes = c('Kcnip4', 'Nkain2')
  413. for (gene in top_genes) {
  414. DefaultAssay(adata) = 'RNA'
  415. FeaturePlot(adata, gene, reduction = 'euclidean_umap', max.cutoff = 1) + NoLegend()
  417. ggsave(paste0('males/decontx/gabaergic/genes/',gene,'_before_decont.jpg'))
  419. DefaultAssay(adata) = 'decontXcounts'
  420. FeaturePlot(adata, gene, reduction = 'euclidean_umap', max.cutoff = 1) + NoLegend()
  421. ggsave(paste0('males/decontx/gabaergic/genes/',gene,'_after_decont.jpg'))
  425. }
  430. qsave(adata, 'high_annotated_males_decontx.qs')
  432. #------------------------------------subclustering GABA------------------------------
  434. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  436. library(Seurat)
  438. adata = readRDS("high_annotated_males_decontx.qs")
  439. gc()
  442. adata = subset(adata, subset = cell_type == "Gabaergic")
  443. gc()
  445. DefaultAssay(adata) = "decontXcounts"
  449. VlnPlot(adata, c('percent.mt',
  450. 'percent.ribo',
  451. 'nCount_decontXcounts',
  452. 'nFeature_decontXcounts'),
  453. group.by = 'SampleID',
  454. pt.size = 0,
  455. ncol = 2)
  459. adata = subset(adata, subset = nFeature_decontXcounts > 1500
  460. & nCount_decontXcounts < 20000
  461. & nCount_decontXcounts > 500
  462. & percent.ribo < 1
  463. & percent.mt < 3
  464. )
  466. VlnPlot(adata, c('percent.mt',
  467. 'percent.ribo',
  468. 'nCount_decontXcounts',
  469. 'nFeature_decontXcounts'),
  470. group.by = 'SampleID',
  471. pt.size = 0,
  472. ncol = 2)
  474. p1 = FeatureScatter(adata, 'nCount_decontXcounts', 'percent.mt')
  475. p2 = FeatureScatter(adata, 'nCount_decontXcounts', 'nFeature_decontXcounts')
  476. p3 = FeatureScatter(adata, 'nCount_decontXcounts', 'percent.ribo')
  478. p1+p2+p3
  479. gc()
  481. adata[["decontXcounts"]] <- split(adata[["decontXcounts"]], f = adata$SampleID) #split the counts into sample-layers
  483. DefaultAssay(adata) = 'decontXcounts'
  484. adata
  485. gc()
  486. adata <- SCTransform(adata, vars.to.regress = c("percent.mt", "percent.ribo"),
  487. assay = 'decontXcounts')
  488. gc()
  489. adata = RunPCA(adata)
  490. gc()
  491. ElbowPlot(adata,ndims = 50)
  493. adata <- FindNeighbors(adata, dims = 1:50, reduction = "pca")
  494. gc()
  498. adata <- FindClusters(adata, resolution = 0.1, cluster.name = "gaba_subclusters_0.1")
  499. gc()
  500. adata <- FindClusters(adata, resolution = 0.3, cluster.name = "gaba_subclusters_0.3")
  501. gc()
  502. adata <- FindClusters(adata, resolution = 0.5, cluster.name = "gaba_subclusters_0.5")
  503. gc()
  504. adata <- FindClusters(adata, resolution = 0.6, cluster.name = "gaba_subclusters_0.6")
  505. gc()
  506. adata <- FindClusters(adata, resolution = 0.7, cluster.name = "gaba_subclusters_0.7")
  507. gc()
  508. adata <- FindClusters(adata, resolution = 0.8, cluster.name = "gaba_subclusters_0.8")
  509. gc()
  511. adata = RunUMAP(adata, dims = 1:50, metric = 'euclidean', reduction.name = 'euclidean_umap_gaba')
  512. DimPlot(adata, reduction = "euclidean_umap_gaba", group.by = "gaba_subclusters_0.5", label = TRUE)
  514. adata
  517. gc()
  519. for (sample in 1:length(unique(adata$SampleID))) {
  521. print(sample)
  522. slot(adata@assays$[email hidden][[sample]], name="umi.assay") = "decontXcounts"
  523. }
  525. gc()
  526. adata = PrepSCTFindMarkers(adata)
  527. gc()
  529. Idents(adata) = adata$gaba_subclusters_0.5
  531. markers = FindAllMarkers(adata, only.pos = TRUE)
  533. write.csv(markers, 'gaba_males_subclustered_markers.csv')
  536. gc()
  538. DimPlot(adata, reduction = "euclidean_umap_gaba", group.by = "gaba_subclusters_0.5", label = TRUE)
  540. current.cluster_ids = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
  541. 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
  542. 32, 33, 34, 35, 36, 37, 38, 39)
  544. VlnPlot(adata, 'Sst', pt.size = 0, assay = 'decontXcounts', layer = 'data') + NoLegend()
  545. FeaturePlot(adata, c('Gad2', 'Slc17a6'), reduction = 'euclidean_umap_gaba', blend = T,
  546. max.cutoff = 1, pt.size = 2)
  547. FeaturePlot(adata, 'Sst', reduction = 'euclidean_umap_gaba')
  550. new.cluster_ids = c(
  551. 'CeL_Prkcd',#0
  552. 'CeM_Xylt1',#1
  553. 'Intercalated_Chst9',#2
  554. 'Drd2_Htr2c',#3
  555. 'Ebf1_1',#4
  556. 'Drd2_Calcrl',#5
  557. 'Satb1_Col5a2',#6
  558. 'Dlk1',#7
  559. 'Pallidum',#8
  560. 'Ndst3_Dach2',#9
  561. 'Nxph1_Tcf4',#10
  562. 'Npffr1_Dock11',#11
  563. 'Ebf1_2',#12
  564. 'Nos1_Plagl1_Sst',#13
  565. 'Egfr_Nxph1',#14
  566. 'Intercalated_Nos1',#15
  567. 'Prlr_tac3',#16
  568. 'Vipr2_Pdyn',#17
  569. 'Slc17a6',#18
  570. 'Nxph1_Nxph2',#19
  571. 'Ebf1_2',#20
  572. 'Intercalated_Glp1r',#21
  573. 'Tafa1_Slc17a8',#22
  574. 'Prlr_Esr1',#23
  575. 'Stac_Ndst4',#24
  576. 'Intercalated_Chrm2',#25
  577. 'Ebf1_3',#26
  578. 'Slc17a6',#27
  579. 'Drd2_Rai14',#28
  580. 'Mast4',#29
  581. 'Intercalated_Dach1',#30
  582. 'Lgr6',#31
  583. 'Vip',#32
  584. 'Zbtb7c',#33
  585. 'Sox6_Gpc5',#34
  586. 'Sox6_Ano4',#35
  587. 'Fign_Npsr1',#36
  588. 'Grm3',#37
  589. 'Mast4_Zfp536', #38
  590. 'Tf_Car2'#39
  592. )
  595. library(plyr)
  596. adata$cell_subtype = plyr::mapvalues(x = adata$gaba_subclusters_0.5, from = current.cluster_ids, to = new.cluster_ids)
  599. dir.create('males/subclustering')
  600. dir.create('males/subclustering/gabaergic/')
  601. p = DimPlot(adata, reduction = "euclidean_umap_gaba", group.by = "cell_subtype", label = TRUE,
  602. repel = T, raster = F, label.size = 6) + NoLegend()
  603. pdf('males/subclustering/gabaergic/umap_cellsubtypes.pdf',
  604. height = 15, width = 15)
  605. print(p)
  606. dev.off()
  608. qsave(adata, 'gabaergic_males_subclustered.qs')
  610. adata = qread('gabaergic_males_subclustered.qs')
  612. markers = FindAllMarkers(adata, only.pos = T)
  614. markers$cell_type = markers$cluster
  616. gc()
  617. DefaultAssay(adata) = 'SCT'
  618. #adata = NormalizeData(adata)
  619. #adata = ScaleData(adata)
  621. library(dplyr)
  623. markers %>%
  624. group_by(cell_type) %>%
  625. filter(avg_log2FC > 1) %>%
  626. slice_head(n = 20) %>%
  627. ungroup -> top10
  629. gc()
  631. Idents(adata) = 'cell_subtype'
  633. p = DoHeatmap(subset(adata, downsample = 500), features = top10$gene) + NoLegend()
  635. pdf('males/subclustering/gabaergic/heatmap_cellsubtypes.pdf',
  636. height = 10, width = 20)
  637. print(p)
  638. dev.off()
  640. DefaultAssay(adata) = 'decontXcounts'
  641. adata = NormalizeData(adata)
  642. adata = ScaleData(adata)
  644. gc()
  646. features = c('Prkcd', 'Nr2f2', 'Isl1', 'Crh',
  647. 'Calcrl', 'Sst','Oxtr', 'Nts', 'Tac3',
  648. 'Tacr1', 'Cck', 'Cartpt', 'Cxcl14',
  649. 'Htr2a', 'Htr2c', 'Drd1', 'Drd2',
  650. 'Drd3', 'Cnr1',
  651. 'Oprm1', 'Oprk1', 'Oprd1',
  652. 'Gria1', 'Gria2', 'Gria3', 'Gria4',
  653. 'Grik1', 'Grik2', 'Grik3', 'Grik4',
  654. 'Grik5', 'Grin1', 'Grin2a', 'Grin2b',
  655. 'Grin2c', 'Grin2d', 'Grin3a', 'Grin3b',
  656. 'Grm1', 'Grm5', 'Grm2', 'Grm3',
  657. 'Grm4', 'Grm6', 'Grm7', 'Grm8',
  658. 'Adra1a', 'Adra1b', 'Adra1c',
  659. 'Adra2a', 'Adra2b', 'Adra2c',
  660. 'Adrb1', 'Adrb2', 'Adrb3', 'Chrna1',
  661. 'Chrnb1', 'Chrnd', 'Chrne',
  662. 'Chrm1', 'Chrm2', 'Chrm3', 'Chrm4',
  663. 'Adora1','Adora1a', 'Adora2a', 'Adora3',
  664. 'Il13ra1', 'Il13ra2', 'Il1r1', 'Il6r',
  665. 'Tnfrsf1a', 'Tnfrsf1b'
  666. )
  668. celltype_counts = table(adata$cell_subtype)
  669. celltypes_to_keep = names(celltype_counts[celltype_counts >= 500])
  672. p = DoHeatmap(subset(subset(adata, cells = WhichCells(adata, expression = cell_subtype %in% celltypes_to_keep)), downsample = 500), features = features) + NoLegend()
  674. pdf('males/subclustering/gabaergic/heatmap_cellsubtypes_interesting_stuff.pdf',
  675. height = 15, width = 25)
  676. print(p)
  677. dev.off()
  681. FeaturePlot(adata, "Pde1c", reduction = "euclidean_umap_gaba")
  684. RidgePlot(adata, features = "percent.ribo")
  685. RidgePlot(adata, features = "percent.mt")
  686. RidgePlot(adata, features = "nCount_RNA")
  687. RidgePlot(adata, features = "nFeature_RNA")
  688. VlnPlot(adata, features = "percent.mt", group.by = "SampleID", pt.size = 0)
  692. # hdWGCNA
  694. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  698. library(Seurat)
  700. library(tidyverse)
  701. library(cowplot)
  702. library(patchwork)
  704. library(WGCNA)
  705. library(hdWGCNA)
  706. library(writexl)
  707. library(readxl)
  708. library(qs)
  711. # using the cowplot theme for ggplot
  712. theme_set(theme_cowplot())
  714. # set random seed for reproducibility
  715. set.seed(12345)
  717. # optionally enable multithreading
  718. allowWGCNAThreads(nThreads = 7)
  721. #gabaergic cells only
  723. adata = qread('gabaergic_males_subclustered.qs')
  725. describe(adata$nCount_decontXcounts)
  726. describe(adata$nFeature_decontXcounts)
  727. adata
  728. describe(adata$cell_subtype)
  729. unique(adata$cell_subtype)
  732. gc()
  733. adata
  734. Idents(adata) = 'cell_subtype'
  735. markers = FindAllMarkers(adata, only.pos = T)
  737. write_xlsx(markers, 'Fgfr1_paper/gabaergic_markers.xlsx')
  739. markers$cell_type = markers$cluster
  743. Idents(adata) = 'cell_type'
  745. DefaultAssay(adata) = 'decontXcounts'
  751. gc()
  752. adata
  754. adata[['decontXcounts']] = JoinLayers(adata[['decontXcounts']])
  757. adata = NormalizeData(adata)
  759. adata = SetupForWGCNA(
  760. adata,
  761. gene_select = "fraction",
  762. fraction = 0.05,
  763. wgcna_name = 'gaba_males_consensus'
  764. )
  766. gc()
  768. adata
  770. adata = MetacellsByGroups(
  771. adata,
  772. group.by = c('cell_type', 'SampleID', 'Groups'),
  773. ident.group = 'cell_type',
  774. k = 25,
  775. max_shared = 12,
  776. min_cells = 50,
  777. reduction = 'pca'
  778. )
  780. gc()
  781. adata = NormalizeMetacells(adata)
  783. gc()
  784. adata = SetMultiExpr(
  785. adata,
  786. group_name = 'Gabaergic',
  787. group.by = 'cell_type',
  788. multi.group.by = "Groups"
  789. )
  792. adata = TestSoftPowersConsensus(adata)
  794. plot_list = PlotSoftPowers(adata)
  797. consensus_groups <- unique(adata$Groups)
  798. p_list <- lapply(1:length(consensus_groups), function(i){
  799. cur_group <- consensus_groups[[i]]
  800. plot_list[[i]][[1]] + ggtitle(paste0('Group: ', cur_group)) + theme(plot.title=element_text(hjust=0.5))
  801. })
  803. dir.create('males/wgcna/consensus')
  804. dir.create('males/wgcna/consensus/gabaergic')
  805. dir.create('males/wgcna/consensus/gabaergic/plots/')
  807. pdf('males/wgcna/consensus/gabaergic/plots/soft_powers.pdf')
  808. wrap_plots(p_list, ncol=2)
  809. dev.off()
  811. gc()
  813. adata = ConstructNetwork(
  814. adata,
  815. soft_power = c(9,8,8),
  816. consensus = TRUE,
  817. tom_name = 'gaba_male_consensus'
  818. )
  820. pdf('males/wgcna/consensus/gabaergic/plots/dendrogram.pdf')
  821. PlotDendrogram(adata, main = 'Groups consensus dendrogram')
  822. dev.off()
  824. adata = ScaleData(adata)
  826. gc()
  828. adata = ModuleEigengenes(
  829. adata,
  830. group.by.vars = "SampleID"
  831. )
  834. hMEs = GetMEs(adata)
  835. MEs = GetMEs(adata, harmonized=FALSE)
  837. adata = ModuleConnectivity(
  838. adata,
  839. group.by = 'cell_type', group_name = 'Gabaergic'
  840. )
  842. adata = ResetModuleNames(
  843. adata,
  844. new_name = "Gaba-M"
  845. )
  847. gc()
  849. p = PlotKMEs(adata, ncol = 6)
  851. pdf('males/wgcna/consensus/gabaergic/plots/kMEs.pdf',
  852. width = 10, height = 10)
  853. p
  854. dev.off()
  856. modules = GetModules(adata) %>% subset(module != 'grey')
  858. head(modules[,1:6])
  860. write_xlsx(modules, 'males/wgcna/consensus/gabaergic/consensus_modules.xlsx')
  862. qsave(adata, 'males/wgcna/consensus/gabaergic/wgcna_consensus_gaba_males.qs')
  864. ############################ network visualization
  866. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  868. library(Seurat)
  869. library(tidyverse)
  870. library(cowplot)
  871. library(patchwork)
  872. library(WGCNA)
  873. library(hdWGCNA)
  874. library(igraph)
  875. library(qs)
  876. theme_set(theme_cowplot())
  877. set.seed(12345)
  878. gc()
  879. adata = qread('males/wgcna/consensus/gabaergic/wgcna_consensus_gaba_males.qs')
  885. hMEs <- GetMEs(adata)
  888. modules = GetModules(adata)
  889. mods = levels(modules$module); mods = mods[mods != 'grey']
  891. [email hidden] = cbind([email hidden], hMEs)
  893. ################### Make files fo rGEO ##################################
  895. metadata = [email hidden]
  897. write.csv(metadata, 'males/wgcna/consensus/gabaergic/metadata.csv')
  899. dir.create('Fgfr1_paper/for_geo')
  901. gc()
  903. library(Matrix)
  905. raw_counts = adata@assays$RNA@layers$counts
  906. write.table(colnames(raw_counts), "Fgfr1_paper/for_geo/barcodes.tsv")
  907. write.table(rownames(raw_counts), "Fgfr1_paper/for_geo/genes.tsv")
  909. gc()
  910. writeMM(raw_counts, file = 'Fgfr1_paper/for_geo/raw_counts.mtx')
  912. rm(raw_counts)
  913. gc()
  915. decont_counts = adata@assays$decontXcounts@layers$counts
  917. writeMM(decont_counts, file = 'Fgfr1_paper/for_geo/decontXcounts_counts.mtx')
  919. rm(decont_counts)
  920. gc()
  922. decont_SCT_data = adata@assays$SCT@data
  923. write.table(rownames(decont_SCT_data), "Fgfr1_paper/for_geo/genes_SCT.tsv")
  926. writeMM(decont_SCT_data, file = 'Fgfr1_paper/for_geo/decontXcounts_SCT_data.mtx')
  930. p = DotPlot(adata, features = mods ,group.by = 'cell_subtype')
  932. p = p +
  933. RotatedAxis() +
  934. scale_color_gradient2(high = 'red', mid = 'grey95', low = 'blue') +
  935. coord_flip() +
  936. theme(
  937. axis.text.x = element_text(size = 16, face = "bold"), # x-axis labels
  938. axis.text.y = element_text(size = 16, face = "bold"), # y-axis labels
  939. axis.title.x = element_text(size = 2, face = "bold"), # x-axis title
  940. axis.title.y = element_text(size = 2, face = "bold") # y-axis title
  941. )
  944. p
  946. pdf('males/wgcna/consensus/gabaergic/plots/modules_dotplot.pdf',
  947. height = 10, width = 20)
  948. print(p)
  949. dev.off()
  951. PlotKMEs_1 <- function(
  952. seurat_obj,
  953. n_hubs=10,
  954. text_size=2,
  955. ncol = 5,
  956. plot_widths = c(3,2),
  957. wgcna_name = NULL
  958. ){
  960. if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
  962. modules <- GetModules(seurat_obj, wgcna_name) %>% subset(module != 'grey')
  963. mods <- levels(modules$module); mods <- mods[mods != 'grey']
  964. mod_colors <- modules %>% subset(module %in% mods) %>%
  965. dplyr::select(c(module, color)) %>%
  966. dplyr::distinct()
  968. # get hub genes:
  969. hub_df <- GetHubGenes(seurat_obj, n_hubs=n_hubs, wgcna_name=wgcna_name)
  971. plot_list <- lapply(mods, function(x){
  972. cur_color <- subset(mod_colors, module == x) %>% .$color
  973. cur_df <- subset(hub_df, module == x)
  974. top_genes <- cur_df %>% dplyr::top_n(n_hubs, wt=kME) %>% .$gene_name
  975. p <- cur_df %>% ggplot(aes(x = reorder(gene_name, kME), y = kME)) +
  976. geom_bar(stat='identity', width=1, color = cur_color, fill=cur_color) +
  977. ggtitle(x) +
  978. ylim(0, 1) +
  979. #xlab(paste0('kME_', x)) +
  980. theme(
  981. axis.ticks.x = element_blank(),
  982. axis.text.x = element_blank(),
  983. plot.title = element_text(hjust=0.5, size = 36, face = 'bold'),
  984. axis.title.x = element_blank(),
  985. axis.line.x = element_blank()
  986. )
  987. p_anno <- ggplot() + annotate(
  988. "label",
  989. x = 0,
  990. y = 0,
  991. label = paste0(top_genes, collapse="\n"),
  992. size=text_size,
  993. fontface = 'bold',
  994. label.size=0
  995. ) + theme_void()
  996. patch <- p + p_anno + plot_layout(widths=plot_widths)
  997. patch
  998. })
  1000. wrap_plots(plot_list, ncol=ncol)
  1002. }
  1006. p = PlotKMEs_1(adata, ncol = 3, text_size = 4)
  1007. p
  1009. pdf('males/wgcna/consensus/gabaergic/plots/kMEs.pdf',
  1010. width = 20, height = 15)
  1011. p
  1012. dev.off()
  1014. plot_list = ModuleFeaturePlot(
  1015. adata,
  1016. features = 'hMEs',
  1017. order=TRUE,
  1018. reduction = 'euclidean_umap_gaba'
  1019. )
  1022. pdf('males/wgcna/consensus/gabaergic/plots/modules_umap.pdf',
  1023. height = 10, width = 10)
  1024. wrap_plots(plot_list, ncol = 6)
  1025. dev.off()
  1027. p = ModuleRadarPlot(
  1028. adata,
  1029. group.by = 'cell_subtype', axis.label.size = 3, grid.label.size = 3
  1030. )
  1033. pdf('males/wgcna/consensus/gabaergic/plots/radar_plot_cell_subtype.pdf',
  1034. height = 10, width = 15)
  1035. p
  1036. dev.off()
  1039. ModuleNetworkPlot(
  1040. adata,
  1041. outdir = 'males/wgcna/consensus/gabaergic/plots/mandala_plots'
  1042. )
  1044. library(readxl)
  1046. ################## DME analysis ##################
  1050. # celltypes
  1051. dir.create('males/wgcna/consensus/gabaergic/DMEs')
  1053. #res vs sen
  1055. DMEs <- data.frame()
  1057. for (cell_type in unique(adata$cell_subtype)) {
  1058. print(cell_type)
  1059. # group1 = [email hidden] %>% subset(cell_subtype == cell_type & Groups == 'Resistant') %>% rownames
  1060. # group2 = [email hidden] %>% subset(cell_subtype == cell_type & Groups == 'Sensitive') %>% rownames
  1062. group1 <- rownames(
  1063. [email hidden][
  1064. as.character([email hidden]$cell_subtype) == cell_type &
  1065. [email hidden]$Groups == "Resistant", ]
  1066. )
  1068. group2 <- rownames(
  1069. [email hidden][
  1070. as.character([email hidden]$cell_subtype) == cell_type &
  1071. [email hidden]$Groups == "Sensitive", ]
  1072. )
  1074. dir.create('males/wgcna/consensus/gabaergic/DMEs/res_vs_sen')
  1075. dir.create(paste0('males/wgcna/consensus/gabaergic/DMEs/res_vs_sen/',cell_type))
  1077. cur_DMEs <- FindDMEs(
  1078. adata,
  1079. barcodes1 = group1,
  1080. barcodes2 = group2,
  1081. test.use='wilcox',
  1082. wgcna_name='gaba_males_consensus'
  1083. )
  1085. cur_DMEs$celltype = cell_type
  1087. DMEs = rbind(DMEs, cur_DMEs)
  1089. p = PlotDMEsLollipop(
  1090. adata,
  1091. cur_DMEs,
  1092. wgcna_name='gaba_males_consensus',
  1093. pvalue = "p_val_adj"
  1094. )
  1096. pdf(paste0('males/wgcna/consensus/gabaergic/DMEs/res_vs_sen/',cell_type,'/lollipop.pdf'))
  1097. print(p)
  1098. dev.off()
  1103. }
  1106. modules <- GetModules(adata)
  1107. mods <- levels(modules$module); mods <- mods[mods != 'grey']
  1109. plot_df <- DMEs
  1111. plot_df$module <- factor(as.character(plot_df$module), levels=mods)
  1113. maxval <- 0.5; minval <- -0.5
  1114. plot_df$avg_log2FC <- ifelse(plot_df$avg_log2FC > maxval, maxval, plot_df$avg_log2FC)
  1115. plot_df$avg_log2FC <- ifelse(plot_df$avg_log2FC < minval, minval, plot_df$avg_log2FC)
  1117. plot_df$Significance <- gtools::stars.pval(plot_df$p_val_adj)
  1119. plot_df$textcolor <- ifelse(plot_df$avg_log2FC > 0.2, 'black', 'white')
  1121. p <- plot_df %>%
  1122. ggplot(aes(y=celltype, x=module, fill=avg_log2FC)) +
  1123. geom_tile()
  1125. p <- p +
  1126. geom_text(label=plot_df$Significance, color=plot_df$textcolor)
  1128. p <- p +
  1129. scale_fill_gradient2(low='purple', mid='black', high='yellow') +
  1130. RotatedAxis() +
  1131. theme(
  1132. panel.border = element_rect(fill=NA, color='black', size=1),
  1133. axis.line.x = element_blank(),
  1134. axis.line.y = element_blank(),
  1135. plot.margin=margin(0,0,0,0)
  1136. ) + xlab('') + ylab('') +
  1137. coord_equal() +
  1138. ggtitle("Resistant vs Sensitive")
  1141. pdf('males/wgcna/consensus/gabaergic/DMEs/res_vs_sen/dmes_heatmap_celltype.pdf',
  1142. height = 10, width = 10)
  1143. p
  1144. dev.off()
  1146. library(writexl)
  1148. write_xlsx(DMEs, 'males/wgcna/consensus/gabaergic/DMEs/res_vs_sen/dmes_celltype.xlsx')
  1151. #res vs yok
  1153. DMEs <- data.frame()
  1155. for (cell_type in unique(adata$cell_subtype)) {
  1156. print(cell_type)
  1157. # group1 = [email hidden] %>% subset(cell_subtype == cell_type & Groups == 'Resistant') %>% rownames
  1158. # group2 = [email hidden] %>% subset(cell_subtype == cell_type & Groups == 'Sensitive') %>% rownames
  1160. group1 <- rownames(
  1161. [email hidden][
  1162. as.character([email hidden]$cell_subtype) == cell_type &
  1163. [email hidden]$Groups == "Resistant", ]
  1164. )
  1166. group2 <- rownames(
  1167. [email hidden][
  1168. as.character([email hidden]$cell_subtype) == cell_type &
  1169. [email hidden]$Groups == "Yoked", ]
  1170. )
  1172. dir.create('males/wgcna/consensus/gabaergic/DMEs/res_vs_yok')
  1173. dir.create(paste0('males/wgcna/consensus/gabaergic/DMEs/res_vs_yok/',cell_type))
  1175. cur_DMEs <- FindDMEs(
  1176. adata,
  1177. barcodes1 = group1,
  1178. barcodes2 = group2,
  1179. test.use='wilcox',
  1180. wgcna_name='gaba_males_consensus'
  1181. )
  1183. cur_DMEs$celltype = cell_type
  1185. DMEs = rbind(DMEs, cur_DMEs)
  1187. p = PlotDMEsLollipop(
  1188. adata,
  1189. cur_DMEs,
  1190. wgcna_name='gaba_males_consensus',
  1191. pvalue = "p_val_adj"
  1192. )
  1194. pdf(paste0('males/wgcna/consensus/gabaergic/DMEs/res_vs_yok/',cell_type,'/lollipop.pdf'))
  1195. print(p)
  1196. dev.off()
  1201. }
  1204. modules <- GetModules(adata)
  1205. mods <- levels(modules$module); mods <- mods[mods != 'grey']
  1207. plot_df <- DMEs
  1209. plot_df$module <- factor(as.character(plot_df$module), levels=mods)
  1211. maxval <- 0.5; minval <- -0.5
  1212. plot_df$avg_log2FC <- ifelse(plot_df$avg_log2FC > maxval, maxval, plot_df$avg_log2FC)
  1213. plot_df$avg_log2FC <- ifelse(plot_df$avg_log2FC < minval, minval, plot_df$avg_log2FC)
  1215. plot_df$Significance <- gtools::stars.pval(plot_df$p_val_adj)
  1217. plot_df$textcolor <- ifelse(plot_df$avg_log2FC > 0.2, 'black', 'white')
  1219. p <- plot_df %>%
  1220. ggplot(aes(y=celltype, x=module, fill=avg_log2FC)) +
  1221. geom_tile()
  1223. p <- p +
  1224. geom_text(label=plot_df$Significance, color=plot_df$textcolor)
  1226. p <- p +
  1227. scale_fill_gradient2(low='purple', mid='black', high='yellow') +
  1228. RotatedAxis() +
  1229. theme(
  1230. panel.border = element_rect(fill=NA, color='black', size=1),
  1231. axis.line.x = element_blank(),
  1232. axis.line.y = element_blank(),
  1233. plot.margin=margin(0,0,0,0)
  1234. ) + xlab('') + ylab('') +
  1235. coord_equal() +
  1236. ggtitle("Resistant vs Yoked")
  1239. pdf('males/wgcna/consensus/gabaergic/DMEs/res_vs_yok/dmes_heatmap_celltype.pdf',
  1240. height = 10, width = 10)
  1241. p
  1242. dev.off()
  1244. library(writexl)
  1246. write_xlsx(DMEs, 'males/wgcna/consensus/gabaergic/DMEs/res_vs_yok/dmes_celltype.xlsx')
  1249. #Enrichment
  1253. setwd("D:/Leon/snRNAseq/heilig_data_transfer/R/high_depth")
  1255. library(Seurat)
  1256. library(tidyverse)
  1257. library(cowplot)
  1258. library(patchwork)
  1259. library(WGCNA)
  1260. library(hdWGCNA)
  1261. library(igraph)
  1262. library(qs)
  1263. theme_set(theme_cowplot())
  1264. set.seed(12345)
  1265. gc()
  1266. adata = qread('males/wgcna/consensus/gabaergic/wgcna_consensus_gaba_males.qs')
  1269. library(enrichR)
  1270. library(GeneOverlap)
  1272. dbs = c('GO_Biological_Process_2025',
  1273. 'GO_Cellular_Component_2025',
  1274. 'GO_Molecular_Function_2025')
  1276. adata = RunEnrichr(
  1277. adata,
  1278. dbs=dbs,
  1279. max_genes = 100
  1280. )
  1282. enrich_df = GetEnrichrTable(adata)
  1284. dir.create('males/wgcna/consensus/gabaergic/enrichment')
  1286. write_xlsx(enrich_df, 'males/wgcna/consensus/gabaergic/enrichment/modules_enrichment.xlsx')
  1288. EnrichrBarPlot_1 <- function(
  1289. seurat_obj,
  1290. outdir = "enrichr_plots",
  1291. n_terms = 25,
  1292. p_cutoff = 0.05,
  1293. p_adj = TRUE,
  1294. plot_size = c(6,15),
  1295. logscale=FALSE,
  1296. plot_bar_color=NULL,
  1297. plot_text_color=NULL,
  1298. wgcna_name=NULL
  1299. ){
  1301. # get data from active assay if wgcna_name is not given
  1302. if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
  1304. # get modules:
  1305. modules <- GetModules(seurat_obj, wgcna_name)
  1306. mods <- levels(modules$module)
  1307. mods <- mods[mods != 'grey']
  1309. # get Enrichr table
  1310. enrichr_df <- GetEnrichrTable(seurat_obj, wgcna_name)
  1311. dbs <- as.character(unique(enrichr_df$db))
  1313. # subset based on significance level:
  1314. if(p_adj){
  1315. enrichr_df <- subset(enrichr_df, Adjusted.P.value <= p_cutoff)
  1316. } else{
  1317. enrichr_df <- subset(enrichr_df, P.value <= p_cutoff)
  1318. }
  1320. # helper function to wrap text
  1321. wrapText <- function(x, len) {
  1322. sapply(x, function(y) paste(strwrap(y, len), collapse = "\n"), USE.NAMES = FALSE)
  1323. }
  1325. # make output dir if it doesn't exist:
  1326. if(!dir.exists(outdir)){dir.create(outdir)}
  1328. # loop through modules:
  1329. for(i in 1:length(mods)){
  1331. cur_mod <- mods[i]
  1332. cur_terms <- subset(enrichr_df, module == cur_mod)
  1333. print(cur_mod)
  1335. # get color for this module:
  1336. cur_color <- modules %>% subset(module == cur_mod) %>% .$color %>% unique %>% as.character
  1337. if(!is.null(plot_bar_color)){
  1338. cur_color <- plot_bar_color
  1339. }
  1341. cur_color <- grDevices::adjustcolor(cur_color, alpha.f = 0.6)
  1343. # skip if there are not any terms for this module:
  1344. if(nrow(cur_terms) == 0){next}
  1345. cur_terms$wrap <- wrapText(cur_terms$Term, 45)
  1347. # plot top n_terms as barplot
  1348. plot_list <- list()
  1349. for(cur_db in dbs){
  1351. plot_df <- subset(cur_terms, db==cur_db) %>%
  1352. slice_max(order_by=Combined.Score, n=n_terms)
  1354. # text color:
  1355. if(is.null(plot_text_color)){
  1356. if(cur_color == 'black'){
  1357. text_color = 'grey'
  1358. } else {
  1359. text_color = 'black'
  1360. }
  1361. } else{
  1362. text_color <- plot_text_color
  1363. }
  1366. # logscale?
  1367. if(logscale){
  1368. plot_df$Combined.Score <- log(plot_df$Combined.Score)
  1369. lab <- 'Enrichment log(combined score)'
  1370. x <- 0.2
  1371. } else{lab <- 'Enrichment (combined score)'; x <- 5}
  1373. # make bar plot:
  1374. plot_list[[cur_db]] <- ggplot(plot_df, aes(x=Combined.Score, y=reorder(wrap, Combined.Score)))+
  1375. geom_bar(stat='identity', position='identity', color='white', fill=cur_color) +
  1376. geom_text(aes(label=wrap), x=x, color=text_color, size=8, hjust='left') +
  1377. scale_x_continuous(expand = c(0, 0), limits = c(0, NA)) +
  1378. xlab(lab) + ylab('') + ggtitle(cur_db) +
  1379. theme(
  1380. panel.grid.major=element_blank(),
  1381. panel.grid.minor=element_blank(),
  1382. legend.title = element_blank(),
  1383. axis.ticks.y=element_blank(),
  1384. axis.text.y=element_blank(),
  1385. plot.title = element_text(hjust = 0.5),
  1386. axis.line.y=element_blank()
  1387. )
  1388. }
  1390. # make pdfs in output dir
  1391. pdf(paste0(outdir, '/', cur_mod, '.pdf'), width=plot_size[1], height=plot_size[2])
  1392. for(plot in plot_list){
  1393. print(plot)
  1394. }
  1395. dev.off()
  1396. }
  1397. }
  1400. EnrichrBarPlot_1(
  1401. adata,
  1402. outdir = 'males/wgcna/consensus/gabaergic/enrichment',
  1403. n_terms = 5,
  1404. plot_size = c(10,10),
  1405. logscale = TRUE, p_cutoff = 0.1
  1406. )
  1408. p = EnrichrDotPlot(
  1409. adata,
  1410. mods = "all", # use all modules (default)
  1411. database = "GO_Biological_Process_2025", # this must match one of the dbs used previously
  1412. n_terms=2, # number of terms per module
  1413. term_size=10, # font size for the terms
  1414. p_adj = FALSE # show the p-val or adjusted p-val?
  1415. ) + scale_color_stepsn(colors=rev(viridis::magma(256)))
  1417. pdf('males/wgcna/consensus/gabaergic/enrichment/enrichment_dotplot_bio_process.pdf',
  1418. height = 15, width = 15)
  1419. print(p)
  1420. dev.off()
  1422. p = EnrichrDotPlot(
  1423. adata,
  1424. mods = "all", # use all modules (default)
  1425. database = "GO_Cellular_Component_2025", # this must match one of the dbs used previously
  1426. n_terms=2, # number of terms per module
  1427. term_size=10, # font size for the terms
  1428. p_adj = FALSE # show the p-val or adjusted p-val?
  1429. ) + scale_color_stepsn(colors=rev(viridis::magma(256)))
  1431. pdf('males/wgcna/consensus/gabaergic/enrichment/enrichment_dotplot_cell_process.pdf',
  1432. height = 15, width = 15)
  1433. print(p)
  1434. dev.off()
  1436. p = EnrichrDotPlot(
  1437. adata,
  1438. mods = "all", # use all modules (default)
  1439. database = "GO_Molecular_Function_2025", # this must match one of the dbs used previously
  1440. n_terms=2, # number of terms per module
  1441. term_size=10, # font size for the terms
  1442. p_adj = FALSE # show the p-val or adjusted p-val?
  1443. ) + scale_color_stepsn(colors=rev(viridis::magma(256)))
  1445. pdf('males/wgcna/consensus/gabaergic/enrichment/enrichment_dotplot_mol_function.pdf',
  1446. height = 15, width = 15)
  1447. print(p)
  1448. dev.off()
  1451. ################### modules in the same direction ##################
  1453. library(readxl)
  1454. res_vs_sen = read_excel('males/wgcna/consensus/gabaergic/DMEs/res_vs_sen/dmes_celltype.xlsx')
  1455. res_vs_yok = read_excel('males/wgcna/consensus/gabaergic/DMEs/res_vs_yok/dmes_celltype.xlsx')
  1457. res_vs_sen$comparison = 'res_vs_sen'
  1458. res_vs_yok$comparison = 'res_vs_yok'
  1460. df1_sig = subset(res_vs_sen, p_val_adj < 0.05)
  1461. df2_sig = subset(res_vs_yok, p_val_adj < 0.05)
  1463. library(dplyr)
  1465. combined_data = merge(df1_sig, df2_sig, c('celltype', 'module'))
  1466. combined_data$avg_log2FC.x <- as.numeric(combined_data$avg_log2FC.x)
  1467. combined_data$avg_log2FC.y <- as.numeric(combined_data$avg_log2FC.y)
  1469. df_filtered <- combined_data %>%
  1470. filter(
  1471. (avg_log2FC.x > 0 & avg_log2FC.y > 0 | avg_log2FC.x < 0 & avg_log2FC.y < 0) &
  1472. abs(avg_log2FC.x) > 0.2 &
  1473. abs(avg_log2FC.y) > 0.2 &
  1474. pct.1.x > 0.5 &
  1475. pct.1.y > 0.5 &
  1476. pct.2.x > 0.5 &
  1477. pct.2.y > 0.5
  1478. )
  1480. write_xlsx(df_filtered, 'Fgfr1_paper/DMEs_for_paper.xlsx')
  1482. long_data <- df_filtered %>%
  1483. pivot_longer(
  1484. cols = c(avg_log2FC.x, avg_log2FC.y),
  1485. names_to = "comparison",
  1486. names_pattern = "avg_log2FC\\.(.)",
  1487. values_to = "logFC"
  1488. ) %>%
  1489. mutate(
  1490. comparison_label = ifelse(comparison == "x", "Res vs Sen", "Res vs Yok"),
  1491. module_label = paste0(module, " (", comparison_label, ")")
  1492. )
  1495. write_xlsx(long_data, 'Fgfr1_paper/DMEs_for_paper.xlsx')
  1497. #################### Distinct #####################
  1499. setwd('outputs/distinct/males/')
  1501. library(Seurat)
  1502. library(SingleCellExperiment)
  1503. library(distinct)
  1504. library(qs)
  1505. library(writexl)
  1506. library(readxl)
  1508. adata = qread('../../../objects/gabaergic_males_subclustered.qs')
  1509. DefaultAssay(adata) = 'decontXcounts'
  1510. adata[['decontXcounts']] = JoinLayers(adata[['decontXcounts']])
  1511. adata
  1512. sce = as.SingleCellExperiment(adata)
  1513. rm(adata)
  1514. sce
  1515. gc()
  1516. sce
  1519. cpm_mat = scater::calculateCPM(counts(sce))
  1521. assay(sce, "cpm") = cpm_mat
  1523. sce = scater::logNormCounts(sce)
  1524. sce
  1526. gc()
  1527. colData(sce)
  1528. gc()
  1529. sce$SampleID = as.factor(sce$SampleID)
  1530. sce$Groups = as.factor(sce$Groups)
  1533. sample_md <- unique(as.data.frame(colData(sce)[, c("SampleID", "Groups")]))
  1534. design <- model.matrix(~Groups, data = sample_md)
  1535. rownames(design) <- sample_md$SampleID
  1536. design
  1539. rownames(design)
  1541. set.seed(1234)
  1545. ###################################
  1548. sce_sub <- sce[, sce$Groups %in% c("Resistant", "Sensitive")]
  1551. unique(sce_sub$Groups)
  1552. sce_sub$Groups = droplevels(sce_sub$Groups)
  1553. unique(sce_sub$Groups)
  1555. set.seed(123)
  1558. sample_md <- unique(as.data.frame(colData(sce_sub)[, c("SampleID", "Groups")]))
  1561. design_sub <- model.matrix(~ Groups, data = sample_md)
  1562. rownames(design_sub) <- sample_md$SampleID
  1565. design_sub
  1566. gc()
  1568. sce_sub
  1570. sum(logcounts(sce_sub)["Sik2", ] > 0)
  1573. res_sen <- distinct_test(
  1574. x = sce_sub,
  1575. name_assays_expression = "cpm",
  1576. name_cluster = "cell_subtype",
  1577. name_sample = "SampleID",
  1578. design = design_sub,
  1579. column_to_test = 2,
  1580. min_non_zero_cells = 2000,
  1581. n_cores = 10,
  1582. P_4 = 20000
  1583. )
  1586. res_sen = log2_FC(res = res_sen,
  1587. x = sce_sub,
  1588. name_assays_expression = 'cpm',
  1589. name_group = "Groups",
  1590. name_cluster = "cell_subtype")
  1594. top = top_results(res_sen,
  1595. significance = 0.05, global = T)
  1597. plot_densities(x = sce_sub,
  1598. gene = 'Cdk10',
  1599. cluster = 'Drd2_Calcrl',
  1600. name_assays_expression = 'cpm',
  1601. name_cluster = 'cell_subtype',
  1602. name_sample = 'SampleID',
  1603. name_group = 'Groups', group_level = F)
  1605. plot_densities(x = sce_sub,
  1606. gene = 'Cdk10',
  1607. cluster = 'Drd2_Calcrl',
  1608. name_assays_expression = 'cpm',
  1609. name_cluster = 'cell_subtype',
  1610. name_sample = 'SampleID',
  1611. name_group = 'Groups', group_level = T)
  1613. write_xlsx(x = top, 'res_sen.xlsx')
  1615. #########################################
  1617. sce_sub <- sce[, sce$Groups %in% c("Resistant", "Yoked")]
  1620. unique(sce_sub$Groups)
  1621. sce_sub$Groups = droplevels(sce_sub$Groups)
  1622. unique(sce_sub$Groups)
  1624. set.seed(123)
  1627. sample_md <- unique(as.data.frame(colData(sce_sub)[, c("SampleID", "Groups")]))
  1630. design_sub <- model.matrix(~ Groups, data = sample_md)
  1631. rownames(design_sub) <- sample_md$SampleID
  1634. design_sub
  1635. gc()
  1637. sce_sub
  1639. sum(logcounts(sce_sub)["Slc17a7", ] > 0)
  1642. res_yok <- distinct_test(
  1643. x = sce_sub,
  1644. name_assays_expression = "cpm",
  1645. name_cluster = "cell_subtype",
  1646. name_sample = "SampleID",
  1647. design = design_sub,
  1648. column_to_test = 2,
  1649. min_non_zero_cells = 2000,
  1650. n_cores = 10,
  1651. P_4 = 20000
  1652. )
  1655. res_yok = log2_FC(res = res_yok,
  1656. x = sce_sub,
  1657. name_assays_expression = 'cpm',
  1658. name_group = "Groups",
  1659. name_cluster = "cell_subtype")
  1663. top = top_results(res_yok,
  1664. significance = 0.05, global = T)
  1666. plot_densities(x = sce_sub,
  1667. gene = 'Cdk10',
  1668. cluster = 'Drd2_Calcrl',
  1669. name_assays_expression = 'cpm',
  1670. name_cluster = 'cell_subtype',
  1671. name_sample = 'SampleID',
  1672. name_group = 'Groups', group_level = F)
  1674. plot_densities(x = sce_sub,
  1675. gene = 'Cdk10',
  1676. cluster = 'Drd2_Calcrl',
  1677. name_assays_expression = 'cpm',
  1678. name_cluster = 'cell_subtype',
  1679. name_sample = 'SampleID',
  1680. name_group = 'Groups', group_level = T)
  1683. write_xlsx(x = top, 'res_yok.xlsx')
  1684. gc()
  1686. ######################### Analysis of results #################################
  1688. res_sen = read_excel('res_sen.xlsx')
  1689. res_yok = read_excel('res_yok.xlsx')
  1692. head(res_sen)
  1693. head(res_yok)
  1696. library(dplyr)
  1698. overlap <- inner_join(
  1699. res_sen %>% select(gene, cluster_id, everything()),
  1700. res_yok %>% select(gene, cluster_id, everything()),
  1701. by = c("gene", "cluster_id"),
  1702. suffix = c("_sen", "_yok")
  1703. )
  1708. overlap_same_direction <- overlap %>%
  1709. filter(
  1710. (`log2FC_Resistant/Sensitive` > 0 & `log2FC_Resistant/Yoked` > 0) |
  1711. (`log2FC_Resistant/Sensitive` < 0 & `log2FC_Resistant/Yoked` < 0)
  1712. )
  1714. overlap_same_direction$average_log2FC = (overlap_same_direction$`log2FC_Resistant/Sensitive`+overlap_same_direction$`log2FC_Resistant/Yoked`)/2
  1715. overlap_same_direction$average_p_adj = (overlap_same_direction$p_adj.glb_sen + overlap_same_direction$p_adj.glb_yok)/2
  1718. write_xlsx(overlap_same_direction, 'DEG_list_males_gabaergic.xlsx')
  1720. for_ipa = overlap_same_direction %>% select(gene, average_log2FC, average_p_adj, cluster_id)
  1722. dir.create('for_ipa')
  1723. celltypes <- unique(df$cluster_id)
  1725. for (ct in celltypes) {
  1726. sub <- for_ipa %>% filter(cluster_id == ct)
  1727. write.csv(sub, paste0("for_ipa/IPA_", ct, ".csv"), row.names = FALSE)
  1728. }
  1733. ################################# Volcanoplots ############################3
  1735. overlap_same_direction = read_xlsx('DEG_list_males_gabaergic.xlsx')
  1738. library(dplyr)
  1739. library(ggpubr)
  1740. library(viridis)
  1741. library(dplyr)
  1742. library(ggplot2)
  1743. library(patchwork)
  1744. library(ggrepel)
  1747. df = overlap_same_direction
  1749. df <- df %>%
  1750. mutate(
  1751. reg_sen = case_when(
  1752. p_adj.glb_sen < 0.05 & `log2FC_Resistant/Sensitive` > 0.1 ~ "Up",
  1753. p_adj.glb_sen < 0.05 & `log2FC_Resistant/Sensitive` < -0.1 ~ "Down",
  1754. TRUE ~ "NS"
  1755. ),
  1756. reg_yok = case_when(
  1757. p_adj.glb_yok < 0.05 & `log2FC_Resistant/Yoked` > 0.1 ~ "Up",
  1758. p_adj.glb_yok < 0.05 & `log2FC_Resistant/Yoked` < -0.1 ~ "Down",
  1759. TRUE ~ "NS"
  1760. )
  1761. )
  1764. make_volcano_pair <- function(df, cluster) {
  1766. df_clust <- df %>% filter(cluster_id == cluster)
  1768. p_sen <- ggplot(df_clust, aes(
  1769. x = `log2FC_Resistant/Sensitive`,
  1770. y = -log10(p_adj.glb_sen)
  1771. )) +
  1772. geom_point(aes(color = reg_sen), size = 1.2, alpha = 0.7) +
  1773. geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") +
  1774. geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  1775. geom_text_repel(
  1776. data = df_clust %>% filter(reg_sen != "NS"),
  1777. aes(label = gene),
  1778. size = 3,
  1779. max.overlaps = 20
  1780. ) +
  1781. scale_color_manual(
  1782. values = c(Up = "red", Down = "blue", NS = "grey70")
  1783. ) +
  1784. labs(
  1785. title = "Resistant vs Sensitive",
  1786. x = "log2FC",
  1787. y = "-log10(adj p-value)"
  1788. ) +
  1789. theme_classic() +
  1790. theme(legend.position = "none")
  1792. p_yok <- ggplot(df_clust, aes(
  1793. x = `log2FC_Resistant/Yoked`,
  1794. y = -log10(p_adj.glb_yok)
  1795. )) +
  1796. geom_point(aes(color = reg_yok), size = 1.2, alpha = 0.7) +
  1797. geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") +
  1798. geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  1799. geom_text_repel(
  1800. data = df_clust %>% filter(reg_yok != "NS"),
  1801. aes(label = gene),
  1802. size = 3,
  1803. max.overlaps = 20
  1804. ) +
  1805. scale_color_manual(
  1806. values = c(Up = "red", Down = "blue", NS = "grey70")
  1807. ) +
  1808. labs(
  1809. title = "Resistant vs Yoked",
  1810. x = "log2FC",
  1811. y = "-log10(adj p-value)"
  1812. ) +
  1813. theme_classic() +
  1814. theme(legend.position = "none")
  1816. (p_sen | p_yok) + plot_annotation(title = cluster)
  1817. }
  1821. clusters <- unique(df$cluster_id)
  1823. volcano_plots <- lapply(clusters, function(cl) {
  1824. make_volcano_pair(df, cl)
  1825. })
  1829. volcano_plots[[1]]
  1831. dir.create('volcano_plots')
  1834. clusters <- unique(df$cluster_id)
  1836. for (cl in clusters) {
  1838. p <- make_volcano_pair(df, cl)
  1840. ggsave(
  1841. filename = paste0("volcano_plots/volcano_", cl, ".pdf"),
  1842. plot = p,
  1843. width = 15,
  1844. height = 10,
  1845. units = "in"
  1846. )
  1847. }
  1851. ################################################################################
  1853. df_plot <- overlap_same_direction %>% count(cluster_id)
  1855. p = ggbarplot(
  1856. df_plot,
  1857. x = "cluster_id",
  1858. y = "n",
  1859. fill = "cluster_id",
  1860. palette = viridis(nrow(df_plot)),
  1861. label = TRUE,
  1862. lab.vjust = -0.5,
  1863. sort.val = "desc", # sort bars by height (largest at top)
  1864. sort.by.groups = FALSE,
  1865. x.text.angle = 45,
  1866. legend = "none" # remove legend
  1867. ) +
  1868. labs(
  1869. title = "Number of DEGs per cluster",
  1870. x = "Cluster ID",
  1871. y = "Gene Count"
  1872. )
  1873. p
  1874. pdf('number_of_degs_per_cluster.pdf')
  1875. print(p)
  1876. dev.off()
  1879. m8 = read_excel('modules/M8_consensus.xlsx', col_names = 'gene')
  1880. m5 = read_excel('modules/M5_consensus.xlsx', col_names = 'gene')
  1881. m4 = read_excel('modules/M4_consensus.xlsx', col_names = 'gene')
  1882. m2 = read_excel('modules/M2_consensus.xlsx', col_names = 'gene')
  1884. degs_in_m8 = overlap_same_direction %>%
  1885. filter(gene %in% m8$gene)
  1887. degs_in_m5 = overlap_same_direction %>%
  1888. filter(gene %in% m5$gene)
  1890. degs_in_m4 = overlap_same_direction %>%
  1891. filter(gene %in% m4$gene)
  1893. degs_in_m2 = overlap_same_direction %>%
  1894. filter(gene %in% m2$gene)
  1899. adata = qread('../../../objects/gabaergic_males_subclustered.qs')
  1902. FeaturePlot(adata, reduction = 'euclidean_umap_gaba', 'Tmcc3')
  1904. p = DotPlot(
  1905. adata,
  1906. c('Pdzrn3', 'Sox5', 'Tmcc3', 'Arhgap6', 'Hdac9', 'Ankfn1', 'Dock10',
  1907. 'Palmd', 'Pde7b', 'Maml2', 'Fgfr1', 'Slc35f3', 'Spata13', 'Epb41l4b',
  1908. 'Kalrn', 'Shank2', 'Sec14l1', 'Camk4', 'Rap1gap', 'Mbnl1', 'Ptpro',
  1909. 'Cnih3', 'Glis3', 'Vat1l', 'Dgkg'),
  1910. group.by = 'cell_subtype'
  1911. )
  1913. p = p + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  1916. dir.create('../../males')
  1917. dir.create('../../males/misc_plots')
  1919. pdf('../../males/misc_plots/dotplot_hubgenes_celltypes.pdf',
  1920. width = 20, height = 20)
  1921. print(p)
  1922. dev.off()
  1925. ################################ Main celltypes ##################################################
  1926. setwd("~/Documents/Projects/Compulsivity/snRNAseq_males_and_females/compulsive_snrnaseq_females_males")
  1928. dir.create('outputs/distinct/males/main_celltypes')
  1929. setwd('outputs/distinct/males/main_celltypes')
  1931. library(Seurat)
  1932. library(SingleCellExperiment)
  1933. library(distinct)
  1934. library(qs)
  1935. library(writexl)
  1936. library(readxl)
  1937. library(dplyr)
  1938. adata = readRDS('../../../../objects/high_annotated_males_decontx.rds')
  1939. DefaultAssay(adata) = 'decontXcounts'
  1940. adata
  1942. adata$cell_type = recode(
  1943. adata$cell_type,
  1944. "Oligendrocytes" = "Oligodendrocytes",
  1945. "Hbb-b1_Hbb-a2_Blood?" = "Blood",
  1946. "Glutaminergic_Slc17a6" = "Glutamatergic_Vglut2",
  1947. "Glutaminergic_Slc17a7" = "Glutamatergic_Vglut1"
  1948. )
  1950. unique(adata$cell_type)
  1951. sce = as.SingleCellExperiment(adata)
  1952. rm(adata)
  1953. sce
  1954. gc()
  1956. cpm_mat = scater::calculateCPM(counts(sce))
  1958. assay(sce, "cpm") = cpm_mat
  1959. gc()
  1960. sce = scater::logNormCounts(sce)
  1961. sce
  1963. gc()
  1964. colData(sce)
  1965. gc()
  1966. sce$SampleID = as.factor(sce$SampleID)
  1967. sce$Groups = as.factor(sce$Groups)
  1968. rm(cpm_mat)
  1969. gc()
  1970. ###################################
  1973. sce_sub <- sce[, sce$Groups %in% c("Resistant", "Sensitive")]
  1976. res_list = list()
  1978. clusters = unique(sce_sub$cell_type)
  1979. #clusters = 'Microglia'
  1980. set.seed(123)
  1981. gc()
  1982. for (cl in clusters) {
  1985. #cl = 'Microglia'
  1986. cat("Processing cluster:", cl, "\n")
  1988. sce_cl <- sce_sub[, sce_sub$cell_type == cl]
  1991. sce_cl$Groups = droplevels(sce_cl$Groups)
  1992. sce_cl$cell_type = droplevels(sce_cl$cell_type)
  1993. sample_md <- unique(as.data.frame(colData(sce_cl)[, c("SampleID", "Groups")]))
  1996. design_sub <- model.matrix(~ Groups, data = sample_md)
  1997. rownames(design_sub) <- sample_md$SampleID
  1998. design_sub
  1999. gc()
  2000. res_cl <- distinct_test(
  2001. x = sce_cl,
  2002. name_assays_expression = "cpm",
  2003. name_cluster = "cell_type", # single cluster, doesn't matter
  2004. name_sample = "SampleID",
  2005. design = design_sub,
  2006. column_to_test = 2,
  2007. min_non_zero_cells = 2000,
  2008. n_cores = 10,
  2009. P_4 = 5000
  2010. )
  2012. # Compute log2 FC
  2013. res_cl <- log2_FC(
  2014. res = res_cl,
  2015. x = sce_cl,
  2016. name_assays_expression = "cpm",
  2017. name_group = "Groups",
  2018. name_cluster = "cell_type"
  2019. )
  2021. # Store results with cluster name
  2022. res_list[[cl]] <- res_cl
  2023. gc()
  2024. }
  2026. res_sen <- bind_rows(res_list, .id = "cell_type")
  2029. top <- top_results(res_sen, significance = 0.05, global = TRUE)
  2030. write_xlsx(top, "res_sen_by_cluster.xlsx")
  2036. gc()
  2037. unique(sce_sub$Groups)
  2038. sce_sub$Groups = droplevels(sce_sub$Groups)
  2039. unique(sce_sub$Groups)
  2047. design_sub
  2048. gc()
  2050. sce_sub
  2052. sum(logcounts(sce_sub)["Sik2", ] > 0)
  2053. table(sce$cell_type)
  2055. res_sen <- distinct_test(
  2056. x = sce_sub,
  2057. name_assays_expression = "cpm",
  2058. name_cluster = "cell_type",
  2059. name_sample = "SampleID",
  2060. design = design_sub,
  2061. column_to_test = 2,
  2062. min_non_zero_cells = 200,
  2063. n_cores = 1,
  2064. P_4 = 20000
  2065. )
  2068. res_sen = log2_FC(res = res_sen,
  2069. x = sce_sub,
  2070. name_assays_expression = 'cpm',
  2071. name_group = "Groups",
  2072. name_cluster = "cell_subtype")
  2076. top = top_results(res_sen,
  2077. significance = 0.05, global = T)
  2079. plot_densities(x = sce_sub,
  2080. gene = 'Cdk10',
  2081. cluster = 'Drd2_Calcrl',
  2082. name_assays_expression = 'cpm',
  2083. name_cluster = 'cell_subtype',
  2084. name_sample = 'SampleID',
  2085. name_group = 'Groups', group_level = F)
  2087. plot_densities(x = sce_sub,
  2088. gene = 'Cdk10',
  2089. cluster = 'Drd2_Calcrl',
  2090. name_assays_expression = 'cpm',
  2091. name_cluster = 'cell_subtype',
  2092. name_sample = 'SampleID',
  2093. name_group = 'Groups', group_level = T)
  2095. write_xlsx(x = top, 'res_sen.xlsx')
  2097. #########################################
  2099. sce_sub <- sce[, sce$Groups %in% c("Resistant", "Yoked")]
  2102. unique(sce_sub$Groups)
  2103. sce_sub$Groups = droplevels(sce_sub$Groups)
  2104. unique(sce_sub$Groups)
  2106. set.seed(123)
  2109. sample_md <- unique(as.data.frame(colData(sce_sub)[, c("SampleID", "Groups")]))
  2112. design_sub <- model.matrix(~ Groups, data = sample_md)
  2113. rownames(design_sub) <- sample_md$SampleID
  2116. design_sub
  2117. gc()
  2119. sce_sub
  2121. sum(logcounts(sce_sub)["Slc17a7", ] > 0)
  2124. res_yok <- distinct_test(
  2125. x = sce_sub,
  2126. name_assays_expression = "cpm",
  2127. name_cluster = "cell_subtype",
  2128. name_sample = "SampleID",
  2129. design = design_sub,
  2130. column_to_test = 2,
  2131. min_non_zero_cells = 2000,
  2132. n_cores = 10,
  2133. P_4 = 20000
  2134. )
  2137. res_yok = log2_FC(res = res_yok,
  2138. x = sce_sub,
  2139. name_assays_expression = 'cpm',
  2140. name_group = "Groups",
  2141. name_cluster = "cell_subtype")
  2145. top = top_results(res_yok,
  2146. significance = 0.05, global = T)
  2148. plot_densities(x = sce_sub,
  2149. gene = 'Cdk10',
  2150. cluster = 'Drd2_Calcrl',
  2151. name_assays_expression = 'cpm',
  2152. name_cluster = 'cell_subtype',
  2153. name_sample = 'SampleID',
  2154. name_group = 'Groups', group_level = F)
  2156. plot_densities(x = sce_sub,
  2157. gene = 'Cdk10',
  2158. cluster = 'Drd2_Calcrl',
  2159. name_assays_expression = 'cpm',
  2160. name_cluster = 'cell_subtype',
  2161. name_sample = 'SampleID',
  2162. name_group = 'Groups', group_level = T)
  2165. write_xlsx(x = top, 'res_yok.xlsx')
  2166. gc()
  2168. ######################### Analysis of results #################################
  2170. res_sen = read_excel('res_sen.xlsx')
  2171. res_yok = read_excel('res_yok.xlsx')
  2174. head(res_sen)
  2175. head(res_yok)
  2178. library(dplyr)
  2180. overlap <- inner_join(
  2181. res_sen %>% select(gene, cluster_id, everything()),
  2182. res_yok %>% select(gene, cluster_id, everything()),
  2183. by = c("gene", "cluster_id"),
  2184. suffix = c("_sen", "_yok")
  2185. )
  2190. overlap_same_direction <- overlap %>%
  2191. filter(
  2192. (`log2FC_Resistant/Sensitive` > 0 & `log2FC_Resistant/Yoked` > 0) |
  2193. (`log2FC_Resistant/Sensitive` < 0 & `log2FC_Resistant/Yoked` < 0)
  2194. )
  2196. overlap_same_direction$average_log2FC = (overlap_same_direction$`log2FC_Resistant/Sensitive`+overlap_same_direction$`log2FC_Resistant/Yoked`)/2
  2197. overlap_same_direction$average_p_adj = (overlap_same_direction$p_adj.glb_sen + overlap_same_direction$p_adj.glb_yok)/2
  2200. write_xlsx(overlap_same_direction, 'DEG_list_males_main_celltypes.xlsx')
  2202. for_ipa = overlap_same_direction %>% select(gene, average_log2FC, average_p_adj, cluster_id)
  2204. dir.create('for_ipa')
  2205. celltypes <- unique(df$cluster_id)
  2207. for (ct in celltypes) {
  2208. sub <- for_ipa %>% filter(cluster_id == ct)
  2209. write.csv(sub, paste0("for_ipa/IPA_", ct, ".csv"), row.names = FALSE)
  2210. }
  2215. library(dplyr)
  2216. library(ggpubr)
  2217. library(viridis)
  2219. ##########################################################################
  2220. overlap_same_direction = read_xlsx('DEG_list_males_main_celltypes.xlsx')
  2223. library(dplyr)
  2224. library(ggpubr)
  2225. library(viridis)
  2226. library(dplyr)
  2227. library(ggplot2)
  2228. library(patchwork)
  2229. library(ggrepel)
  2232. df = overlap_same_direction
  2234. df <- df %>%
  2235. mutate(
  2236. reg_sen = case_when(
  2237. p_adj.glb_sen < 0.05 & `log2FC_Resistant/Sensitive` > 0.1 ~ "Up",
  2238. p_adj.glb_sen < 0.05 & `log2FC_Resistant/Sensitive` < -0.1 ~ "Down",
  2239. TRUE ~ "NS"
  2240. ),
  2241. reg_yok = case_when(
  2242. p_adj.glb_yok < 0.05 & `log2FC_Resistant/Yoked` > 0.1 ~ "Up",
  2243. p_adj.glb_yok < 0.05 & `log2FC_Resistant/Yoked` < -0.1 ~ "Down",
  2244. TRUE ~ "NS"
  2245. )
  2246. )
  2249. make_volcano_pair <- function(df, cluster) {
  2251. df_clust <- df %>% filter(cluster_id == cluster)
  2253. p_sen <- ggplot(df_clust, aes(
  2254. x = `log2FC_Resistant/Sensitive`,
  2255. y = -log10(p_adj.glb_sen)
  2256. )) +
  2257. geom_point(aes(color = reg_sen), size = 1.2, alpha = 0.7) +
  2258. geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") +
  2259. geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  2260. geom_text_repel(
  2261. data = df_clust %>% filter(reg_sen != "NS"),
  2262. aes(label = gene),
  2263. size = 3,
  2264. max.overlaps = 20
  2265. ) +
  2266. scale_color_manual(
  2267. values = c(Up = "red", Down = "blue", NS = "grey70")
  2268. ) +
  2269. labs(
  2270. title = "Resistant vs Sensitive",
  2271. x = "log2FC",
  2272. y = "-log10(adj p-value)"
  2273. ) +
  2274. theme_classic() +
  2275. theme(legend.position = "none")
  2277. p_yok <- ggplot(df_clust, aes(
  2278. x = `log2FC_Resistant/Yoked`,
  2279. y = -log10(p_adj.glb_yok)
  2280. )) +
  2281. geom_point(aes(color = reg_yok), size = 1.2, alpha = 0.7) +
  2282. geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") +
  2283. geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  2284. geom_text_repel(
  2285. data = df_clust %>% filter(reg_yok != "NS"),
  2286. aes(label = gene),
  2287. size = 3,
  2288. max.overlaps = 20
  2289. ) +
  2290. scale_color_manual(
  2291. values = c(Up = "red", Down = "blue", NS = "grey70")
  2292. ) +
  2293. labs(
  2294. title = "Resistant vs Yoked",
  2295. x = "log2FC",
  2296. y = "-log10(adj p-value)"
  2297. ) +
  2298. theme_classic() +
  2299. theme(legend.position = "none")
  2301. (p_sen | p_yok) + plot_annotation(title = cluster)
  2302. }
  2306. clusters <- unique(df$cluster_id)
  2308. volcano_plots <- lapply(clusters, function(cl) {
  2309. make_volcano_pair(df, cl)
  2310. })
  2314. volcano_plots[[1]]
  2316. dir.create('volcano_plots')
  2319. clusters <- unique(df$cluster_id)
  2321. for (cl in clusters) {
  2323. p <- make_volcano_pair(df, cl)
  2325. ggsave(
  2326. filename = paste0("volcano_plots/volcano_", cl, ".pdf"),
  2327. plot = p,
  2328. width = 15,
  2329. height = 10,
  2330. units = "in"
  2331. )
  2332. }
  2336. ################################################################################
  2338. df_plot <- overlap_same_direction %>% count(cluster_id)
  2340. p = ggbarplot(
  2341. df_plot,
  2342. x = "cluster_id",
  2343. y = "n",
  2344. fill = "cluster_id",
  2345. palette = viridis(nrow(df_plot)),
  2346. label = TRUE,
  2347. lab.vjust = -0.5,
  2348. sort.val = "desc", # sort bars by height (largest at top)
  2349. sort.by.groups = FALSE,
  2350. x.text.angle = 45,
  2351. legend = "none" # remove legend
  2352. ) +
  2353. labs(
  2354. title = "Number of DEGs per cluster",
  2355. x = "Cluster ID",
  2356. y = "Gene Count"
  2357. )
  2358. p
  2359. pdf('number_of_degs_per_cluster.pdf')
  2360. print(p)
  2361. dev.off()

master.R at commit 115c830, under MIT · at the source

Overview

Authors: Estelle Barbier1, Leon Höglund1, Li Xu1, Nihal A Salem2, Elizabeth A Osterndorff-Kahanek2, Riccardo Barchiesi2, Antonio Lacorte3,4, Lief E Fenno2, Robert O Messing2, R Dayne Mayfield2, Markus Heilig1
ORCID iDs: Estelle Barbier, Li Xu, Nihal A Salem, Elizabeth A Osterndorff-Kahanek, Riccardo Barchiesi, Antonio Lacorte, Lief E Fenno, Robert O Messing, R Dayne Mayfield, Markus Heilig
  1. Center for Social and Affective Neuroscience, BKV, Linköping University, Linköping, Sweden
  2. Waggoner Center for Alcohol and Addiction Research and Departments of Neuroscience and Neurology, University of Texas at Austin, Austin, TX USA
  3. Department of Pharmacy and Biotechnology, Alma Mater Studiorum – University of Bologna, Bologna, Italy
  4. Center for Neuroscience, University of Camerino, Camerino, Italy
Institutions: Linköping University (Sweden); The University of Texas at Austin (United States); Università di Camerino (Italy); University of Bologna (Italy)
Journal: Nature communications, volume 17, issue 1, article 9407
Dates: received 7 April 2026; accepted 20 August 2026; published online 2 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-77302-9 · PMID 42686758 · PMCID PMC13538410 · OpenAlex W7206173153
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), rat (organism), other condition (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials, Spectral & time-frequency
Keywords: Molecular neuroscience, Addiction
MeSH: Alcoholism*, Central Amygdaloid Nucleus*, Receptor, Fibroblast Growth Factor, Type 1*, Animals, Disease Models, Animal, Ethanol, Male, Neurons, Protein Kinase C-delta, Pyrimidines, Rats, Self Administration (* major topic)
Topic: Neurotransmitter Receptor Influence on Behavior (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: Vetenskapsrådet (Swedish Research Council) (2019-01138, 2013-07434, 2022-06568, 2024-02537); NIAAA NIH HHS (U01 AA013520, K00 AA029955, R01 AA012404, R01 AA028807, U01 AA020926); U.S. Department of Health & Human Services | National Institutes of Health (NIH) (K00AA029955, K00AA029955 (NAS), R01AA012404, U01AA020926, U01AA020926(RDM), R01AA012404(RDM)); U.S. Department of Health & Human Services | NIH | National Institute on Alcohol Abuse and Alcoholism (NIAAA) (AA013520, AA028807)
Citations: not cited yet (Europe PMC); 37 references in the paper
Research resources: The University of Texas at Austin RRID:SCR_021713

Abstract

A significant minority of alcohol users develop alcohol addiction, characterized by continued use despite negative consequences, referred to as compulsive-like. We previously showed that vulnerability to compulsive-like alcohol use can be modeled in rats using punished alcohol self-administration and in male rats is mediated by PKCδ+ neurons in the central nucleus of the amygdala (CeA). Here, we used cell-type-specific transcriptomics to identify molecular mechanisms underlying individual differences in this behavior. Transcriptional changes were restricted to a limited number of CeA neuronal populations, including PKCδ+ neurons, where weighted Gene Co-expression Network Analysis identified an upregulated co-expression module in punishment-resistant rats with FGFR1 as a druggable upstream regulator. Selective silencing of FgfR1 in PKCδ+ neurons normalized elevated PKCδ+ expression, and reduced punishment-resistant alcohol self-administration. This effect was recapitulated by systemic administration of the FgfR1-antagonist PD173074. These findings identify distinct CeA circuits that promote addiction vulnerability, and position FGFR1 as potential therapeutic targets.

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

leohog/Fgfr1_analysis

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 115c830f9e669f8eae77446eaca0d9d75e1e7ab0, 20 April 2026
Languages: R (1)
Size: 4 files, 1 script
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: cowplot (1 file), ggplot2 (1 file), ggpubr (1 file), igraph (1 file), patchwork (1 file), Seurat (1 file), SingleCellExperiment (1 file), tidyverse (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
3 files
At the source: github.com/leohog/Fgfr1_analysis

Code availability

The code used to perform these analyses is publicly available on GitHub at https://github.com/leohog/Fgfr1_analysis. For package version information, see supplementary data 3.

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

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;
  • 1 script, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

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

Data

Datasets cited

Data Availability Statement

The RNA sequencing data generated in this study have been deposited in the GEO database under accession code GSE310229 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE310229) https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE310229. Source data for all figures is provided in this paper. Source data are provided with this paper.

The code used to perform these analyses is publicly available on GitHub at https://github.com/leohog/Fgfr1_analysis. For package version information, see supplementary data 3.

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, 11 authors, 2 keywords, 12 MeSH terms, 4 funders, 36 references, 1 RRID.

Cite

This paper

Barbier, E., Höglund, L., Xu, L., Salem, N. A., Osterndorff-Kahanek, E. A., Barchiesi, R., Lacorte, A., Fenno, L. E., Messing, R. O., Mayfield, R. D., & Heilig, M. (2026). Identification and validation of central amygdala FGFR1 as a therapeutic target for alcohol use disorder using single-nucleus sequencing in rats. Nature communications, 17(1), 9407. https://doi.org/10.1038/s41467-026-77302-9

BibTeX

@article{barbier2026identification,
author = {Barbier, Estelle and Höglund, Leon and Xu, Li and Salem, Nihal A and Osterndorff-Kahanek, Elizabeth A and Barchiesi, Riccardo and Lacorte, Antonio and Fenno, Lief E and Messing, Robert O and Mayfield, R Dayne and Heilig, Markus},
title = {{Identification and validation of central amygdala FGFR1 as a therapeutic target for alcohol use disorder using single-nucleus sequencing in rats}},
journal = {Nature communications},
year = {2026},
month = sep,
volume = {17},
number = {1},
pages = {9407},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-77302-9},
url = {https://doi.org/10.1038/s41467-026-77302-9},
pmid = {42686758},
pmcid = {PMC13538410}
}

RIS

TY - JOUR
AU - Barbier, Estelle
AU - Höglund, Leon
AU - Xu, Li
AU - Salem, Nihal A
AU - Osterndorff-Kahanek, Elizabeth A
AU - Barchiesi, Riccardo
AU - Lacorte, Antonio
AU - Fenno, Lief E
AU - Messing, Robert O
AU - Mayfield, R Dayne
AU - Heilig, Markus
TI - Identification and validation of central amygdala FGFR1 as a therapeutic target for alcohol use disorder using single-nucleus sequencing in rats
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/09/02
VL - 17
IS - 1
SP - 9407
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-77302-9
UR - https://doi.org/10.1038/s41467-026-77302-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-77302-9",
"type": "article-journal",
"title": "Identification and validation of central amygdala FGFR1 as a therapeutic target for alcohol use disorder using single-nucleus sequencing in rats",
"container-title": "Nature communications",
"author": [
{
"family": "Barbier",
"given": "Estelle"
},
{
"family": "Höglund",
"given": "Leon"
},
{
"family": "Xu",
"given": "Li"
},
{
"family": "Salem",
"given": "Nihal A"
},
{
"family": "Osterndorff-Kahanek",
"given": "Elizabeth A"
},
{
"family": "Barchiesi",
"given": "Riccardo"
},
{
"family": "Lacorte",
"given": "Antonio"
},
{
"family": "Fenno",
"given": "Lief E"
},
{
"family": "Messing",
"given": "Robert O"
},
{
"family": "Mayfield",
"given": "R Dayne"
},
{
"family": "Heilig",
"given": "Markus"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "9407",
"DOI": "10.1038/s41467-026-77302-9",
"PMID": "42686758",
"PMCID": "PMC13538410",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-77302-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
2
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/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: WGCNA, SingleCellExperiment, igraph, 6 other tools, genetics / omics, cellular / molecular, 2 references
[2] 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: WGCNA, SingleCellExperiment, igraph, 6 other tools, genetics / omics, other condition, cellular / molecular, 1 reference
[3] doi:10.1038/s44318-026-00806-z [code]
Interspecific diversity in the neuronal composition of the mammalian cortex arises from heterochrony in neurogenesis.
Journal: The EMBO journal
In common: WGCNA, SingleCellExperiment, igraph, 6 other tools, rat, 1 reference
[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: WGCNA, SingleCellExperiment, igraph, 6 other tools, cellular / molecular, 1 reference
[5] 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, SingleCellExperiment, igraph, 6 other tools, cellular / molecular, 1 reference
[6] doi:10.1038/s41467-026-71542-5 [code]
Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder.
Journal: Nature communications
In common: WGCNA, SingleCellExperiment, igraph, 5 other tools, genetics / omics, cellular / molecular, 1 reference
[7] doi:10.1038/s41593-026-02384-z [code]
cGAS-mediated type I IFN signaling contributes to disease progression in drug-refractory epilepsy.
Journal: Nature neuroscience
In common: WGCNA, SingleCellExperiment, igraph, 5 other tools, cellular / molecular, 1 reference
[8] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: WGCNA, SingleCellExperiment, Seurat, 5 other tools, genetics / omics, cellular / molecular, 1 reference
[9] doi:10.1016/j.isci.2026.115196 [code]
Transcriptional and cellular maturation of the chick spinal cord in the context of distinct neuromuscular circuits.
Journal: iScience
In common: WGCNA, SingleCellExperiment, igraph, 6 other tools
[10] doi:10.1038/s41380-026-03629-w [code]
Maternal fasting during early gestation induces epigenetic alterations and schizophrenia-related phenotypes.
Journal: Molecular psychiatry
In common: WGCNA, igraph, Seurat, 4 other tools, genetics / omics, cellular / molecular, 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.