add figures

run.sh

@@ -2,4 +2,4 @@
 
 # R -e "devtools::document()"
 # R -e "pkgdown::build_site(preview = FALSE, lazy=TRUE)"
-# R -e "pkgdown::build_article(name='myofibroblast-GRN', quiet=FALSE)"
+R -e "pkgdown::build_article(name='myofibroblast-GRN', quiet=FALSE)"

vignettes/myofibroblast-GRN.Rmd

@@ -23,22 +23,22 @@ The script of cleaning the data and preparing these objects is found
 [here](https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/01_prepare_data.html).
 
 Run the following commands to download the data:
-```{r download_data}
-# if(!dir.exists('./Fibroblast')){
-#     dir.create('./Fibroblast')
-# }
-
-# download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/snRNA.rds', 
-#               destfile = './Fibroblast/snRNA.rds', 
-#               method = 'wget', extra = '--no-check-certificate')
-
-# download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/snATAC.rds', 
-#               destfile = './Fibroblast/snATAC.rds', 
-#               method = 'wget', extra = '--no-check-certificate')
+```{r download_data, eval=TRUE}
+if(!dir.exists('./Fibroblast')){
+    dir.create('./Fibroblast')
+}
+
+download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/snRNA.rds', 
+              destfile = './Fibroblast/snRNA.rds', 
+              method = 'wget', extra = '--no-check-certificate')
+
+download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/snATAC.rds', 
+              destfile = './Fibroblast/snATAC.rds', 
+              method = 'wget', extra = '--no-check-certificate')
               
-# download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/gene.activity.rds', 
-#               destfile = './Fibroblast/gene.activity.rds', 
-#               method = 'wget', extra = '--no-check-certificate')
+download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/Fibroblast/gene.activity.rds', 
+              destfile = './Fibroblast/gene.activity.rds', 
+              method = 'wget', extra = '--no-check-certificate')
 ```
 
 
@@ -555,123 +555,14 @@ p2 <- PseudotimePlot(object = obj, tf.use = "RUNX1")
 p1 + p2
 ```
 
-The x-axis in above plots present pseudotime point along the trajectory, and the y-axis represent TF binding acitivty, TF expression, and TF target expression after z-score transformation.
-
-### GRN in space
-
-As an additional step, if the spatial transcriptome (ST) data (e.g., generated by using 10X Visium assay) are available, we can then visualize the target gene expression for some interesting TFs to understand the regulon activity in spatial space. First let's download the ST data
-
-```{r, eval=FALSE}
-dir.create('./VisiumSpatial')
-
-download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/VisiumSpatial/AKK004_157772.rds', 
-              destfile = './VisiumSpatial/AKK004_157772.rds', 
-              method = 'wget', extra = '--no-check-certificate')
-```
-
-Of note, each spot in the 10X Visium assay contains multiple cells, meaning that we usually cannot identify cell types by clustering the data. To address this issue, people developed various computation tools to estimate cell-type-specific abundance in each spot. Here we have perform such analysis using [cell2location](https://www.nature.com/articles/s41587-021-01139-4) based on our snRNA-seq data. Let's first visualize the fibroblast proportion.
-
-
-```{r, fig.height = 5, fig.width = 5, fig.align = "center", eval=TRUE}
-obj.spatial <- readRDS("./VisiumSpatial/AKK004_157772.rds")
-
-DefaultAssay(obj.spatial) <- "c2l_props"
-
-p <- SpatialFeaturePlot(obj.spatial, features = "Fib", min.cutoff = "q5",
-                       max.cutoff = "q95", image.alpha = 1, 
-                       pt.size.factor = 400, keep.scale = NULL) +
-    scale_fill_viridis(option = "D") +
-    ggtitle("") +
-    labs(fill = "") +
-    theme(legend.position = "bottom")
-
-print(p)
-```
-
-We next can check the gene expression in space. Here, we use the above TFs.
-
-
-```{r, fig.height = 5, fig.width = 10, fig.align = "center", eval=TRUE}
-DefaultAssay(obj.spatial) <- "SCT"
-
-p1 <- SpatialFeaturePlot(object = obj.spatial, features = "NR3C2", image.alpha = 1, pt.size.factor = 400, keep.scale = NULL) +
-    ggtitle("") +
-    theme(legend.position = "bottom") +
-scale_fill_viridis(option = "C")
-
-p2 <- SpatialFeaturePlot(object = obj.spatial, features = "RUNX1", image.alpha = 1, pt.size.factor = 400, keep.scale = NULL) +
-    ggtitle("") +
-    theme(legend.position = "bottom") +
-scale_fill_viridis(option = "C")
-
-p1 + p2
-```
-
-Then, we visualize the regulon activity of NR3C2 and RUNX1 by using the GRN that we predicted. The major idea is that if the target genes of a TF are up-regulated, then this TF should have high activity. This principle has been used to develop computation methods to estimate TF activity from gene expression data, such as [DoRothEA](https://genome.cshlp.org/content/29/8/1363.short).
-
-```{r, fig.height = 5, fig.width = 10, fig.align = "center", eval=TRUE}
-DefaultAssay(obj.spatial) <- "SCT"
-
-p1 <- GRNSpatialPlot(object = obj.spatial, assay = "SCT", df.grn = df.grn2,
-                   tf.use = "NR3C2", vis.option = "B") +
-    ggtitle("") +
-    theme(legend.position = "bottom") 
-
-p2 <- GRNSpatialPlot(object = obj.spatial, assay = "SCT", df.grn = df.grn2,
-                   tf.use = "RUNX1", vis.option = "B") +
-    ggtitle("") +
-    theme(legend.position = "bottom")
-
-p1 + p2
-```
-
-The above plots demonstrate a clear spatial-distributed pattern for the regulon activity of NR3C2 and RUNX1.
+The x-axis in above plots present pseudotime point along the trajectory, 
+and the y-axis represent TF binding acitivty, TF expression, and TF target expression after z-score transformation.
 
 Save data
 ```{r, fig.height = 5, fig.width = 6, fig.align = "center", eval=TRUE}
 saveRDS(obj, "./Fibroblast/coembed.trajectory.rds")
 ```
 
-## Network analysis
-
-As an addition step, we next will explore the gene regulatory networks generated by scMEGA in fibroblasts during myocardial infarction by using network analysis methods. 
-
-This will allow us to 
-
-* explore regulators regarding their network topological properties
-* to vizualise regulatory networks centred in selected transcription factors.  
-
-We first load the network produced by scMEGA:
-```{r, fig.height = 5, fig.width = 6, fig.align = "center", eval=TRUE}
-dfgrn <- readRDS("./Fibroblast/final_grn.Rds")
-netobj <- graph_from_data_frame(dfgrn,directed = TRUE)
-V(netobj)$type <- ifelse(V(netobj)$name %in% dfgrn$tf,"TF/Gene","Gene")
-```
-
-Next, we will estimate distint network topological measures as node degree, node centrality and node page rank score and plot these in a principal component based embedding with a focus on transcription factors. We plot bellow the 3 first PCs. 
-
-### Network topological measures embedding: PC1 and PC2
-```{r,fig.width=12, fig.height=12, eval=TRUE}
-p <- TopEmbGRN(df.grn=netobj)
-```
-
-### Network topological measures embedding: PC2 and PC3
-```{r,fig.width=12, fig.height=12, eval=TRUE}
-p <- TopEmbGRN(df.grn=netobj,axis=c(2,3))
-```
-
-This reveals that RUNX1 is a central TF in this network, as it has high in/out-degree values and high mediator score. 
-
-### RUNX1 centric network. 
-
-To further inspect genes and TFs associated with RUNX1, we create a centric layout embedding.
-
-```{r,fig.width=12,fig.height=12, eval=TRUE}
-NetCentPlot(netobj,"RUNX1")
-```
-
-This network reveals direct targets of RUNX1 such as TEAD4 and HIF1A and indirect targets as TEAD1 and TEAD2. If you want, you can explore a highlight parameter to personalize the plot by only showing the desired genes/TFs. By default, the TFs with a betweenness score higher than an average is depicted in the plot.
-
 
 ```{r}
 # Check session information