OSCR

Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain.

Code ↔ Paper

9 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 9 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › Sample collection ↔ 01b.ent_process_outerMerge_0807.ipynb, lines 47–71 · score 0.77 · superior septum, superior turbinate, middle turbinate, history, surgery, metadata
  2. [2] § Results › Distinct stages of OSNs in adult olfactory epithelium ↔ 02a.traj_integrate_600-6-40-50-graph_k7_FDR_BG-Cov.ipynb, lines 497–515 · score 0.75 · receptor transformations, GBCinp, olfactory neurogenesis, SOX4, precursor, NEUROG1
  3. [3] § Results › Distinct stages of OSNs in adult olfactory epithelium ↔ 01b.ent_process_outerMerge_0807.ipynb, lines 47–71 · score 0.68 · right superior turbinate, superior septum, middle turbinate, surgery
  4. [4] § Methods › Trajectory inference and trajDEG identification ↔ 02a.traj_integrate_600-6-40-50-graph_k7_FDR_BG-Cov.ipynb, lines 76–143 · score 0.65 · graph_test, Moran, fitting, FDR, age, Monocle3
  5. [5] § Results › Distinct stages of OSNs in adult olfactory epithelium ↔ 01d.statstitic.ipynb, lines 22–46 · score 0.61 · cell cycle, l_iOSN, e_iOSN, phase, scores, GBCs
  6. [6] § Methods › snRNA-seq data preprocessing and clustering ↔ 02a.traj_integrate_600-6-40-50-graph_k7_FDR_BG-Cov.ipynb, lines 429–475 · score 0.57 · cell clustering, graph, PCA, PCs, neighborhood, filtered
  7. [7] § Results › Similarities and differences between OSNs and CENs ↔ 02a.traj_integrate_600-6-40-50-graph_k7_FDR_BG-Cov.ipynb, lines 1000–1017 · score 0.56 · trajDEG, DEG trends, biological processes, enriched, genes
  8. [8] § Methods › Genetic enrichment ↔ 01e.scdrs_downstream_stage.R, the whole file · a weak match · score 0.55 · ALS, stroke, ADHD, MDD, MS, PD
  9. [9] § Methods › Data integration ↔ 02a.traj_integrate_600-6-40-50-graph_k7_FDR_BG-Cov.ipynb, lines 429–475 · score 0.53 · IntegrateLayers, anchor, PCs, neighbors, CCA, variables

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

Jupyter notebook · 1,017 lines · 37 KB · no license · 5 matches

  1. # %%
  2. options(repr.plot.width = 10, repr.plot.height = 6, repr.plot.res = 100)
  4. setwd('/sc/arion/projects/roussp01a/liting/Olf')
  6. library(rlang)
  7. library(RColorBrewer)
  8. library(scales)
  9. library(tradeSeq)
  10. library(schard)
  11. library(SeuratWrappers)
  12. library(ggVennDiagram)
  13. library(stringr)
  14. library(rrvgo)
  15. library(slingshot)
  16. library(ComplexHeatmap)
  17. library(circlize)
  18. library(dplyr)
  19. library(gprofiler2)
  20. library(monocle3)
  21. library(ggplot2)
  23. # %%
  24. # functions
  25. scanpy2seurat <- function(file_name){
  26. h5ad=paste0('/sc/arion/projects/CommonMind/roussp01a/ENT/snRNAseq/qc_scanpy/',file_name,'.h5ad')
  27. rds=paste0('/sc/arion/projects/CommonMind/roussp01a/ENT/snRNAseq/qc_scanpy/',file_name,'.rds')
  28. ent = schard::h5ad2seurat(h5ad)
  29. saveRDS(ent,rds)
  30. }
  32. read_obj <- function(file_name, hvg){
  33. ent_N <- readRDS(paste0('/sc/arion/projects/CommonMind/roussp01a/ENT/snRNAseq/qc_scanpy/',file_name,'.rds'))
  34. names(ent_N@reductions) <- gsub('^X|_$','',names(ent_N@reductions))
  35. ent_N <- FindVariableFeatures(ent_N,nfeatures = hvg) # for identify DEGs using slingshot
  36. return(ent_N)
  37. }
  39. #N_color <- c('#1f77b4','#ff7f0e',"#D32B29")
  40. N_color <- c('#1f77b4',"#ea801c","#D32B29")
  41. names(N_color) <- c('GBC','iOSN','mOSN')
  44. get_slingshot_traj <- function(obj){
  45. dimred <- obj@reductions$[email hidden]
  46. #clustering <- as.numeric(obj$N_leiden_res0_25)+1
  47. clustering <- obj$N_types
  49. set.seed(1)
  50. pto <- slingshot(dimred, clustering, start.clus = 'GBC')
  51. obj$sling_pseudotime <- slingPseudotime(pto)
  52. FeaturePlot(object = obj, features = 'sling_pseudotime',reduction = "umap.cca")
  54. png(filename = "./figures/outermerge_slingshot_pseudotime.png",width = 8, height = 8, units = "cm", res=300)
  56. plot(dimred[, 1:2], col = N_color[clustering], cex = 0.5, pch = 16)
  57. lines(SlingshotDataSet(pto), lwd=2, col='black')
  58. dev.off()
  59. return(obj)
  60. }
  62. tradeseq_DEG <- function(data_ns, traj_method){
  64. data_ns$cca_clusters <- obj$cca_clusters
  65. counts <- as.matrix(data_ns@assays$RNA@counts[data_ns@assays$[email hidden], ])
  66. filt_counts <- counts[rowSums(counts > 0) > ncol(counts)/100, ]
  67. Pseudotime <- as.matrix(data_ns$sling_pseudotime)
  69. batch <- data_ns$batch
  70. U <- model.matrix(~batch)
  72. sce <- fitGAM(counts = as.matrix(filt_counts),
  73. pseudotime = Pseudotime,U=U,
  74. cellWeights = as.matrix(rep(1,ncol(data_ns))))
  76. # test for dynamic expression
  77. pseudotime_association <- associationTest(sce)
  79. pseudotime_association$fdr <- p.adjust(pseudotime_association$pvalue, method = "fdr")
  80. pseudotime_association <- pseudotime_association[order(pseudotime_association$pvalue), ]
  81. pseudotime_association$feature_id <- rownames(pseudotime_association)
  82. pseudotime_association_sig <- subset(pseudotime_association, fdr < 0.05)
  84. return(pseudotime_association_sig)
  85. #return(pseudotime_association)
  86. }
  90. get_monocle3_traj <- function(scdata,lab,mm, res_traj){
  91. # input: scanpy data
  92. # get root node from plot_cells
  94. # as Monocle3 data
  95. mnc3_data <- SeuratWrappers::as.cell_data_set(scdata)
  96. mnc3_data <- estimate_size_factors(mnc3_data)
  98. # Cluster your cells
  99. mnc3_data <- cluster_cells(mnc3_data,resolution=res_traj)
  100. mnc3_data <- learn_graph(mnc3_data) # # Learn the trajectory graph
  101. plot_cells(mnc3_data, color_cells_by = "cca_N_types_stage", cell_size = 1, label_principal_points = F, label_leaves=F,
  102. label_branch_points=F)
  104. print('-------------check root cell based on umap--------------------')
  105. if(lab=='x' & mm=='max'){rootcell = names(which.max(subset(data_ns, cca_N_types =='GBC')@reductions$[email hidden][,1])) }
  106. if(lab=='x' & mm=='min'){rootcell = names(which.min(subset(data_ns, cca_N_types =='GBC')@reductions$[email hidden][,1])) }
  107. if(lab=='y' & mm=='max'){rootcell = names(which.max(subset(data_ns, cca_N_types =='GBC')@reductions$[email hidden][,2])) }
  108. if(lab=='y' & mm=='min'){rootcell = names(which.min(subset(data_ns, cca_N_types =='GBC')@reductions$[email hidden][,2])) }
  110. mnc3_data <- order_cells(mnc3_data, root_cells= rootcell )
  111. plot_cells(mnc3_data, color_cells_by = "pseudotime", cell_size = 1)
  113. return(mnc3_data)
  114. }
  116. get_monocle3_DEG <- function(mnc3_data){
  118. pr_test_res <- monocle3:::graph_test(mnc3_data, neighbor_graph="principal_graph", cores=6)
  120. pr_test_res <- pr_test_res[order(pr_test_res$q_value),]
  121. pr_test_res_sig <- subset(pr_test_res ) # get_monocle3_DEG #& morans_I > 0.1 & morans_I > 0.05
  122. return(pr_test_res_sig)
  123. #return(pr_test_res)
  124. }
  126. get_monocle3_DEG_contBatch <- function(mnc3_data){
  127. mnc3_data$Age <- scale(mnc3_data$Age)
  128. pr_test_res <- fit_models(mnc3_data,model_formula_str = "~monocle3_pseudotime + batch + Sex + Age")#expression_family="negbinomial",
  129. pr_test_res <- coefficient_table(pr_test_res)
  130. pr_test_res <- pr_test_res %>% filter(term == "monocle3_pseudotime") %>%
  131. select(gene_id, term, q_value, estimate)
  132. pr_test_res <- pr_test_res[order(pr_test_res$q_value),]
  133. pr_test_res_sig <- subset(pr_test_res, q_value < 0.05 )
  134. pr_test_res_sig <- as.data.frame(pr_test_res_sig)
  135. rownames(pr_test_res_sig) <- pr_test_res_sig$gene_id
  138. return(pr_test_res_sig)
  139. #return(pr_test_res)
  140. }
  143. get_heatmap_re_cluster_kmeans <- function(scdata, method, DEGs, n_split, recluster_id,n_REsplit){
  145. #if (method=='monocle') { pseudoTime <- pseudotime(scdata) }
  146. if (method=='monocle3') { pseudoTime <- scdata$monocle3_pseudotime}
  147. if (method=='slingshot') { pseudoTime <- scdata$sling_pseudotime }
  148. if (method=='palantir') { pseudoTime = scdata$palantir_pseudotime}
  149. if (method=='paga') { pseudoTime = scdata$dpt_pseudotime}
  151. pt.matrix1 <- as.matrix(scdata@assays$RNA@counts[rownames(DEGs),order(pseudoTime)])
  152. #pt.matrix1 <- as.matrix(scdata[["RNA"]]$data[rownames(DEGs),order(pseudoTime)])
  153. #pt.matrix1 <- as.matrix(scdata[["RNA"]]$scale.data[rownames(DEGs),order(pseudoTime)])
  155. #Can also use "normalized_counts" instead of "exprs" to use various normalization methods, for example:
  156. #normalized_counts(cds, norm_method = "log")
  157. pt.matrix <- pt.matrix1
  158. pt.matrix <- t(apply(pt.matrix,1,function(x){smooth.spline(x,df=3)$y}))
  159. pt.matrix <- t(apply(pt.matrix,1,function(x){(x-mean(x))/sd(x)}))
  161. pt.matrix[pt.matrix > 4] = 4
  162. pt.matrix[pt.matrix < -4] = -4
  163. #
  165. rownames(pt.matrix) <- rownames(DEGs);
  166. colnames(pt.matrix) <- colnames(scdata)[order(pseudoTime)]
  168. set.seed(42)
  170. ks_1st <- kmeans(pt.matrix, centers = n_split , iter.max = 100)
  171. C_K1 <- ks_1st$cluster
  172. C_K2 <- ''
  173. if (recluster_id!=''){
  174. ks_2st <- kmeans(pt.matrix[ks_1st$cluster==recluster_id,], centers = n_REsplit, iter.max = 100)
  175. C_K2 <- ks_2st$cluster[names(ks_1st$cluster)]
  176. C_K2 <- ifelse(is.na(C_K2),'',paste0('-',C_K2))
  178. }
  181. rowsplit <- paste0(C_K1,C_K2)
  182. dict_cluster <- c('iOSN','GBC','mOSN')
  183. names(dict_cluster) <- names(sort(table(rowsplit)))
  184. rowsplit=dict_cluster[rowsplit]
  187. hthc=ComplexHeatmap::draw(Heatmap(
  188. pt.matrix,
  189. row_split = rowsplit,
  190. col = colorRamp2(seq(from=-2,to=2,length=11),rev(brewer.pal(11, "Spectral"))),
  191. #row_title = "cluster_%s",
  192. row_gap = unit(c(2.5), "mm"),
  193. cluster_row_slices = FALSE,
  194. cluster_columns = F,
  195. show_row_names = FALSE,
  196. show_column_names = FALSE,
  197. show_row_dend = F,
  198. top_annotation = HeatmapAnnotation(
  199. cell_identity=[email hidden][colnames(pt.matrix),'N_types'],
  200. cell_subcluster=[email hidden][colnames(pt.matrix),'leiden'],
  201. #batch=[email hidden][colnames(pt.matrix),'batch'],
  202. #subc=[email hidden][colnames(pt.matrix),'N_leiden_res0_2'],
  203. pseudotime=sort(pseudoTime),
  204. col = list(cell_identity =N_color)#,
  205. # pseudotime= colorRamp2(c(0, 0.5, 1, 1.5), c("#35008C", "#8700A8", "#D3546F",'#E7F92D'))
  206. #)
  207. )
  208. ))
  211. deg_cl <- as.data.frame(cbind(as.character(rowsplit), rownames(DEGs),method))
  212. colnames(deg_cl) <- c('cluster','DEG','method')
  214. return(deg_cl)
  216. }
  218. library(enrichR)
  219. get_bp_enrichr <- function(query_gene){
  220. dbs <- "GO_Biological_Process_2023"
  221. enriched <- enrichr(query_gene, dbs)
  222. enrr <- enriched[['GO_Biological_Process_2023']]
  223. return(enrr)
  224. }
  228. get_bp_enrichr <- function(cl, md){
  229. query_gene <- subset(DEG_ent, cluster==cl & method==md)[,'DEG']
  230. dbs <- dbs <- c( "GO_Biological_Process_2023",
  231. "KEGG_2021_Human")
  233. enriched <- enrichr(query_gene, dbs)
  234. enrr <- rbind(enriched[['GO_Biological_Process_2023']])
  235. enrr$cluster <- cl
  236. enrr$method <- md
  237. return(enrr)
  238. }
  240. get_bp_enrichr_reducedTerms <- function(DEG_ent, cl,md){
  242. query_gene <- subset(DEG_ent, cluster==cl & method==md)[,'DEG']
  243. dbs <- c( "GO_Biological_Process_2023")
  244. #"KEGG_2021_Human")
  246. enriched <- enrichr(query_gene, dbs)
  247. enrr <- enriched[['GO_Biological_Process_2023']]
  248. enrr$genecluster <- cl
  249. enrr$method <- md
  250. go_analysis <- subset(enrr, p_value < 0.05 )
  251. go_analysis$GOID <- str_split(go_analysis$Term,'\\(|\\)',simplify = T)[,2]
  252. simMatrix <- calculateSimMatrix(go_analysis$GOID,
  253. orgdb="org.Hs.eg.db",
  254. ont="BP",
  255. method="Rel")
  258. scores <- setNames(-log10(as.numeric(go_analysis$p_value)), go_analysis$GOID)
  259. reducedTerms <- reduceSimMatrix(simMatrix,
  260. scores,
  261. orgdb="org.Hs.eg.db")
  263. treemapPlot(reducedTerms)
  265. reducedTerms_sumscore <- aggregate(score~parentTerm,reducedTerms, sum)
  266. reducedTerms_sumscore$genecluster <- cl
  267. reducedTerms_sumscore$method <- md
  268. reducedTerms_sumscore$term <- reducedTerms_sumscore$parentTerm
  271. reducedTerms <- reducedTerms%>%group_by(cluster)%>%top_n(1,score)
  272. reducedTerms$genecluster <- cl
  273. reducedTerms$method <- md
  274. reducedTerms$term <- reducedTerms$parentTerm
  277. return(list(enrr=enrr,reducedTerms=reducedTerms,reducedTerms_sumscore=reducedTerms_sumscore))
  279. }
  282. #table(ent_nn$batch)
  284. # %%
  285. #colData(int_N)
  287. # %%
  288. graph_test_lm_ENT <- function (cds, neighbor_graph = c("knn", "principal_graph"),
  289. reduction_method = "UMAP", k = 25, method = c("Moran_I"),
  290. alternative = "greater", expression_family = "quasipoisson",
  291. cores = 1, verbose = FALSE)
  292. {
  294. #nn_control_default <- get_global_variable('nn_control_annoy_euclidean')
  295. nn_control_default <- list(method='annoy', metric='euclidean', n_trees=50, M=48, ef_construction=200, ef=150, grain_size=1, cores=1)
  297. nn_control <- monocle3:::set_nn_control(mode=3,
  299. nn_control_default=nn_control_default,
  300. nn_index=NULL,
  301. k=k,
  302. verbose=verbose)
  304. neighbor_graph <- match.arg(neighbor_graph)
  306. lw <- monocle3:::calculateLW(cds=cds,
  307. k = k,
  308. neighbor_graph = neighbor_graph,
  309. reduction_method = reduction_method,
  310. verbose = verbose,
  311. nn_control = nn_control_default
  312. )
  315. exprs_mat <- SingleCellExperiment::counts(cds)[, attr(lw, "region.id"), drop = FALSE]
  316. sz <- size_factors(cds)[attr(lw, "region.id")]
  317. wc <- spdep::spweights.constants(lw, zero.policy = TRUE,
  318. adjust.n = TRUE)
  319. test_res <- pbmcapply::pbmclapply(row.names(exprs_mat), FUN = function(x,
  320. sz, alternative, method, expression_family) {
  321. exprs_val <- exprs_mat[x, ]
  322. if (expression_family %in% c("uninormal", "binomialff")) {
  323. exprs_val <- exprs_val
  324. }
  325. else {
  326. exprs_val <- log10(exprs_val/sz + 0.1)
  327. }
  328. df = cbind(as.data.frame(exprs_val), colData(cds)$Sex, scale(colData(cds)$Age),colData(cds)$batch ,log(colData(cds)$nCount_RNA))
  329. colnames(df) <- c("exp", "Sex", "Age",'batch','log1p_total_counts')
  330. test_res <- tryCatch({
  331. if (method == "Moran_I") {
  332. mt <- suppressWarnings(monocle3:::my.moran.test(df, lw, wc, alternative = alternative))
  333. data.frame(status = "OK", p_value = mt$p.value,
  334. morans_test_statistic = mt$statistic, morans_I = mt$estimate[["Moran I statistic"]])
  335. }
  336. else if (method == "Geary_C") {
  337. gt <- suppressWarnings(my.geary.test(exprs_val,
  338. lw, wc, alternative = alternative))
  339. data.frame(status = "OK", p_value = gt$p.value,
  340. geary_test_statistic = gt$statistic, geary_C = gt$estimate[["Geary C statistic"]])
  341. }
  342. }, error = function(e) {
  343. data.frame(status = "FAIL", p_value = NA, morans_test_statistic = NA,
  344. morans_I = NA)
  345. })
  346. }, sz = sz, alternative = alternative, method = method, expression_family = expression_family,
  347. mc.cores = cores, ignore.interactive = TRUE)
  348. if (verbose) {
  349. message("returning results: ...")
  350. }
  351. test_res <- do.call(rbind.data.frame, test_res)
  352. row.names(test_res) <- row.names(cds)
  353. test_res <- merge(test_res, rowData(cds), by = "row.names")
  354. row.names(test_res) <- test_res[, 1]
  355. test_res[, 1] <- NULL
  356. test_res$q_value <- 1
  357. test_res$q_value[which(test_res$status == "OK")] <- stats::p.adjust(subset(test_res,
  358. status == "OK")[, "p_value"], method = "BH")
  359. test_res$status = as.character(test_res$status)
  360. test_res[row.names(cds), ]
  361. }
  363. my.moran.test_lm_ENT <- function (x, listw, wc, alternative = "greater", randomisation = TRUE)
  364. {
  366. zero.policy = TRUE
  367. adjust.n = TRUE
  368. na.action = stats::na.fail
  369. drop.EI2 = FALSE
  370. xname <- deparse(substitute(x))
  371. wname <- deparse(substitute(listw))
  372. NAOK <- deparse(substitute(na.action)) == "na.pass"
  373. x <- na.action(x)
  374. na.act <- attr(x, "na.action")
  375. if (!is.null(na.act)) {
  376. subset <- !(1:length(listw$neighbours) %in% na.act)
  377. listw <- subset(listw, subset, zero.policy = zero.policy)
  378. }
  379. n <- length(listw$neighbours)
  380. S02 <- wc$S0 * wc$S0
  381. model <- lm(exp ~ Age + Sex + batch + log1p_total_counts, data = x) # batch + Sex + Age +as.numeric(log1p_total_counts)
  382. res <- spdep::lm.morantest(model, listw, zero.policy = zero.policy, alternative = alternative)
  383. statistic = as.numeric(res[1])
  384. names(statistic) <- "Moran I statistic standard deviate"
  385. PrI = as.numeric(res[2])
  386. vec <- c(res[3]$estimate[1], res[3]$estimate[2], res[3]$estimate[3])
  387. names(vec) <- c("Moran I statistic", "Expectation", "Variance")
  388. method <- paste("Moran I test under", ifelse(randomisation,
  389. "randomisation", "normality"))
  390. res <- list(statistic = statistic, p.value = PrI, estimate = vec)
  391. if (!is.null(na.act))
  392. attr(res, "na.action") <- na.act
  393. class(res) <- "htest"
  394. res
  395. }
  398. assignInNamespace(x = "graph_test", value = graph_test_lm_ENT, ns = "monocle3")
  400. assignInNamespace(x = "my.moran.test", value = my.moran.test_lm_ENT, ns = "monocle3")
  404. # %%
  405. #xtabs(~batch+N_types,[email hidden])
  407. # %% [markdown]
  408. # ### CCA data integration
  410. # %%
  411. # 1. read tata
  412. # from 2_integrate_Neuron_nn_ent_outer
  414. scanpy2seurat('ent_nn_merge_rawcount')
  415. ent_nn <- read_obj('ent_nn_merge_rawcount', hvg = 3000)
  418. ent_nn=ent_nn[!grepl('^RPS|^RPL|^LINC|^MT',rownames(ent_nn)),]
  420. #data <- subset(data,batch%in%c(setdiff(unique(data$batch),c('Set4_C1','Set4_C2','Set1_C1','Set2_C1','Set2_C2') )))
  421. #ent_nn$batch <- str_split(ent_nn$batch,'_',simplify = T)[,1]
  422. ent_nn$batch <- ent_nn$Set
  424. ent_nn[["RNA"]] <- split(ent_nn[["RNA"]], f = ent_nn$batch)
  425. #data[["RNA"]] <- split(data[["RNA"]], f = data$dataset)
  426. ent_nn <- subset(ent_nn,batch%in%c(setdiff(unique(ent_nn$batch),c('Set2', 'Set3', 'Set4') )))#
  429. # %%
  431. n_features <- c(600)
  432. n_pcs <- c(6)
  433. n_neighbors <- c(40)
  434. n_kweight <- c(50)
  436. for (n_pc in n_pcs){
  437. for (n_feature in n_features){
  438. for (n_neighbor in n_neighbors){
  439. for(kw in n_kweight){
  442. data <- NormalizeData(ent_nn)
  443. data <- FindVariableFeatures(data,selection.method='mean.var.plot',nfeatures=n_feature)
  444. data <- ScaleData(data,features=VariableFeatures(data))
  445. #run PCA. Select significant PCs based on a scree plot. Look for the last point before the plot becomes flat
  446. data <- RunPCA(data,features = VariableFeatures(data),verbose=F)
  448. # 4. batch correction
  449. ## 4.1 cca
  450. obj <- IntegrateLayers(
  451. object = data, method = CCAIntegration,k.weight =kw,
  452. orig.reduction = "pca", new.reduction = "integrated.cca",
  453. verbose = FALSE, dims = 1:n_pc
  454. # k.weight = 10,
  455. # k.anchor = 10,
  456. # k.filter = 10,
  457. # k.score = 10
  458. )#
  459. # 5 identify cell clusters
  460. ## 5.1 cca
  461. obj <- FindNeighbors(obj, reduction = "integrated.cca", dims = 1:n_pc)
  462. obj <- FindClusters(obj, resolution = 0.25, cluster.name = "cca_clusters")
  464. obj <- RunUMAP(obj, reduction = "integrated.cca", dims = 1:n_pc, reduction.name = "umap.cca", n.neighbors=n_neighbor)
  466. p1 <- DimPlot(
  467. obj,
  468. reduction = "umap.cca",
  469. group.by = c("batch" ,'dataset','cca_clusters'),label.size = 2
  470. )
  471. # ggsave(p1, file=paste0('./figures/integrated_pcs/',n_pc,'_',n_feature,"_",n_neighbor,"_",kw,'inte_umap.pdf'), width=14, height=4)
  473. }
  474. }
  475. }}
  480. # %%
  481. # label
  482. DimPlot(
  483. obj,
  484. reduction = "umap.cca",
  485. group.by = c('cca_clusters','N_types'),label.size = 2)
  487. # %%
  488. #DEG_trend_Olf['CALB1',]
  489. #FeaturePlot(obj, features = c('CALB1','CALB2') , cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_minimal()+theme_void()+theme(legend.position = '')
  491. #FeaturePlot(obj, features = DEG_trend_Olf$gene[DEG_trend_Olf$cl=='7'][14:18] , cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_minimal()+theme_void()+theme(legend.position = '')
  492. #scz_risk_genes <- c("COMT", "DISC1", "ZNF804A", "NRG1", "DTNBP1", "G72", "MAOA", "SLC6A4", "CACNA1C", "TTC28")
  494. #FeaturePlot(obj, pt.size = 1,features =scz_risk_genes[4], cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_minimal()+theme_void()+theme(legend.position = '')
  497. # %%
  498. # Single-cell transcriptomics reveals receptor transformations during olfactory neurogenesis
  499. # for progenitors, Ascl1 (achaete-scute complex homolog 1); for precursors, Neurog1 (neurogenin 1) and/or Neurod1 (neurogenic differentiation 1);
  500. GBCinp <- c("HES6","CXCR4","NEUROD1","NEUROG1")
  501. GBCprogenitors <- c('ASCL1','MKI67','TOP2A')
  502. iOSN <- c('GNG8', 'GAP43','LHX2','SOX4')#c('EBF2', 'EMX2')
  503. p1 <- FeaturePlot(obj, features = GBCinp[1], cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_void()+theme(legend.position = '', aspect.ratio = 1)
  504. p2 <- FeaturePlot(obj, features = GBCinp[2], cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_void()+theme(legend.position = '', aspect.ratio = 1)
  505. p3 <- FeaturePlot(obj, features = GBCinp[3], cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_void()+theme(legend.position = '', aspect.ratio = 1)
  506. p4 <- FeaturePlot(obj, features = GBCinp[4], cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_void()+theme(legend.position = '', aspect.ratio = 1)
  507. p5 <- FeaturePlot(obj, features = "LHX2", cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_void()+theme(legend.position = '', aspect.ratio = 1)
  508. p6 <- FeaturePlot(obj, features = "GAP43", cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_void()+theme(legend.position = '', aspect.ratio = 1)
  510. px <- cowplot::plot_grid(p1,p2,p3,p4,nrow=2)
  511. pX2 <- cowplot::plot_grid(p3,p4,p5,p6,nrow=2)
  513. #FeaturePlot(obj, features = GBCprogenitors , cols = c("#FFF5F0",'#F75D42', "#6A010D") )+theme_minimal()+theme_void()+theme(legend.position = '')
  514. ggsave(px, file=paste0('./figures/03GBC_markers.pdf'), width=5, height=5)
  515. ggsave(pX2, file=paste0('./figures/03GBC_markers2.pdf'), width=5, height=5)
  517. # %%
  518. # mesenchymal cells
  519. # iOSN <- c('ACTA2','MAP1B', 'COL1A2','FZD2', 'FZD7', 'ROR2', 'SFRP1' , 'SFRP2', 'CTNNB1' , 'JAG1', 'PSEN1' , 'APH1A')#c('EBF2', 'EMX2')
  520. # FeaturePlot(obj, features = iOSN, cols = c("#FFF5F0",'#F75D42', "#6A010D") )
  522. # iOSN <- c('TAGLN', 'COL1A2', 'COL1A1', 'CALD1', 'TPM2', 'COL3A1', 'TPM1', 'LGALS1')
  523. # FeaturePlot(obj, features = iOSN, cols = c("#FFF5F0",'#F75D42', "#6A010D") ,pt.size = 0.2)
  525. # iOSN <- c('UCHL1', 'MAP1A', 'MAP1B', 'TUBB3', 'INA', 'NRP1', 'MKI67', 'ACTB','NES')
  526. # FeaturePlot(obj, features = iOSN, cols = c("#FFF5F0",'#F75D42', "#6A010D") ,pt.size = 0.2)
  528. iOSN <- mkx <- c('UCHL1', 'MAP1A', 'MAP1B', 'TUBB3', 'INA', 'NRP1', 'MKI67', 'ACTB','NES')
  529. FeaturePlot(obj, features = iOSN, cols = c("#FFF5F0",'#F75D42', "#6A010D") ,pt.size = 0.2)
  530. STK6), PLK1, E2F1, FOXM1, MKI67
  532. # %%
  533. ent_N <- readRDS('/sc/arion/projects/CommonMind/roussp01a/ENT/snRNAseq/qc_scanpy/ent_nn_merge_cca.rds')
  535. # %%
  536. # # #obj <- FindNeighbors(obj, reduction = "integrated.cca", dims = 1:n_pc)
  537. # obj <- FindClusters(obj, resolution = 0.25, cluster.name = "cca_clusters")
  539. # obj <- RunUMAP(obj, reduction = "integrated.cca", dims = 1:n_pc, reduction.name = "umap.cca", n.neighbors=n_neighbor)
  541. # DimPlot(
  542. # obj,
  543. # reduction = "umap.cca",
  544. # group.by = c('cca_clusters'),label.size = 2
  545. # )
  547. cca_label <- c( 'GBC','e_iOSN','l_iOSN','mOSN')
  548. names(cca_label) <- c('3','2','0','1')
  549. obj <- RenameIdents(obj, cca_label)
  550. [email hidden]$cca_N_types_stage <- cca_label[as.character(obj$cca_clusters)]
  552. [email hidden]$cca_N_types <- ifelse([email hidden]$cca_N_types_stage%in%c('e_iOSN','l_iOSN'),'iOSN',[email hidden]$cca_N_types_stage)
  554. DimPlot(
  555. obj,
  556. reduction = "umap.cca",
  557. group.by = c('cca_N_types_stage', 'cca_N_types'),label.size = 2
  558. )
  560. saveRDS(obj,'/sc/arion/projects/CommonMind/roussp01a/ENT/snRNAseq/qc_scanpy/ent_nn_merge_cca.rds')
  562. # %%
  563. options(repr.plot.width = 10, repr.plot.height = 8, repr.plot.res = 100)
  565. FeaturePlot(obj, features = iOSN, cols = c("#FFF5F0",'#F75D42', "#6A010D") )#+theme( aspect.ratio = 1)
  567. DotPlot(obj, features = c('NEUROD1','NEUROG1','SOX4','LHX2','GNG8','GAP43','GNAL','GNG13'),
  568. cols = c("white", "#A40F14"),scale = FALSE) + RotatedAxis()+theme_bw()+
  569. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),legend.position = 'bottom')+ylab('Cell types')
  571. dev.print(pdf,file = "/sc/arion/projects/roussp01a/liting/Olf/figures/marker_dot.pdf", width = 3.5, height = 5)
  574. # %%
  575. write.csv([email hidden], file='./data/ent_nn_merge_cca.metadata.csv')
  577. # %%
  578. obj_combined <- JoinLayers(object = obj)
  579. #markers_2 <- FindMarkers(CB30_combined, ident.1 = "CD4+", ident.2 = "CD4+ UT", verbose = FALSE)
  582. Idents(obj_combined) <- "cca_N_types_stage"
  583. GBC.de.markers <- FindMarkers(obj_combined, ident.1 = "GBC", ident.2 = NULL, only.pos = TRUE)%>%subset(p_val_adj < 0.05)
  584. e_iOSN.de.markers <- FindMarkers(obj_combined, ident.1 = "e_iOSN", ident.2 = NULL, only.pos = TRUE)%>%subset(p_val_adj < 0.05)
  585. l_iOSN.de.markers <- FindMarkers(obj_combined, ident.1 = "l_iOSN", ident.2 = NULL, only.pos = TRUE)%>%subset(p_val_adj < 0.05)
  586. mOSN.de.markers <- FindMarkers(obj_combined, ident.1 = "mOSN", ident.2 = NULL, only.pos = TRUE)%>%subset(p_val_adj < 0.05)
  588. mOSN.de.markers[1:3,]
  590. # %%
  591. # subset(xx,grepl('apopto',term_name))[1:5,]
  593. # %%
  594. # xx <- gost(rownames(mOSN.de.markers), source='GO:BP' ,
  595. # correction_method = "fdr",user_threshold=0.1,significant=F)$result %>%subset(term_size < 3000)
  596. # xx
  598. # %%
  599. #FeaturePlot(obj, pt.size =0.5, features = rownames(e_iOSN.de.markers)[1:4] ,
  600. # cols = c("#FFF5F0",'#F75D42', "#6A010D") , )+theme(legend.position = "none")
  601. #e_iOSN.de.markers[1:15,]
  603. # %%
  604. #png(filename = "./figures/outermerge_mnc3_pseudotime.png",width = 10, height = 8, units = "cm", res=300)
  605. options(repr.plot.width = 12, repr.plot.height = 12, repr.plot.res = 100)
  607. library(scattermore)#scattermore
  608. UMAPCCA <- as.data.frame(obj@reductions$[email hidden])
  609. UMAPCCA$dataset <- obj$dataset
  610. UMAPCCA$N_types <- obj$N_types
  611. UMAPCCA$N_types_stage <- obj$cca_N_types_stage
  613. brewer.pal(n = 12, name = "Paired")
  615. p1 <- ggplot(UMAPCCA, aes(x=umapcca_1,y=umapcca_2, color=N_types ))+
  616. geom_scattermore(pointsize=5)+theme_bw()+theme(aspect.ratio = 1)+
  617. scale_color_manual(values = c('#1f77b4','#ff7f0e',"#D32B29"),breaks = c('GBC','iOSN','mOSN'))+xlab('UMAP1')+ylab('UMAP2')
  619. p1 <- ggplot(UMAPCCA, aes(x=umapcca_1,y=umapcca_2, color=factor(N_types_stage,levels = c('GBC','e_iOSN','l_iOSN','mOSN') )))+
  620. geom_scattermore(pointsize=5)+theme_bw()+theme(aspect.ratio = 1)+labs(col='Cell types')+
  621. scale_color_manual(values = c('#1f77b4','#FDBF6F','#ff7f0e',"#D32B29"),
  622. breaks = c('GBC','e_iOSN','l_iOSN','mOSN'))+xlab('UMAP1')+ylab('UMAP2')+
  623. theme(panel.grid.major = element_blank(),
  624. panel.grid.minor = element_blank(),
  625. axis.text = element_blank(), # Remove axis text
  626. axis.ticks = element_blank()
  627. )
  629. p2 <- ggplot(UMAPCCA, aes(x=umapcca_1,y=umapcca_2, color=dataset ))+xlab('UMAP1')+ylab('UMAP2')+geom_scattermore(pointsize=5)+theme_bw()+theme(aspect.ratio = 1)+
  630. theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
  631. axis.text = element_blank(), # Remove axis text
  632. axis.ticks = element_blank()
  633. )+scale_color_manual(values = c('#2D81BC','#B8B8B8'))
  635. pdf(file='./figures/03ent_nn_N_merge.umap.pdf',width=8, height=3.5)
  636. p2+p1
  637. dev.off()
  639. # %% [markdown]
  640. # ### slingshot
  642. # %%
  643. #
  644. data_ns <- read_obj('ent_nn_merge_rawcount', hvg = 8000)
  645. data_ns=data_ns[!grepl('^RPS|^RPL|^LINC|^MT',rownames(data_ns)),]
  647. data_ns$batch <- data_ns$Set
  649. #data_ns <- JoinLayers(data_ns)
  651. data_ns <- subset(data_ns, batch%in%c(setdiff(unique(data_ns$batch),c('Set2','Set3','Set4') )))
  652. data_ns <- NormalizeData(data_ns)
  653. data_ns <- FindVariableFeatures(data_ns, selection.method='mean.var.plot')
  654. data_ns@reductions$umap <- obj@reductions$umap.cca
  656. sling_merge <- get_slingshot_traj(obj=obj)
  657. [email hidden]$sling_pseudotime <- sling_merge$sling_pseudotime
  658. [email hidden]$cca_N_types <- [email hidden][rownames([email hidden]),'cca_N_types']
  659. [email hidden]$cca_N_types_stage <- [email hidden][rownames([email hidden]),'cca_N_types_stage']
  660. DEG_sling_merge <- tradeseq_DEG(data_ns=data_ns, 'slingshot')
  663. # %%
  664. get_bg <- function(celltype){
  665. data_sub <- subset(data_ns, cca_N_types_stage=='GBC')
  666. custom_olfbg <- rownames(data_sub)[rowSums(as.data.frame(data_sub@assays$RNA@counts > 0)) > 0.05*ncol(data_sub)]
  667. return(custom_olfbg)
  668. }
  670. custom_olfbg <- lapply(c('GBC','e_iOSN','l_iOSN','mOSN'),get_bg)
  671. custom_olfbg <- unique(unlist(custom_olfbg))
  673. custom_olfbg <- rownames(data_ns)[rowSums(as.data.frame(data_ns@assays$RNA@counts > 0)) > 0.05*ncol(data_ns)]
  674. custom_olfbg <- rownames(data_ns)[rowSums(as.data.frame(data_ns@assays$RNA@counts > 0)) >= 5]
  676. save(custom_olfbg, file='./data/custom_olfbg.RData')
  679. # %% [markdown]
  680. # ### monocle3
  682. # %%
  683. int_N <- get_monocle3_traj(scdata=data_ns,'y','min',res_traj=0.001)
  684. save(int_N, file='./figures/mnc3_ent_nn.RData')
  686. p1 <- plot_cells(int_N, color_cells_by = "pseudotime", cell_size = 0.8,label_roots=F,label_branch_points=F,label_leaves=F)+theme(legend.position = 'top', aspect.ratio = 1)
  687. p2 <- plot_cells(int_N, color_cells_by = "cca_N_types_stage", cell_size = 0.8,label_roots=F,label_branch_points=F,label_leaves=F)+
  688. #scale_color_manual(values = c('#1f77b4','#ff7f0e',"#D32B29"),breaks = c('GBC','iOSN','mOSN'))+
  689. scale_color_manual(values = c('#1f77b4','#FDBF6F','#ff7f0e',"#D32B29"),breaks = c('GBC','e_iOSN','l_iOSN','mOSN'))+
  690. theme(legend.position = 'top', aspect.ratio = 1)
  691. cowplot::plot_grid(p2,p1,nrow = 1)
  693. pdf( "./figures/03outermerge_mnc3_pseudotime_p1.pdf",width = 4, height = 4)
  694. p1
  695. dev.off()
  697. pdf( "./figures/03outermerge_mnc3_pseudotime_p2.pdf",width = 4, height = 4)
  698. p2
  699. dev.off()
  701. # %%
  702. save(int_N, file='./figures/mnc3_ent_nn.RData')
  705. # %%
  708. # %%
  709. data_ns$monocle3_pseudotime <- pseudotime(int_N)
  710. int_N$monocle3_pseudotime <- pseudotime(int_N)
  711. #DEG_mnc3_int_glm <- get_monocle3_DEG_contBatch(int_N)
  712. DEG_mnc3_int_graph <- get_monocle3_DEG(int_N)
  715. # %%
  716. library(ggplot2)
  717. library(ggsignif) # For adding significance tests
  718. library(ggpubr)
  719. N_stage_color <- c('#1f77b4','#FDBF6F','#ff7f0e',"#D32B29")
  720. names(N_stage_color) <- c('GBC','e_iOSN','l_iOSN','mOSN')
  723. [email hidden]$cca_N_types_stage <- factor([email hidden]$cca_N_types_stage, levels = c('GBC','e_iOSN','l_iOSN','mOSN'))
  724. ggboxplot([email hidden], x = "cca_N_types_stage", y = "monocle3_pseudotime",
  725. color = "dataset",
  726. )+stat_compare_means(aes(group = dataset), method = "anova",label = "p.signif")+ylab('Pseudotime') +xlab('')+
  727. scale_color_manual(values = c('#2D81BC','#B8B8B8'))+
  728. theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  730. dev.print(pdf, file='./figures/pseudo_boxplot.pdf',height=3.8,width=3.5)
  732. # %%
  733. write.table([email hidden], file='./data/integrate_entnn_pseduotime_meta.txt',sep='\t',quote=F)
  735. # %%
  736. DEGs1 <- rownames(subset(DEG_mnc3_int_graph, morans_I > 0.05 & q_value < 0.05))
  737. #DEGs2 <- rownames(subset(DEG_mnc3_int_glm, q_value < 0.05))
  738. #DEGs3 <- unique(c(rownames(GBC.de.markers), rownames(e_iOSN.de.markers), rownames(l_iOSN.de.markers), rownames(mOSN.de.markers)))
  741. # %%
  742. save(DEG_mnc3_int_graph,file='./data/DEG_mnc3_int_graph.RData' )
  744. # %%
  745. length(DEGs1)
  746. #length(DEGs2)
  748. # %% [markdown]
  749. # ### Gene clusters
  751. # %%
  753. DEGs <- DEGs1
  754. options(repr.plot.width = 8, repr.plot.height = 8, repr.plot.res = 100)
  756. scdata <- data_ns
  757. save(scdata, file='./data/scdata.RData')
  758. #pseudoTime = scdata$sling_pseudotime#
  759. pseudoTime = scdata$monocle3_pseudotime
  761. pt.matrix1 <- as.matrix(scdata@assays$RNA@counts[DEGs,order(pseudoTime)])
  762. pt.matrix <- t(apply(pt.matrix1,1,function(x){smooth.spline(x,df=3)$y}))
  763. #pt.matrix <- t(apply(pt.matrix,1,function(x){(x-mean(x))/sd(x)}))
  764. pt.matrix <- t(apply(pt.matrix,1,function(x){(x-min(x))/(max(x)-min(x))}))
  766. rownames(pt.matrix) <- DEGs;
  767. colnames(pt.matrix) <- colnames(scdata)[order(pseudoTime)]
  769. n_cluster <- 7
  771. set.seed(123)
  773. heatmap_cluster <- Heatmap(
  774. pt.matrix,
  775. row_km =n_cluster,
  776. col = colorRamp2(seq(from=0,to=1,length=11),rev(brewer.pal(11, "Spectral"))),
  777. #row_title = "cluster_%s",
  778. row_gap = unit(c(2.5), "mm"),
  779. cluster_row_slices = F,
  780. cluster_columns = F,
  781. show_row_names = FALSE,
  782. show_column_names = FALSE,
  783. show_row_dend = F,
  784. top_annotation = HeatmapAnnotation(
  785. cell_identity=[email hidden][colnames(pt.matrix),'cca_N_types'],
  786. cell_subcluster=[email hidden][colnames(pt.matrix),'cca_N_types_stage'],
  787. #batch=[email hidden][colnames(pt.matrix),'batch'],
  788. #subc=[email hidden][colnames(pt.matrix),'N_leiden_res0_2'],
  789. pseudotime=sort(pseudoTime),
  790. col = list(cell_identity = N_color)
  791. # pseudotime= colorRamp2(c(0, 10, 20, 30), c("#35008C", "#8700A8", "#D3546F",'#E7F92D'))
  792. )
  793. )
  794. hthc=ComplexHeatmap::draw(heatmap_cluster)
  795. hthc
  797. save(hthc,file = './figures/heatmap_cluster_graphK7.RData')
  799. # %%
  803. # %%
  804. #load('./figures/heatmap_cluster_graphK7.RData')
  805. DEG_trend_Olf <- c()
  807. down_cl <- c(1,3)
  808. trans_up_cl <- c(6,5,7)
  809. up_cl <- c(4,2)
  810. trans_down_cl <- c()
  812. # down_cl <- c(1,2,4)
  813. # trans_up_cl <- c(7,6)
  814. # up_cl <- c(5,3)
  815. # trans_down_cl <- c()
  817. for(cl in 1:n_cluster){
  818. deg_td <- cbind(DEGs[row_order(hthc)[[cl]]], cl )
  820. DEG_trend_Olf <- rbind(DEG_trend_Olf, deg_td )
  821. }
  822. DEG_trend_Olf <- as.data.frame(DEG_trend_Olf)
  823. colnames(DEG_trend_Olf) <- c('gene','cl')
  826. DEG_trend_Olf <- within(DEG_trend_Olf,{
  827. trend_class=''
  828. trend_class[cl%in%down_cl] <- 'down'
  829. trend_class[cl%in%trans_up_cl] <- 'trans_up'
  830. trend_class[cl%in%up_cl] <- 'up'
  831. trend_class[cl%in%trans_down_cl] <- 'trans_down'
  833. })
  836. order_c <- as.character(c(down_cl,trans_up_cl,up_cl,trans_down_cl ))
  837. #order_c <- as.character(c(6,2,7,1,4,5,3))
  838. order_n <- c(1:length(order_c))
  839. names(order_c) <- order_n
  841. order_c <- sort(order_c)
  842. order_n <- names(order_c)
  844. rownames(DEG_trend_Olf) <- DEG_trend_Olf$gene
  845. DEG_trend_Olf$g_cluster <- as.character(order_n)[as.numeric(DEG_trend_Olf$cl)]
  846. trend_color <- brewer.pal(n = 8, name = "Paired")[c(1,3,5,7)]
  847. names(trend_color) <- c('down','trans_down','up','trans_up')
  849. Nstage_color <- c('#1f77b4','#FDBF6F','#ff7f0e',"#D32B29")
  850. names(Nstage_color) <- c('GBC','e_iOSN','l_iOSN','mOSN')
  854. # %%
  855. library(RColorBrewer)
  856. brewer.pal(n = 8, name = "Paired")[c(1,3,5,7)]
  858. # %%
  859. save(DEG_trend_Olf, file='/sc/arion/projects/roussp01a/liting/Olf/data/DEG_trend_Olf_k7_graph.RData')
  861. # %%
  862. load('/sc/arion/projects/roussp01a/liting/Olf/data/DEG_trend_Olf_k7_graph.RData')
  864. # %%
  865. supp_degs <- subset(DEG_mnc3_int_graph, morans_I > 0.05 & q_value < 0.05)
  866. supp_degs$trajDEG_cluster <- DEG_trend_Olf[rownames(supp_degs),'g_cluster']
  867. supp_degs$trend_class <- DEG_trend_Olf[rownames(supp_degs),'trend_class']
  868. write.table(supp_degs[,c('morans_I','q_value','trajDEG_cluster','trend_class')],file='./data/supp_DEGs.txt',sep='\t')
  870. # %%
  873. # %%
  874. write.table(DEG_trend_Olf, file='./data/DEG_trend_Olf.txt',sep='\t')
  875. table(DEG_trend_Olf$trend_class)
  877. # %%
  878. DEG_trend_Olf[1:3,]
  880. # %%
  881. # get_bp_enrichr <- function(x){
  882. # query_gene <- subset(DEG_trend_Olf, g_cluster%in%as.character(x))[,'gene']
  883. # dbs <- "GO_Biological_Process_2023"
  884. # enriched <- enrichr(query_gene, dbs)
  885. # enrr <- enriched[['GO_Biological_Process_2023']]
  886. # return(enrr)
  887. # }
  890. # enr <- get_bp_enrichr(c(6,7))
  891. # enr[1:35,]
  894. # %%
  895. #load('/sc/arion/projects/roussp01a/liting/Olf/data/DEG_trend_Olf_k7_graph.RData')
  897. lister_hvg <- read.csv('lister_hvg.csv')
  898. #
  899. table(DEG_trend_Olf$gene%in%lister_hvg$X[lister_hvg$highly_variable=='True'])
  900. #DEG_trend_Olf$gene
  902. # %%
  903. library(viridis)
  904. ptm <- pt.matrix[unlist(row_order(hthc)[order_c] ),]
  907. hthc <- ComplexHeatmap::draw(Heatmap(
  908. ptm, name = "Expression",
  909. col = colorRamp2(seq(from=0,to=1,length=9),rev(brewer.pal(9, "Spectral"))),
  910. row_gap = unit(c(2.5), "mm"),
  911. cluster_rows = FALSE,
  912. cluster_columns = F,
  913. show_row_names = FALSE,
  914. show_column_names = FALSE,
  915. use_raster = TRUE, raster_quality = 5,
  916. show_row_dend = F,
  917. row_split = DEG_trend_Olf[rownames(ptm),'g_cluster'],
  918. #heatmap_legend_param = list(legend_width = unit(20, "cm"),title_gap = unit(10, "cm")),
  920. right_annotation = rowAnnotation(Trend = DEG_trend_Olf[rownames(ptm),'trend_class'],
  921. col=list(Trend = trend_color)),
  922. top_annotation = HeatmapAnnotation(
  923. Neuron = [email hidden][colnames(pt.matrix),'cca_N_types'],
  924. Stage = factor([email hidden][colnames(pt.matrix),'cca_N_types_stage'],levels = c('GBC','e_iOSN','l_iOSN','mOSN')),
  925. Pseudotime=sort(pseudoTime),
  926. col = list(Neuron = N_color,
  927. Stage=Nstage_color ,
  928. Pseudotime=colorRamp2(seq(from=0,to=30,length=50), plasma(50)) ),simple_anno_size = unit(0.3, "cm"))),
  929. heatmap_legend_side = "bottom", annotation_legend_side = "bottom",merge_legend = TRUE )
  933. # %%
  934. #png(filename = './figures/03heatmap_DEG_clusters_k7.png', width = 15, height = 20, units = "cm", res=300)
  935. pdf( './figures/03heatmap_DEG_clusters_k7.pdf', width = 5, height = 10)
  936. hthc
  937. dev.off()
  939. # %%
  940. n_cl <- aggregate(DEG_trend_Olf,cl~g_cluster,length)
  943. # %%
  944. #pt.matrix
  946. # %%
  947. library(splines)
  949. # %%
  950. options(repr.plot.width = 8, repr.plot.height = 8, repr.plot.res = 100)
  951. g_mean <- c()
  952. for (i in unique(DEG_trend_Olf$g_cluster)){
  953. Gene <- subset(DEG_trend_Olf,g_cluster==i)[,'gene']
  954. g1_mean <- as.data.frame(cbind(colMeans(pt.matrix[Gene,]), c(sort(pseudoTime)) ))%>%mutate(g_cluster=i)
  955. g_mean <- rbind(g1_mean,g_mean)
  956. }
  958. g_mean <- merge(g_mean, unique(DEG_trend_Olf[,c('trend_class','g_cluster')]), by='g_cluster')
  959. n_cl <- aggregate(DEG_trend_Olf,cl~g_cluster,length)
  960. g_mean_cl <- merge(g_mean, n_cl)
  961. g_mean_cl$gcl <- paste0(g_mean_cl$g_cluster,' ( N = ',g_mean_cl$cl ,')')
  962. #pdf('./figures/03meanexp_DEG_clusters_k7.pdf', width = 4.5, height = 8)
  963. pdf('./figures/03meanexp_DEG_clusters_k7.pdf', width = 3, height = 7)
  965. ggplot(subset(g_mean_cl), aes(V2, V1, color=trend_class)) +
  966. #geom_smooth(method = "loess")+
  967. geom_smooth(method = "lm", formula = y ~ ns(x, df = 3),se = T)+
  968. xlab('Pseudotime')+ylab('Expression')+theme_light() +
  969. theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position='bottom',
  970. legend.title = element_blank())+
  971. facet_wrap(~gcl,nrow = length(order_c)) +
  972. scale_color_manual(values = brewer.pal(n = 8, name = "Paired")[c(1,7,5)])+
  973. scale_y_continuous(breaks = seq(0, 1, by = 1))
  974. dev.off()
  976. # %%
  978. library(splines)
  979. Gene <- subset(DEG_trend_Olf,g_cluster==1)[,'gene']
  981. xx <- reshape2::melt(pt.matrix[Gene,])
  982. xx$pseudoTime <- scdata$monocle3_pseudotime[as.character(xx$Var2)]
  983. xxx <- aggregate(xx,value~Var2+pseudoTime , mean)
  986. ggplot(subset(xx), aes(x=pseudoTime, y=value, group=Var1)) +
  987. #geom_smooth(method = "loess")+
  988. geom_smooth(method = "lm", formula = y ~ ns(x, df = 3),se = F)+
  989. xlab('Pseudotime')+ylab('Expression')+theme_light() +
  990. theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = '')
  994. ggplot(subset(xxx), aes(pseudoTime, value)) +
  995. #geom_smooth(method = "loess")+
  996. geom_smooth(method = "lm", formula = y ~ ns(x, df = 3),se = F)+
  997. xlab('Pseudotime')+ylab('Expression')+theme_light() +
  998. theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = '')
  1000. # %%
  1001. load('/sc/arion/projects/roussp01a/liting/Olf/data/DEG_trend_Olf_k7.RData')
  1003. goterms <- gost(split(DEG_trend_Olf$gene,DEG_trend_Olf$g_cluster ), source='GO:BP' ,#custom_bg = rownames(data) ,
  1004. correction_method = "fdr",user_threshold=0.05,significant=F)$result %>%subset(term_size < 2000)
  1006. sig_goterms <- subset(goterms, p_value < 0.05)#%>%group_by(query)#%>%top_n(20,-p_value)
  1007. sig_goterms <- sig_goterms[,c(1,3:6,9,11)]
  1008. colnames(sig_goterms)[1:2] <- c('trajDEG cluster','FDR')
  1010. write.table(sig_goterms, file='./figures/supple_DE_cluster_GO_all.txt',row.names = F, sep='\t')
  1012. # for(i in 1:7){print(goterms%>%subset(query==i & term_size < 2000 & term_size > 30)%>%top_n(20,-p_value)%>%select(c(query,term_name,p_value)))}
  1015. # dbs <- "GO_Biological_Process_2023"
  1016. # enriched <- lapply(split(DEG_trend_Olf$gene,DEG_trend_Olf$g_cluster ), function(x)enrichr(x, dbs)[['GO_Biological_Process_2023']])
  1017. # enriched

02a.traj_integrate_600-6-40-50-graph_k7_FDR_BG-Cov.ipynb at commit 59edf81, no license · at the source

Overview

Authors: Liting Song1,2,3,4, John F Fullard1,2,3,4, Claire Coleman1,2,3,4, Evelyn Hennigan1,2,3,4, Clara Casey1,2,3,4, Xinyi Wang1,2,3,4, Ayako Kawatake-Kuno1,2,3,4, Stathis Argyriou1,2,3,4, Steven Kleopoulos1,2,3,4, Samuel DeMaria5, Jaroslav Bendl1,2,3,4, Alfred Marc Iloreta6, Pengfei Dong1,2,3,4, Panos Roussos1,2,3,4,7,8
ORCID iDs: Liting Song, John F Fullard, Claire Coleman, Evelyn Hennigan, Xinyi Wang, Samuel DeMaria, Jaroslav Bendl, Pengfei Dong, Panos Roussos
  1. Center for Disease Neurogenomics, Icahn School of Medicine at Mount Sinai, New York, NY USA
  2. Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, NY USA
  3. Department of Psychiatry, Icahn School of Medicine at Mount Sinai, New York, NY USA
  4. Department of Genetics and Genomic Science, Icahn School of Medicine at Mount Sinai, New York, NY USA
  5. Department of Anesthesiology, Icahn School of Medicine at Mount Sinai, New York, NY USA
  6. Department of Otolaryngology, Icahn School of Medicine at Mount Sinai, New York, NY USA
  7. Mental Illness Research Education and Clinical Center (MIRECC), James J. Peters VA Medical Center, Bronx, New York, NY USA
  8. Center for Precision Medicine and Translational Therapeutics, James J. Peters VA Medical Center, Bronx, New York, NY USA
Institutions: Icahn School of Medicine at Mount Sinai (United States); James J. Peters VA Medical Center (United States); Mental Illness Research, Education and Clinical Centers (United States)
Journal: Nature communications, volume 17, issue 1, article 8744
Dates: received 6 August 2025; accepted 9 July 2026; published online 17 July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-75722-1 · PMID 42463696 · PMCID PMC13494014 · OpenAlex W7169028888
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), developmental (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials
Keywords: Olfactory system, Genetics of the nervous system
MeSH: Brain*, Neurogenesis*, Olfactory Mucosa*, Olfactory Receptor Neurons*, Adult, Cell Differentiation, Cell Nucleus, Female, Gene Expression Profiling, Humans, Male, Neural Stem Cells, Neurodevelopment, RNA-Seq, Single-Cell Analysis, Single-Cell Gene Expression Analysis, Transcription Factors, Transcriptome (* major topic)
Topic: Olfactory and Sensory Function Studies (Sensory Systems, Neuroscience), according to OpenAlex
Funding: NIA NIH HHS (R01 AG067025, R01 AG050986, R01 AG065582, R01 AG082185, U24 AG087563); U.S. Department of Health & Human Services | NIH | National Institute on Aging (U.S. National Institute on Aging) (R01-AG067025, R01-AG065582, R01-AG082185, U24-AG087563, R01-AG050986); NIMH NIH HHS (R01 MH110921, RF1 MH133703, R01 MH125246, U01 MH116442); U.S. Department of Health & Human Services | NIH | National Institute of Mental Health (NIMH) (R01-MH110921, R01-MH125246, RF1-MH133703, U01-MH116442); NINDS NIH HHS (U01 NS125580); U.S. Department of Health & Human Services | NIH | National Institute of Neurological Disorders and Stroke (NINDS) (U01NS125580)
Citations: cited by 1 paper (Europe PMC); 93 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

well.ox.ac.uk/~wrayner/tools

License: none: the authors keep all their rights
State: the link is dead, verified on 27 September 2026
Evidence: found in the paper
Software Heritage: not checked
Found in: the text, “SNP array genotyping workflow”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link is dead (HTTP 404)
  • 27 September 2026: the link is dead (HTTP 404)
At the source: well.ox.ac.uk/~wrayner/tools/

Zenodo 20867452

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: anndata (6 files), pandas (6 files), Scanpy (6 files), seaborn (4 files), tidyverse (4 files), ggplot2 (3 files), Matplotlib (3 files), cowplot (2 files), ggpubr (2 files), NumPy (2 files), Seurat (2 files), circlize (1 file), ComplexHeatmap (1 file), Monocle 3 (1 file), reshape2 (1 file), scikit-learn (1 file), SciPy (1 file), SingleCellExperiment (1 file), UMAP (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
10 files
At the source: zenodo.org/records/20867452

Zenodo 20867451

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: anndata (6 files), pandas (6 files), Scanpy (6 files), seaborn (4 files), tidyverse (4 files), ggplot2 (3 files), Matplotlib (3 files), cowplot (2 files), ggpubr (2 files), NumPy (2 files), Seurat (2 files), circlize (1 file), ComplexHeatmap (1 file), Monocle 3 (1 file), reshape2 (1 file), scikit-learn (1 file), SciPy (1 file), SingleCellExperiment (1 file), UMAP (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
10 files
At the source:

litingsong/oe_nn

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 59edf81e4b2c1655b39e1afa1873dd6968f2ea07, 26 June 2026
Languages: Jupyter (21), R (5), Python (2)
Size: 29 files, 28 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, 21 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: anndata (6 files), pandas (6 files), Scanpy (6 files), seaborn (4 files), tidyverse (4 files), ggplot2 (3 files), Matplotlib (3 files), cowplot (2 files), ggpubr (2 files), NumPy (2 files), Seurat (2 files), circlize (1 file), ComplexHeatmap (1 file), Monocle 3 (1 file), reshape2 (1 file), scikit-learn (1 file), SciPy (1 file), SingleCellExperiment (1 file), UMAP (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
10 files
At the source: github.com/LitingSong/OE_NN/tree/v1.0.0

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-75722-1.

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:

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

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

Data

Datasets cited

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-75722-1.

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, 14 authors, 2 keywords, 18 MeSH terms, 6 funders, 90 references.

Cite

This paper

Song, L., Fullard, J. F., Coleman, C., Hennigan, E., Casey, C., Wang, X., Kawatake-Kuno, A., Argyriou, S., Kleopoulos, S., DeMaria, S., Bendl, J., Iloreta, A. M., Dong, P., & Roussos, P. (2026). Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain. Nature communications, 17(1), 8744. https://doi.org/10.1038/s41467-026-75722-1

BibTeX

@article{song2026single,
author = {Song, Liting and Fullard, John F and Coleman, Claire and Hennigan, Evelyn and Casey, Clara and Wang, Xinyi and Kawatake-Kuno, Ayako and Argyriou, Stathis and Kleopoulos, Steven and DeMaria, Samuel and Bendl, Jaroslav and Iloreta, Alfred Marc and Dong, Pengfei and Roussos, Panos},
title = {{Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {8744},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75722-1},
url = {https://doi.org/10.1038/s41467-026-75722-1},
pmid = {42463696},
pmcid = {PMC13494014}
}

RIS

TY - JOUR
AU - Song, Liting
AU - Fullard, John F
AU - Coleman, Claire
AU - Hennigan, Evelyn
AU - Casey, Clara
AU - Wang, Xinyi
AU - Kawatake-Kuno, Ayako
AU - Argyriou, Stathis
AU - Kleopoulos, Steven
AU - DeMaria, Samuel
AU - Bendl, Jaroslav
AU - Iloreta, Alfred Marc
AU - Dong, Pengfei
AU - Roussos, Panos
TI - Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/17
VL - 17
IS - 1
SP - 8744
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75722-1
UR - https://doi.org/10.1038/s41467-026-75722-1
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75722-1",
"type": "article-journal",
"title": "Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain",
"container-title": "Nature communications",
"author": [
{
"family": "Song",
"given": "Liting"
},
{
"family": "Fullard",
"given": "John F"
},
{
"family": "Coleman",
"given": "Claire"
},
{
"family": "Hennigan",
"given": "Evelyn"
},
{
"family": "Casey",
"given": "Clara"
},
{
"family": "Wang",
"given": "Xinyi"
},
{
"family": "Kawatake-Kuno",
"given": "Ayako"
},
{
"family": "Argyriou",
"given": "Stathis"
},
{
"family": "Kleopoulos",
"given": "Steven"
},
{
"family": "DeMaria",
"given": "Samuel"
},
{
"family": "Bendl",
"given": "Jaroslav"
},
{
"family": "Iloreta",
"given": "Alfred Marc"
},
{
"family": "Dong",
"given": "Pengfei"
},
{
"family": "Roussos",
"given": "Panos"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8744",
"DOI": "10.1038/s41467-026-75722-1",
"PMID": "42463696",
"PMCID": "PMC13494014",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75722-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
17
]
]
}
}

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: Monocle 3, SingleCellExperiment, UMAP, 16 other tools, genetics / omics, 4 references
[2] doi:10.1038/s41467-026-71595-6 [code]
A single-cell and spatial atlas of early human olfactory development.
Journal: Nature communications
In common: SingleCellExperiment, anndata, circlize, 12 other tools, genetics / omics, 6 references
[3] doi:10.1038/s41467-026-71803-3 [code]
Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: anndata, Scanpy, Seurat, 10 other tools, genetics / omics, 9 references
[4] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: Monocle 3, UMAP, anndata, 14 other tools, genetics / omics, 3 references
[5] doi:10.1038/s41467-026-71525-6 [code]
Single-nucleus brain transcriptomics reveals microglia dysfunction in multiple system atrophy.
Journal: Nature communications
In common: UMAP, anndata, circlize, 11 other tools, genetics / omics, 5 references
[6] doi:10.1038/s42003-026-10034-0 [code]
Region- and cell type-specific changes in gene expression in the cerebellum after classical fear conditioning.
Journal: Communications biology
In common: SingleCellExperiment, anndata, circlize, 13 other tools, genetics / omics, 3 references
[7] doi:10.1038/s41586-026-10490-y [code]
Lineage and organ signals sequentially build organ intrinsic nervous systems.
Journal: Nature
In common: Monocle 3, UMAP, anndata, 13 other tools, developmental, 3 references
[8] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Monocle 3, SingleCellExperiment, anndata, 15 other tools, genetics / omics
[9] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Monocle 3, SingleCellExperiment, anndata, 15 other tools, genetics / omics
[10] doi:10.1038/s41514-026-00391-9 [code]
Region-specific transcriptional signatures of brain aging in the absence of neuropathology at the single-cell level.
Journal: npj aging
In common: anndata, circlize, Scanpy, 13 other tools, genetics / omics, 3 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.