3-R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "16cd4b0e",
   "metadata": {},
   "source": [
    "## R: Differential expression analysis, Monocle/trajectory analysis, cluster 6 analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "007232f4",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(dplyr)\n",
    "library(Seurat)\n",
    "library(patchwork)\n",
    "library(tidyverse)\n",
    "library(SeuratWrappers)\n",
    "library(monocle3)\n",
    "library(biomaRt)\n",
    "library(ggrepel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb785bf",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "data = readRDS(\"seurat_final.rds\")\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57695927",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "DimPlot(data, reduction = \"umap\", group.by =\"seurat_clusters\", pt.size = 0.1, cols =\n",
    "          c('0' = '#0078B9', '1' = '#FF7C00', '2'='#00A163', '3'='#E7001A', '4'='#B62BFE', '5'='#92563A', \n",
    "            '6'='#F56DC2','7'='#ADC152', '8'='#00C1D3', '9'='#A6C9ED', '10'='#FFB86B', '11'='#81E580', \n",
    "            '12'='#FF9191')) + \n",
    "  ggtitle(\"\") + xlab(\"UMAP 1\") + ylab(\"UMAP 2\")\n",
    "table(data@meta.data$patient)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ada909fa",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "DimPlot(data, reduction = \"umap\", group.by =\"type\", pt.size = 0.1, \n",
    "        cols = c('Neonate'='#FF7300','Child'='#006EB6')) + ggtitle(\"\") + \n",
    "  xlab(\"UMAP 1\") + ylab(\"UMAP 2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a4616af",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "counts <- group_by(data@meta.data, type, seurat_clusters) %>% summarise(count = n())\n",
    "counts$count_norm <- ave(counts$count, counts$type, FUN=function(x) x/sum(x))\n",
    "ggplot(counts, aes(type, count_norm, fill = seurat_clusters)) + \n",
    "                         geom_bar(stat = 'identity', colour=\"black\", size = 0.1) + \n",
    "                         theme_classic() + theme(text = element_text(size = 16),legend.position = \"none\") + \n",
    "                         xlab(\"\") + ylab(\"Fraction of Cells\") +\n",
    "                         scale_fill_manual(\"legend\" ,values = c('0' = '#0078B9', '1' = '#FF7C00', '2'='#00A163', \n",
    "                                                                '3'='#E7001A', '4'='#B62BFE', '5'='#92563A', \n",
    "                                                                '6'='#F56DC2', '7'='#ADC152', '8'='#00C1D3', \n",
    "                                                                '9'='#A6C9ED', '10'='#FFB86B', '11'='#81E580', \n",
    "                                                                '12'='#FF9191'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00fd4070",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "DotPlot(data, features = c(\"PDGFRA\",\"DDR2\",\"COL1A1\",\"POSTN\",\n",
    "                            \"ACKR3\",\"THY1\",\"AXL\",\"CD34\",\n",
    "                            \"APOE\",\"BMP4\",\"ADM\",\n",
    "                            \"WIF1\",\"WNT5A\",\"WNT16\",\"DKK3\",\"SFRP2\",\"FRZB\",\"AXIN2\",\"PTN\",\"TGFB3\",\"CCL19\",\"CXCL12\",\n",
    "                            \"ELN\",\"MFAP4\",\"COL3A1\",\"COL14A1\",\"FBLN1\",\"NPPC\",\"CTGF\",\"NOX4\",\"ACTA2\",\"PLOD2\"\n",
    "                         )) + \n",
    "  theme(axis.text.x = element_text(angle = 90,hjust=1)) + coord_flip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9339a787",
   "metadata": {
    "collapsed": true,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c0cd404",
   "metadata": {},
   "source": [
    "### Get cluster markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34c23950",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "markers = FindAllMarkers(data, logfc.threshold = 0.25, \n",
    "                        test.use = \"wilcox\", \n",
    "                        min.pct = 0.1, min.diff.pct = 0.1,\n",
    "                        only.pos = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b6bad41",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "DotPlot(obj, features = c(\"ACTA2\",\"IGFBP3\",\"IGFBP5\",\n",
    "                         \"SOX17\",\"TCF4\",\"SOX4\",\n",
    "                         \"\")) + \n",
    "  theme(axis.text.x = element_text(angle = 90,hjust=1)) + coord_flip()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f9135c4",
   "metadata": {},
   "source": [
    "# Bulk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d03733b",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "data@meta.data %>% head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b2ef7d2",
   "metadata": {},
   "source": [
    "### PCA neonate vs. child"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b2e411",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "sub_agg = AggregateExpression(data, slot = \"scale.data\", group.by = \"patient\", verbose = TRUE)\n",
    "# png(\"PCApseudobulk.png\", width = 2200, height = 2000, res = 300)\n",
    "a = prcomp(sub_agg$SCT)$rotation %>% as.data.frame() %>% rownames_to_column(\"patient\") \n",
    "right_join(a, data@meta.data[,c(3,4)] %>% group_by(patient,type) %>% summarise(count = n())) %>% \n",
    "    ggplot(aes(x = PC1, y = PC2, label = patient)) + \n",
    "        geom_point(aes(color = type), size = 5) + theme_classic() +\n",
    "        scale_color_manual(values = c('Neonate' = '#F47C20', 'Child' = '#2E7BBD')) + \n",
    "        theme(text = element_text(size = 20))       \n",
    "# dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2d10cab",
   "metadata": {},
   "source": [
    "## nCPC vs cCPC each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8c35523",
   "metadata": {
    "collapsed": true,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "rna2 = data\n",
    "rna2$edgeR = paste0(rna2@meta.data$patient, \"_\", rna2@meta.data$seurat_clusters)\n",
    "columns = as.character(levels(factor(rna2@meta.data$edgeR)))\n",
    "columns[1:5]\n",
    "length(columns)\n",
    "\n",
    "var = rna2\n",
    "Idents(var) = \"edgeR\"\n",
    "countmatrix = data.frame(rep(0,36601))\n",
    "for (i in 1:length(columns)){\n",
    "    print(columns[i])\n",
    "    var_subset = subset(var, idents = columns[i])\n",
    "    print(dim(var_subset))\n",
    "    var_umi = GetAssayData(var_subset, slot = \"counts\")\n",
    "    var_umi_mat = as.matrix(var_umi)\n",
    "    var_sums = rowSums(var_umi_mat)\n",
    "    countmatrix = cbind(countmatrix, var_sums)}\n",
    "countmtx = countmatrix[,2:ncol(countmatrix)] \n",
    "colnames(countmtx) = columns\n",
    "countmtx%>% tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab2ae658",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "meta = rna2@meta.data[,3:6] %>% group_by_all() %>% tally()\n",
    "meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15d40db8",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "p = str_extract(colnames(countmtx), \"[^_]+\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6262c9ad",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "patient = str_extract(colnames(countmtx), \"[^_]+\")\n",
    "passage = c(); for(j in 1:length(patient)){passage = c(passage,unique(meta[meta$passage == patient[j],3])[[1]])}\n",
    "patient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "856b1fd1",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "left_join(as.data.frame(patient),meta) %>% head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ace239b5",
   "metadata": {
    "collapsed": true,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(edgeR)\n",
    "for(i in 0:12){\n",
    "    cluster = i \n",
    "    pseudo = countmtx[,word(colnames(countmtx), 2, sep = \"_\") == cluster]\n",
    "    patient = str_extract(colnames(pseudo), \"[^_]+\")\n",
    "    passage = left_join(as.data.frame(patient),meta)$passage\n",
    "    group = left_join(as.data.frame(patient),meta)$type\n",
    "    #disease = left_join(as.data.frame(patient),meta)$disease\n",
    "    y = DGEList(counts = pseudo, group = group)\n",
    "    design = model.matrix(~passage+group)\n",
    "    rownames(design) = colnames(y)\n",
    "    keep = filterByExpr(y)\n",
    "    DGEy = y[keep, , keep.lib.sizes=FALSE]\n",
    "    table(keep)\n",
    "    y = calcNormFactors(y)\n",
    "    y = estimateDisp(y, design)\n",
    "    fitlr = glmFit(y, design) #Fit negative binomial generalized log-linear model to read counts for each gene \n",
    "    lr =  glmLRT(fitlr) #\n",
    "    fdr = p.adjust(lr$table$PValue,method=\"BH\")\n",
    "    tab = topTags(lr,n = nrow(lr))\n",
    "    write.csv(tab, paste0(\"cluster\",i,\"_DEGs.csv\"))\n",
    "    \n",
    "    tab$table = tab$table %>% rownames_to_column(\"hgnc_symbol\")\n",
    "    tab$table$p = as.factor(ifelse(tab$table$FDR < 0.05, \"p\", \"no\"))\n",
    "    tab$table$dir = as.factor(ifelse((tab$table$FDR < 0.05 & abs(tab$table$logFC)>1), \n",
    "                                       if_else(tab$table$logFC > 1, \"Neonate\",\"Child\"), \"NS\"))\n",
    "    print(table(tab$table$dir))\n",
    "    tab$table$RNA = \"\"\n",
    "    #get protein_coding only to label\n",
    "    mart = useMart(\"ENSEMBL_MART_ENSEMBL\", host = \"useast.ensembl.org\")\n",
    "    mart = useDataset(\"hsapiens_gene_ensembl\", mart)\n",
    "    pc = left_join(tab$table, getBM(mart = mart,attributes = c(\"hgnc_symbol\",\"gene_biotype\"),\n",
    "                                            filter = \"hgnc_symbol\",\n",
    "                                            values = tab$table$hgnc_symbol,uniqueRows=TRUE)) %>% \n",
    "                                                filter(gene_biotype == \"protein_coding\")\n",
    "    toptags = pc %>% group_by(dir) %>% filter(dir!=\"NS\") %>% \n",
    "                slice_max(n = 10, order_by = -log10(FDR))\n",
    "    #plot\n",
    "    tab$table$RNA[which(tab$table$hgnc_symbol %in% toptags$hgnc_symbol)] = tab$table$hgnc_symbol[tab$table$hgnc_symbol %in% toptags$hgnc_symbol]\n",
    "    \n",
    "    assign(paste0(\"cluster\",i,\"p\"),tab$table)\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caf3133a",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "head(cluster0p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b93886cb",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "png(\"cluster7.png\",res = 200, height = 1500, width = 1500)\n",
    "cluster7p %>% ggplot(aes(x = logFC, y = -log10(FDR), label = RNA)) + \n",
    "        geom_point(aes(color = dir)) +\n",
    "        theme_classic() + theme(text = element_text(size = 14),legend.title=element_blank()) +\n",
    "        geom_hline(yintercept = -log10(0.05), color = \"grey\") + \n",
    "        geom_vline(xintercept = c(-1,1), color = \"grey\") +\n",
    "        scale_color_manual(values = c(\"Neonate\" = \"#F47C20\", 'Child' = '#2E7BBD', \"NS\" = \"black\"))+\n",
    "        geom_text_repel()+ggtitle(\"Cluster 7\")+labs(caption = \"Model: ~ passage + age group\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e1a9b7b",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "c6 = subset(data, cells = colnames(data)[Idents(data)==6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c66543a9",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "dim(c6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4b75314",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#c6 = RunPCA(c6, verbose = FALSE)\n",
    "#c6 = RunUMAP(c6, dims = 1:30, verbose = FALSE)\n",
    "#c6 = FindNeighbors(c6, dims = 1:30, verbose = FALSE)\n",
    "c6 = FindClusters(c6, resolution = 0.05, verbose = FALSE)\n",
    "DimPlot(c6, reduction = \"umap\", label = TRUE)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6cfe73",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#png(paste0(maindir,\"plots/Dotplot_Clusters_3.png\"), width = 1800, height = 2200, res = 300)\n",
    "DoHeatmap(c6 %>% FindVariableFeatures(nfeatures = 50), label = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7edbc923",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "DotPlot(c6,features = c(\"COL1A1\",\"PDGFRA\",\"THY1\",\n",
    "             \"POSTN\",\"PLOD2\",\"ACTA2\",\n",
    "            \"WIF1\",\"WNT5A\",\"DKK3\",\"CXCL12\",# Wntx\n",
    "            \"COL3A1\",\"COL14A1\",\"FBLN1\"       # F-trans\n",
    "            \n",
    "                         )) + \n",
    "  theme(axis.text.x = element_text(angle = 90,hjust=1)) + coord_flip()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4a25fef",
   "metadata": {},
   "source": [
    "# Trajectory Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fbb59fa",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#trace(\"calculateLW\", edit = T, where = asNamespace(\"monocle3\"))\n",
    "cds <- as.cell_data_set(data)\n",
    "cds <- estimate_size_factors(cds)\n",
    "rowData(cds)$gene_name <- rownames(cds)\n",
    "rowData(cds)$gene_short_name <- rowData(cds)$gene_name\n",
    "cds <- cluster_cells(cds = cds, reduction_method = \"UMAP\")\n",
    "cds <- learn_graph(cds, learn_graph_control=list(ncenter=500))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "801ed0c5",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87cf245f",
   "metadata": {},
   "source": [
    "# Coexpressed genes along trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37c71d52",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "trace(\"calculateLW\", edit = T, where = asNamespace(\"monocle3\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ecc9905",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "de_res <- graph_test(cds, neighbor_graph = \"principal_graph\", cores = 48)\n",
    "#saveRDS(de_res, \"de_res_09082022_revisions.rds\")\n",
    "head(de_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "987ff3dc",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(data.table)\n",
    "de_res_filt <- subset(de_res, q_value < 0.05)\n",
    "de_res_dt <- data.table(de_res_filt)\n",
    "plot_cells(cds, genes=as.matrix(de_res_dt[order(-morans_I),][1:20,\"gene_name\"])[,1],\n",
    "           show_trajectory_graph=FALSE,\n",
    "           label_cell_groups=FALSE,\n",
    "           label_leaves=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6e5e8c3",
   "metadata": {},
   "source": [
    "# Use Monocle's built-in batch correction before finding gene modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13f7938d",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cds <- preprocess_cds(cds)\n",
    "cds <- align_cds(cds, alignment_group = \"batch\")\n",
    "gene_module_df <- find_gene_modules(cds[rownames(de_res_filt),], resolution = c(10^seq(-6,-1)), cores = 48)b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f110c3c",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(RColorBrewer)\n",
    "cell_group_df <- tibble::tibble(cell=row.names(colData(cds)), \n",
    "                                cell_group=colData(cds)$seurat_clusters)\n",
    "agg_mat <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)\n",
    "row.names(agg_mat) <- stringr::str_c(\"Module \", row.names(agg_mat))\n",
    "pheatmap::pheatmap(agg_mat,\n",
    "                   scale=\"column\", clustering_method=\"ward.D2\")b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "858c2ffd",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(ComplexHeatmap)\n",
    "clust_mod <- agg_mat\n",
    "cell_group_df <- tibble::tibble(cell=row.names(colData(cds)), \n",
    "                                cell_group=colData(cds)$type)\n",
    "age_mod <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)\n",
    "row.names(age_mod) <- stringr::str_c(\"Module \", row.names(age_mod))\n",
    "\n",
    "scaled_mat = scale(as.matrix(clust_mod))\n",
    "\n",
    "Child <- circlize::colorRamp2(seq(min(age_mod[,1]), max(age_mod[,1]), length = 7), brewer.pal(n = 7, name = \"Blues\"))\n",
    "Neonate <- circlize::colorRamp2(seq(min(age_mod[,2]), max(age_mod[,2]), length = 7), brewer.pal(n = 7, name = \"Oranges\"))\n",
    "\n",
    "colnames(scaled_mat) <- stringr::str_c(\"Cluster \", colnames(scaled_mat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0c04571",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "color = list(Child = Child,\n",
    "            Neonate = Neonate)\n",
    "row_ha = rowAnnotation(Child = age_mod[,1], Neonate = age_mod[,2], col=color)\n",
    "col_ha = HeatmapAnnotation(Child = anno_barplot(counts[1:13,4], ylim = c(0,0.35)), Neonate = anno_barplot(counts[14:26,4], ylim = c(0,0.35)))\n",
    "Heatmap(scaled_mat, \n",
    "        right_annotation = row_ha, \n",
    "        top_annotation = col_ha, \n",
    "        col=colorRampPalette(rev(brewer.pal(n = 7, name = \"RdYlBu\")))(100),\n",
    "        name=\"Module\\nExpression\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}