diff --git a/.github/workflows/render-rmarkdown.yaml b/.github/workflows/render-rmarkdown.yaml
index 3a4ad81..5621cb4 100644
--- a/.github/workflows/render-rmarkdown.yaml
+++ b/.github/workflows/render-rmarkdown.yaml
@@ -45,9 +45,13 @@ jobs:
         run: |
           devtools::install_github("CostaLab/scMEGA")
         shell: Rscript {0}
+      # - name: Render Rmarkdown files
+      #   run: |
+      #     pkgdown::build_site()
+      #   shell: Rscript {0}
       - name: Render Rmarkdown files
-        run: |
-          pkgdown::build_site()
+        run : |
+          pkgdown::build_article(name='myofibroblast-GRN', quiet=FALSE)
         shell: Rscript {0}
       - name: Commit files
         run: |
diff --git a/vignettes/myofibroblast-GRN.Rmd b/vignettes/myofibroblast-GRN.Rmd
index a867e46..c281afe 100644
--- a/vignettes/myofibroblast-GRN.Rmd
+++ b/vignettes/myofibroblast-GRN.Rmd
@@ -483,7 +483,9 @@ tf.gene.cor <- GetTFGeneCorrelation(object = obj,
                                     trajectory.name = "Trajectory")
 ```
 
-We can then visualize this correlation matrix by heatmap. Also, we can cluster the genes and TFs to identify different regulatory modules for the predefined sub-populations.
+We can then visualize this correlation matrix by heatmap. 
+Also, we can cluster the genes and TFs to identify different regulatory modules for 
+the predefined sub-populations.
 
 ```{r, fig.height = 6, fig.width = 12, fig.align = "center"}
 ht <- GRNHeatmap(tf.gene.cor, 
@@ -491,7 +493,9 @@ ht <- GRNHeatmap(tf.gene.cor,
 ht
 ```
 
-To associate genes to TFs, we will use the peak-to-gene links and TF binding sites information. Specifically, if a gene is regulated by a peak and this peak is bound by a TF, then we consider this gene to be a target of this TF.
+To associate genes to TFs, we will use the peak-to-gene links and TF binding sites information. 
+Specifically, if a gene is regulated by a peak and this peak is bound by a TF, then we consider 
+this gene to be a target of this TF.
 
 ```{r, fig.height = 6, fig.width = 12, fig.align = "center"}
 motif.matching <- obj@assays$ATAC@motifs@data
@@ -508,7 +512,7 @@ df.grn <- GetGRN(motif.matching = motif.matching,
 Next, we can visualize our network as the last step of this analysis
 
 ```{r viz_network, fig.height = 15, fig.width = 15, fig.align = "center"}
-define colors for nodes representing TFs (i.e., regulators)
+# define colors for nodes representing TFs (i.e., regulators)
 df.cor <- df.cor[order(df.cor$time_point), ]
 tfs.timepoint <- df.cor$time_point
 names(tfs.timepoint) <- df.cor$tfs
@@ -543,12 +547,12 @@ Once we generated the gene regulatory network, we can visualize individual TFs i
 Here we select two TFs for visualization.
 
 ```{r, fig.height = 3.5, fig.width = 12, fig.align = "center"}
-# obj <- AddTargetAssay(object = obj, df.grn = df.grn2)
+obj <- AddTargetAssay(object = obj, df.grn = df.grn2)
 
-# p1 <- PseudotimePlot(object = obj, tf.use = "NR3C2")
-# p2 <- PseudotimePlot(object = obj, tf.use = "RUNX1")
+p1 <- PseudotimePlot(object = obj, tf.use = "NR3C2")
+p2 <- PseudotimePlot(object = obj, tf.use = "RUNX1")
 
-# p1 + p2
+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.
@@ -560,7 +564,6 @@ As an additional step, if the spatial transcriptome (ST) data (e.g., generated b
 ```{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')
@@ -569,7 +572,7 @@ download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/VisiumSpatial
 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=FALSE}
+```{r, fig.height = 5, fig.width = 5, fig.align = "center", eval=TRUE}
 obj.spatial <- readRDS("./VisiumSpatial/AKK004_157772.rds")
 
 DefaultAssay(obj.spatial) <- "c2l_props"
@@ -587,7 +590,7 @@ 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=FALSE}
+```{r, fig.height = 5, fig.width = 10, fig.align = "center", eval=TRUE}
 DefaultAssay(obj.spatial) <- "SCT"
 
 p1 <- SpatialFeaturePlot(object = obj.spatial, features = "NR3C2") +
@@ -605,7 +608,7 @@ 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=FALSE}
+```{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,
@@ -624,7 +627,7 @@ p1 + p2
 The above plots demonstrate a clear spatial-distributed pattern for the regulon activity of NR3C2 and RUNX1.
 
 Save data
-```{r, fig.height = 5, fig.width = 6, fig.align = "center", eval=FALSE}
+```{r, fig.height = 5, fig.width = 6, fig.align = "center", eval=TRUE}
 saveRDS(obj, "./Fibroblast/coembed.trajectory.rds")
 ```
 
@@ -638,7 +641,7 @@ This will allow us to
 * 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=FALSE}
+```{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")
@@ -647,12 +650,12 @@ 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=FALSE}
+```{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=FALSE}
+```{r,fig.width=12, fig.height=12, eval=TRUE}
 p <- TopEmbGRN(df.grn=netobj,axis=c(2,3))
 ```
 
@@ -662,7 +665,7 @@ This reveals that RUNX1 is a central TF in this network, as it has high in/out-d
 
 To further inspect genes and TFs associated with RUNX1, we create a centric layout embedding.
 
-```{r,fig.width=12,fig.height=12, eval=FALSE}
+```{r,fig.width=12,fig.height=12, eval=TRUE}
 NetCentPlot(netobj,"RUNX1")
 ```