diff --git a/.nojekyll b/.nojekyll index d96fcf2c..a4376e72 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -1a667c30 \ No newline at end of file +3fdd61d7 \ No newline at end of file diff --git a/develop/04c_preprocessing_setup.html b/develop/04c_preprocessing_setup.html index 20ddcea9..e4b99a87 100644 --- a/develop/04c_preprocessing_setup.html +++ b/develop/04c_preprocessing_setup.html @@ -362,23 +362,23 @@

Underst

The --input parameter points to the samplesheet.csv file that contains all the info regarding our samples. The file looks like this:

-
- @@ -1513,7 +1513,7 @@

Saved results

- - + @@ -299,6 +299,12 @@
@@ -1255,7 +1261,7 @@

Ex
- + Summarised workflow
diff --git a/develop/09_summarized_workflow.html b/develop/09_summarized_workflow.html index c4ec6a48..24f2fe8f 100644 --- a/develop/09_summarized_workflow.html +++ b/develop/09_summarized_workflow.html @@ -64,6 +64,7 @@ + @@ -106,7 +107,7 @@ - +
@@ -143,7 +144,7 @@ Data processing
  • - + Data analyses
  • @@ -164,23 +165,164 @@ +
    • + +
    -
    +
    -
    +

    Summarised workflow

    -
    +
    @@ -1045,6 +1187,15 @@

    +
    + + Functional class scoring + +
    +
    +
    +

    diff --git a/develop/index.html b/develop/index.html index 23017464..ae7a2ee6 100644 --- a/develop/index.html +++ b/develop/index.html @@ -383,7 +383,7 @@

    Authors

    -
    +

    @@ -400,7 +400,7 @@

    -
    +

    @@ -415,7 +415,7 @@
    -
    +

    @@ -430,7 +430,7 @@
    -
    +

    diff --git a/search.json b/search.json index 308f7274..0e5d19c4 100644 --- a/search.json +++ b/search.json @@ -366,77 +366,132 @@ "href": "develop/09_summarized_workflow.html", "title": "Summarised workflow", "section": "", - "text": "We have detailed the various steps in a differential expression analysis workflow, providing theory with example code. To provide a more succinct reference for the code needed to run a DGE analysis, we have summarized the steps in an analysis below:" + "text": "We have detailed the various steps in a differential expression analysis workflow, providing theory with example code. To provide a more succinct reference for the code needed to run a DGE analysis, we have summarized the steps in an analysis below:", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#libraries", "href": "develop/09_summarized_workflow.html#libraries", "title": "Summarised workflow", "section": "Libraries", - "text": "Libraries\nlibrary(tidyverse)\nlibrary(DESeq2)\nlibrary(ggrepel)\nlibrary(pheatmap)\nlibrary(annotables)\nlibrary(clusterProfiler)\nlibrary(DOSE)\nlibrary(pathview)\nlibrary(org.Hs.eg.db)\nlibrary(tximport)\nlibrary(RColorBrewer)" + "text": "Libraries\nlibrary(tidyverse)\nlibrary(DESeq2)\nlibrary(ggrepel)\nlibrary(pheatmap)\nlibrary(annotables)\nlibrary(clusterProfiler)\nlibrary(DOSE)\nlibrary(pathview)\nlibrary(org.Hs.eg.db)\nlibrary(tximport)\nlibrary(RColorBrewer)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#obtaining-gene-level-counts-from-your-preprocessing-and-create-deseq-object", "href": "develop/09_summarized_workflow.html#obtaining-gene-level-counts-from-your-preprocessing-and-create-deseq-object", "title": "Summarised workflow", "section": "Obtaining gene-level counts from your preprocessing and create DESeq object", - "text": "Obtaining gene-level counts from your preprocessing and create DESeq object\n\nIf you have a traditional raw count matrix\nLoad data and metadata\ndata <- read_table(\"../Data/Vampirium_counts_traditional.tsv\") \n\nmeta <- read_table(\"../Data/samplesheet.csv\")\nmeta$condition <- factor(meta$condition, levels = c(\"vampirium\", \"control\", \"garlicum\"))\nCheck that the row names of the metadata equal the column names of the raw counts data\nall(colnames(data)[-1] == meta$sample)\nCreate DESeq2Dataset object\ndds <- DESeqDataSetFromMatrix(countData = data %>% column_to_rownames(\"GeneSymbol\"), \n colData = meta %>% column_to_rownames(\"sample\"), \n design = ~ condition)\n\n\nIf you have pseudocounts\nLoad samplesheet with all our metadata from our pipeline\n# Load data, metadata and tx2gene and create a txi object\nmeta <- read_csv(\"/work/Intro_bulkRNAseq/Data/samplesheet.csv\")\nmeta$condition <- factor(meta$condition, levels = c(\"vampirium\", \"control\", \"garlicum\"))\nCreate a list of salmon results\ndir <- \"/work/Intro_bulkRNAseq/Data/salmon\"\ntx2gene <- read_table(file.path(dir,\"salmon_tx2gene.tsv\"), col_names = c(\"transcript_ID\",\"gene_ID\",\"gene_symbol\"))\n\n# Get all salmon results files\nfiles <- file.path(dir, meta$sample, \"quant.sf\")\nnames(files) <- meta$sample\nCreate txi object\ntxi <- tximport(files, type=\"salmon\", tx2gene=tx2gene, countsFromAbundance = \"lengthScaledTPM\", ignoreTxVersion = TRUE)\nCreate dds object\ndds <- DESeqDataSetFromTximport(txi,\n colData = meta %>% column_to_rownames(\"sample\"), \n design = ~ condition)" + "text": "Obtaining gene-level counts from your preprocessing and create DESeq object\n\nIf you have a traditional raw count matrix\nLoad data and metadata\ndata <- read_table(\"../Data/Vampirium_counts_traditional.tsv\") \n\nmeta <- read_table(\"../Data/samplesheet.csv\")\nmeta$condition <- factor(meta$condition, levels = c(\"vampirium\", \"control\", \"garlicum\"))\nCheck that the row names of the metadata equal the column names of the raw counts data\nall(colnames(data)[-1] == meta$sample)\nCreate DESeq2Dataset object\ndds <- DESeqDataSetFromMatrix(countData = data %>% column_to_rownames(\"GeneSymbol\"), \n colData = meta %>% column_to_rownames(\"sample\"), \n design = ~ condition)\n\n\nIf you have pseudocounts\nLoad samplesheet with all our metadata from our pipeline\n# Load data, metadata and tx2gene and create a txi object\nmeta <- read_csv(\"/work/Intro_bulkRNAseq/Data/samplesheet.csv\")\nmeta$condition <- factor(meta$condition, levels = c(\"vampirium\", \"control\", \"garlicum\"))\nCreate a list of salmon results\ndir <- \"/work/Intro_bulkRNAseq/Data/salmon\"\ntx2gene <- read_table(file.path(dir,\"salmon_tx2gene.tsv\"), col_names = c(\"transcript_ID\",\"gene_ID\",\"gene_symbol\"))\n\n# Get all salmon results files\nfiles <- file.path(dir, meta$sample, \"quant.sf\")\nnames(files) <- meta$sample\nCreate txi object\ntxi <- tximport(files, type=\"salmon\", tx2gene=tx2gene, countsFromAbundance = \"lengthScaledTPM\", ignoreTxVersion = TRUE)\nCreate dds object\ndds <- DESeqDataSetFromTximport(txi,\n colData = meta %>% column_to_rownames(\"sample\"), \n design = ~ condition)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#exploratory-data-analysis", "href": "develop/09_summarized_workflow.html#exploratory-data-analysis", "title": "Summarised workflow", "section": "Exploratory data analysis", - "text": "Exploratory data analysis\nPrefiltering low count genes + PCA & hierarchical clustering - identifying outliers and sources of variation in the data:\n\nPrefiltering low count genes\nkeep <- rowSums(counts(dds)) > 0\ndds <- dds[keep,]\n\n\nRlog transformation\n# Transform counts for data visualization\nrld <- rlog(dds, \n blind=TRUE)\n\n\nPCA\nPlot PCA\nplotPCA(rld, \n intgroup=\"condition\")\n\n\nHeatmaps\nExtract the rlog matrix from the object\nrld_mat <- assay(rld)\n# Pearson correlation betweeen samples\nrld_cor <- cor(rld_mat) \n# distances are computed by rows, so we need to transponse the matrix\nrld_dist <- as.matrix(dist(t(assay(rld))))\nPlot heatmap of correlations\npheatmap(rld_cor, \n annotation = meta %>% column_to_rownames(\"sample\") %>% select(\"condition\"))\nPlot heatmap of distances with a new color range\n# Colors from the RColorBrewer package (only 6)\nheat.colors <- brewer.pal(6, \"Blues\") \n# Interpolate 100 colors\nheat.colors <- colorRampPalette(heat.colors)(100) \n\npheatmap(rld_dist, annotation = meta %>% \n column_to_rownames(\"sample\") %>% \n select(\"condition\"), color = heat.colors)" + "text": "Exploratory data analysis\nPrefiltering low count genes + PCA & hierarchical clustering - identifying outliers and sources of variation in the data:\n\nPrefiltering low count genes\nkeep <- rowSums(counts(dds)) > 0\ndds <- dds[keep,]\n\n\nRlog transformation\n# Transform counts for data visualization\nrld <- rlog(dds, \n blind=TRUE)\n\n\nPCA\nPlot PCA\nplotPCA(rld, \n intgroup=\"condition\")\n\n\nHeatmaps\nExtract the rlog matrix from the object\nrld_mat <- assay(rld)\n# Pearson correlation betweeen samples\nrld_cor <- cor(rld_mat) \n# distances are computed by rows, so we need to transponse the matrix\nrld_dist <- as.matrix(dist(t(assay(rld))))\nPlot heatmap of correlations\npheatmap(rld_cor, \n annotation = meta %>% column_to_rownames(\"sample\") %>% select(\"condition\"))\nPlot heatmap of distances with a new color range\n# Colors from the RColorBrewer package (only 6)\nheat.colors <- brewer.pal(6, \"Blues\") \n# Interpolate 100 colors\nheat.colors <- colorRampPalette(heat.colors)(100) \n\npheatmap(rld_dist, annotation = meta %>% \n column_to_rownames(\"sample\") %>% \n select(\"condition\"), color = heat.colors)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#run-deseq2", "href": "develop/09_summarized_workflow.html#run-deseq2", "title": "Summarised workflow", "section": "Run DESeq2:", - "text": "Run DESeq2:\nOptional step - Re-create DESeq2 dataset if the design formula has changed after QC analysis in include other sources of variation using\n# dds <- DESeqDataSetFromMatrix(data, colData = meta, design = ~ covariate + condition)\nRun DEseq2\n# Run DESeq2 differential expression analysis\ndds <- DESeq(dds)\nOptional step - Output normalized counts to save as a file to access outside RStudio using\nnormalized_counts <- counts(dds, normalized=TRUE)" + "text": "Run DESeq2:\nOptional step - Re-create DESeq2 dataset if the design formula has changed after QC analysis in include other sources of variation using\n# dds <- DESeqDataSetFromMatrix(data, colData = meta, design = ~ covariate + condition)\nRun DEseq2\n# Run DESeq2 differential expression analysis\ndds <- DESeq(dds)\nOptional step - Output normalized counts to save as a file to access outside RStudio using\nnormalized_counts <- counts(dds, normalized=TRUE)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#check-the-fit-of-the-dispersion-estimates", "href": "develop/09_summarized_workflow.html#check-the-fit-of-the-dispersion-estimates", "title": "Summarised workflow", "section": "Check the fit of the dispersion estimates", - "text": "Check the fit of the dispersion estimates\nPlot dispersion estimates\nplotDispEsts(dds)" + "text": "Check the fit of the dispersion estimates\nPlot dispersion estimates\nplotDispEsts(dds)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#create-contrasts-to-perform-wald-testing-or-the-shrunken-log2-fold-changes-between-specific-conditions", "href": "develop/09_summarized_workflow.html#create-contrasts-to-perform-wald-testing-or-the-shrunken-log2-fold-changes-between-specific-conditions", "title": "Summarised workflow", "section": "Create contrasts to perform Wald testing or the shrunken log2 fold changes between specific conditions", - "text": "Create contrasts to perform Wald testing or the shrunken log2 fold changes between specific conditions\nFormal LFC calculation\n# Specify contrast for comparison of interest\ncontrast <- c(\"condition\", \"control\", \"vampirium\")\n\n# Output results of Wald test for contrast\nres <- results(dds, \n contrast = contrast, \n alpha = 0.05)\nShrinkage\n# Get name of the contrast you would like to use\nresultsNames(dds)\n\n# Shrink the log2 fold changes to be more accurate\nres <- lfcShrink(dds, \n coef = \"condition_control_vs_vampirium\", \n type = \"apeglm\")" + "text": "Create contrasts to perform Wald testing or the shrunken log2 fold changes between specific conditions\nFormal LFC calculation\n# Specify contrast for comparison of interest\ncontrast <- c(\"condition\", \"control\", \"vampirium\")\n\n# Output results of Wald test for contrast\nres <- results(dds, \n contrast = contrast, \n alpha = 0.05)\nShrinkage\n# Get name of the contrast you would like to use\nresultsNames(dds)\n\n# Shrink the log2 fold changes to be more accurate\nres <- lfcShrink(dds, \n coef = \"condition_control_vs_vampirium\", \n type = \"apeglm\")", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#output-significant-results", "href": "develop/09_summarized_workflow.html#output-significant-results", "title": "Summarised workflow", "section": "Output significant results:", - "text": "Output significant results:\n# Set thresholds\npadj.cutoff <- 0.05\n\n# Turn the results object into a tibble for use with tidyverse functions\nres_tbl <- res %>%\n data.frame() %>%\n rownames_to_column(var=\"gene\") %>% \n as_tibble()\n\n# Subset the significant results\nsig_res <- filter(res_tbl, \n padj < padj.cutoff)" + "text": "Output significant results:\n# Set thresholds\npadj.cutoff <- 0.05\n\n# Turn the results object into a tibble for use with tidyverse functions\nres_tbl <- res %>%\n data.frame() %>%\n rownames_to_column(var=\"gene\") %>% \n as_tibble()\n\n# Subset the significant results\nsig_res <- filter(res_tbl, \n padj < padj.cutoff)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#visualize-results-volcano-plots-heatmaps-normalized-counts-plots-of-top-genes-etc.", "href": "develop/09_summarized_workflow.html#visualize-results-volcano-plots-heatmaps-normalized-counts-plots-of-top-genes-etc.", "title": "Summarised workflow", "section": "Visualize results: volcano plots, heatmaps, normalized counts plots of top genes, etc.", - "text": "Visualize results: volcano plots, heatmaps, normalized counts plots of top genes, etc.\nFunction to get gene_IDs based on gene names. The function will take as input a vector of gene names of interest, the tx2gene dataframe and the dds object that you analyzed.\nlookup <- function(gene_name, tx2gene, dds){\n hits <- tx2gene %>% select(gene_symbol, gene_ID) %>% distinct() %>% \n filter(gene_symbol %in% gene_name & gene_ID %in% rownames(dds))\n return(hits)\n}\n\nlookup(gene_name = \"TSPAN7\", tx2gene = tx2gene, dds = dds)\nPlot expression for single gene\nplotCounts(dds, gene=\"ENSG00000156298\", intgroup=\"condition\")\n\nMAplot\nplotMA(res)\n\n\nVolcano plot with labels (top N genes)\n## Obtain logical vector where TRUE values denote padj values < 0.05 and fold change > 1.5 in either direction\nres_tbl <- res_tbl %>% \n mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 0.58)\n## Create an empty column to indicate which genes to label\nres_tbl <- res_tbl %>% mutate(genelabels = \"\")\n\n## Sort by padj values \nres_tbl <- res_tbl %>% arrange(padj)\n\n## Populate the genelabels column with contents of the gene symbols column for the first 10 rows, i.e. the top 10 most significantly expressed genes\nres_tbl$genelabels[1:10] <- as.character(res_tbl$gene[1:10])\n\nhead(res_tbl)\nggplot(res_tbl, aes(x = log2FoldChange, y = -log10(padj))) +\n geom_point(aes(colour = threshold)) +\n geom_text_repel(aes(label = genelabels)) +\n ggtitle(\"Mov10 overexpression\") +\n xlab(\"log2 fold change\") + \n ylab(\"-log10 adjusted p-value\") +\n theme(legend.position = \"none\",\n plot.title = element_text(size = rel(1.5), hjust = 0.5),\n axis.title = element_text(size = rel(1.25))) \n\n\nHeatmap of differentially expressed genes\n# filter significant results from normalized counts\nnorm_sig <- normalized_counts %>% as_tibble(rownames = \"gene\") %>%\n dplyr::filter(gene %in% sig_res$gene) %>% column_to_rownames(var=\"gene\")\npheatmap(norm_sig, \n cluster_rows = T, #cluster by expression pattern\n scale = \"row\", # scale by gene so expression pattern is visible\n treeheight_row = 0, # dont show the row dendogram\n show_rownames = F, # remove rownames so it is more clear\n annotation = meta %>% column_to_rownames(var = \"sample\") %>% \n dplyr::select(\"condition\")\n)" + "text": "Visualize results: volcano plots, heatmaps, normalized counts plots of top genes, etc.\nFunction to get gene_IDs based on gene names. The function will take as input a vector of gene names of interest, the tx2gene dataframe and the dds object that you analyzed.\nlookup <- function(gene_name, tx2gene, dds){\n hits <- tx2gene %>% select(gene_symbol, gene_ID) %>% distinct() %>% \n filter(gene_symbol %in% gene_name & gene_ID %in% rownames(dds))\n return(hits)\n}\n\nlookup(gene_name = \"TSPAN7\", tx2gene = tx2gene, dds = dds)\nPlot expression for single gene\nplotCounts(dds, gene=\"ENSG00000156298\", intgroup=\"condition\")\n\nMAplot\nplotMA(res)\n\n\nVolcano plot with labels (top N genes)\n## Obtain logical vector where TRUE values denote padj values < 0.05 and fold change > 1.5 in either direction\nres_tbl <- res_tbl %>% \n mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 0.58)\n## Create an empty column to indicate which genes to label\nres_tbl <- res_tbl %>% mutate(genelabels = \"\")\n\n## Sort by padj values \nres_tbl <- res_tbl %>% arrange(padj)\n\n## Populate the genelabels column with contents of the gene symbols column for the first 10 rows, i.e. the top 10 most significantly expressed genes\nres_tbl$genelabels[1:10] <- as.character(res_tbl$gene[1:10])\n\nhead(res_tbl)\nggplot(res_tbl, aes(x = log2FoldChange, y = -log10(padj))) +\n geom_point(aes(colour = threshold)) +\n geom_text_repel(aes(label = genelabels)) +\n ggtitle(\"Mov10 overexpression\") +\n xlab(\"log2 fold change\") + \n ylab(\"-log10 adjusted p-value\") +\n theme(legend.position = \"none\",\n plot.title = element_text(size = rel(1.5), hjust = 0.5),\n axis.title = element_text(size = rel(1.25))) \n\n\nHeatmap of differentially expressed genes\n# filter significant results from normalized counts\nnorm_sig <- normalized_counts %>% as_tibble(rownames = \"gene\") %>%\n dplyr::filter(gene %in% sig_res$gene) %>% column_to_rownames(var=\"gene\")\npheatmap(norm_sig, \n cluster_rows = T, #cluster by expression pattern\n scale = \"row\", # scale by gene so expression pattern is visible\n treeheight_row = 0, # dont show the row dendogram\n show_rownames = F, # remove rownames so it is more clear\n annotation = meta %>% column_to_rownames(var = \"sample\") %>% \n dplyr::select(\"condition\")\n)", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#perform-analysis-to-extract-functional-significance-of-results-go-or-kegg-enrichment-gsea-etc.", "href": "develop/09_summarized_workflow.html#perform-analysis-to-extract-functional-significance-of-results-go-or-kegg-enrichment-gsea-etc.", "title": "Summarised workflow", "section": "Perform analysis to extract functional significance of results: GO or KEGG enrichment, GSEA, etc.", - "text": "Perform analysis to extract functional significance of results: GO or KEGG enrichment, GSEA, etc.\n\nAnnotate with annotables\nids <- grch38 %>% dplyr::filter(ensgene %in% res_tbl$gene) \nres_ids <- inner_join(res_tbl, ids, by=c(\"gene\"=\"ensgene\"))\n\n\nPerform enrichment analysis of GO terms (can be done as well with KEGG pathways)\n# Create background dataset for hypergeometric testing using all genes tested for significance in the results\nall_genes <- dplyr::filter(res_ids, !is.na(gene)) %>% \n pull(gene) %>% \n as.character()\n\n# Extract significant results\nsig <- dplyr::filter(res_ids, padj < 0.05 & !is.na(gene))\n\nsig_genes <- sig %>% \n pull(gene) %>% \n as.character()\n# Perform enrichment analysis\nego <- enrichGO(gene = sig_genes, \n universe = all_genes,\n keyType = \"ENSEMBL\",\n OrgDb = org.Hs.eg.db, \n ont = \"BP\", \n pAdjustMethod = \"BH\", \n qvalueCutoff = 0.05, \n readable = TRUE)\n\nego <- enrichplot::pairwise_termsim(ego)\n\n\nVisualize result\ndotplot(ego, showCategory=50)\nemapplot(ego, showCategory = 50)\n\n\nCnetplot\n## To color genes by log2 fold changes, we need to extract the log2 fold changes from our results table creating a named vector\nsig_foldchanges <- sig$log2FoldChange\n\nnames(sig_foldchanges) <- sig$gene\n## Cnetplot details the genes associated with one or more terms - by default gives the top 5 significant terms (by padj)\ncnetplot(ego, \n categorySize=\"pvalue\", \n showCategory = 5, \n foldChange=sig_foldchanges, \n vertex.label.font=6)\n\n\nPerform GSEA analysis of KEGG pathways (can be done as well with GO terms)\n# Extract entrez IDs. IDs should not be duplicated or NA\nres_entrez <- dplyr::filter(res_ids, entrez != \"NA\" & entrez != \"NULL\" & duplicated(entrez)==F)\n\n## Extract the foldchanges\nfoldchanges <- res_entrez$log2FoldChange\n\n## Name each fold change with the corresponding Entrez ID\nnames(foldchanges) <- res_entrez$entrez\n\n## Sort fold changes in decreasing order\nfoldchanges <- sort(foldchanges, decreasing = TRUE)\n# Run GSEA of KEGG\ngseaKEGG <- gseKEGG(geneList = foldchanges, # ordered named vector of fold changes (Entrez IDs are the associated names)\n organism = \"hsa\", # supported organisms listed below\n pvalueCutoff = 0.05, # padj cutoff value\n verbose = FALSE)\n\ngseaKEGG_results <- gseaKEGG@result\nhead(gseaKEGG_results)\n## Plot the GSEA plot for a single enriched pathway:\ngseaplot(gseaKEGG, geneSetID = gseaKEGG_results$ID[1])\n## Output images for a single significant KEGG pathway\npathview(gene.data = foldchanges,\n pathway.id = gseaKEGG_results$ID[1],\n species = \"hsa\",\n limit = list(gene = 2, # value gives the max/min limit for foldchanges\n cpd = 1))\nknitr::include_graphics(paste0(\"./\",gseaKEGG_results$ID[1],\".png\"))" + "text": "Perform analysis to extract functional significance of results: GO or KEGG enrichment, GSEA, etc.\n\nAnnotate with annotables\nids <- grch38 %>% dplyr::filter(ensgene %in% res_tbl$gene) \nres_ids <- inner_join(res_tbl, ids, by=c(\"gene\"=\"ensgene\"))\n\n\nPerform enrichment analysis of GO terms (can be done as well with KEGG pathways)\n# Create background dataset for hypergeometric testing using all genes tested for significance in the results\nall_genes <- dplyr::filter(res_ids, !is.na(gene)) %>% \n pull(gene) %>% \n as.character()\n\n# Extract significant results\nsig <- dplyr::filter(res_ids, padj < 0.05 & !is.na(gene))\n\nsig_genes <- sig %>% \n pull(gene) %>% \n as.character()\n# Perform enrichment analysis\nego <- enrichGO(gene = sig_genes, \n universe = all_genes,\n keyType = \"ENSEMBL\",\n OrgDb = org.Hs.eg.db, \n ont = \"BP\", \n pAdjustMethod = \"BH\", \n qvalueCutoff = 0.05, \n readable = TRUE)\n\nego <- enrichplot::pairwise_termsim(ego)\n\n\nVisualize result\ndotplot(ego, showCategory=50)\nemapplot(ego, showCategory = 50)\n\n\nCnetplot\n## To color genes by log2 fold changes, we need to extract the log2 fold changes from our results table creating a named vector\nsig_foldchanges <- sig$log2FoldChange\n\nnames(sig_foldchanges) <- sig$gene\n## Cnetplot details the genes associated with one or more terms - by default gives the top 5 significant terms (by padj)\ncnetplot(ego, \n categorySize=\"pvalue\", \n showCategory = 5, \n foldChange=sig_foldchanges, \n vertex.label.font=6)\n\n\nPerform GSEA analysis of KEGG pathways (can be done as well with GO terms)\n# Extract entrez IDs. IDs should not be duplicated or NA\nres_entrez <- dplyr::filter(res_ids, entrez != \"NA\" & entrez != \"NULL\" & duplicated(entrez)==F)\n\n## Extract the foldchanges\nfoldchanges <- res_entrez$log2FoldChange\n\n## Name each fold change with the corresponding Entrez ID\nnames(foldchanges) <- res_entrez$entrez\n\n## Sort fold changes in decreasing order\nfoldchanges <- sort(foldchanges, decreasing = TRUE)\n# Run GSEA of KEGG\ngseaKEGG <- gseKEGG(geneList = foldchanges, # ordered named vector of fold changes (Entrez IDs are the associated names)\n organism = \"hsa\", # supported organisms listed below\n pvalueCutoff = 0.05, # padj cutoff value\n verbose = FALSE)\n\ngseaKEGG_results <- gseaKEGG@result\nhead(gseaKEGG_results)\n## Plot the GSEA plot for a single enriched pathway:\ngseaplot(gseaKEGG, geneSetID = gseaKEGG_results$ID[1])\n## Output images for a single significant KEGG pathway\npathview(gene.data = foldchanges,\n pathway.id = gseaKEGG_results$ID[1],\n species = \"hsa\",\n limit = list(gene = 2, # value gives the max/min limit for foldchanges\n cpd = 1))\nknitr::include_graphics(paste0(\"./\",gseaKEGG_results$ID[1],\".png\"))", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/09_summarized_workflow.html#make-sure-to-output-the-versions-of-all-tools-used-in-the-de-analysis", "href": "develop/09_summarized_workflow.html#make-sure-to-output-the-versions-of-all-tools-used-in-the-de-analysis", "title": "Summarised workflow", "section": "Make sure to output the versions of all tools used in the DE analysis:", - "text": "Make sure to output the versions of all tools used in the DE analysis:\nsessionInfo()" + "text": "Make sure to output the versions of all tools used in the DE analysis:\nsessionInfo()", + "crumbs": [ + "Data analyses", + "Workflow summary", + "Summarised workflow" + ] }, { "objectID": "develop/rmd/05b_count_matrix.html", diff --git a/sitemap.xml b/sitemap.xml index 85721f25..654cfc7d 100644 --- a/sitemap.xml +++ b/sitemap.xml @@ -2,154 +2,154 @@ https://hds-sandbox.github.io/bulk_RNAseq_course/develop/05c_count_normalization.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/cards/AdrijaKalvisa.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/cards/extraCards/JenniferBartell.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/cards/HenrikeZschach.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/07c_DEA_visualization.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/08c_FA_GSEA.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/07b_hypothesis_testing.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/07_extra_contrast_design.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/09_summarized_workflow.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/05b_count_matrix.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/05c_count_normalization.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/08c_FA_GSEA.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/06_exploratory_analysis.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/07c_DEA_visualization.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/08a_FA_genomic_annotation.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/08a_FA_genomic_annotation.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/summary.html - 2024-11-12T15:14:07.183Z + 2024-11-12T15:30:09.451Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/04c_preprocessing_setup.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/05a_data_analysis_setup.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/06_exploratory_analysis.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/07a_DEA.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/04b_pipelines.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/02_experimental_planning.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/09_summarized_workflow.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/08b_FA_overrepresentation.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/07a_DEA.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/07b_hypothesis_testing.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/modify_md.html - 2024-11-12T15:14:07.183Z + 2024-11-12T15:30:09.451Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/rmd/07_extra_contrast_design.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/04a_preprocessing.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/index.html - 2024-11-12T15:14:07.059Z + 2024-11-12T15:30:09.327Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/03_data_explanation.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/05b_count_matrix.html - 2024-11-12T15:14:06.495Z + 2024-11-12T15:30:08.767Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/workshop_RNAseq_nov2024.html - 2024-11-12T15:14:07.183Z + 2024-11-12T15:30:09.451Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/cards/thildeTerkelsen.html - 2024-11-12T15:14:06.503Z + 2024-11-12T15:30:08.775Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/cards/AlbaMartinez.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/cards/extraCards/SamueleSoraggi.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.771Z https://hds-sandbox.github.io/bulk_RNAseq_course/develop/08b_FA_overrepresentation.html - 2024-11-12T15:14:06.499Z + 2024-11-12T15:30:08.767Z