Which file and data structure is used to create volcano plot from scRNA-Seq experiment?

Question

[Chris]

Good morning all,

In RNA-Seq, we use a raw count matrix to create a volcano plot, so which data and how we can generate a volcano plot from a scRNA-Seq experiment? I use Seurat but don't know how to extract the data used as input for EnhancedVolcano(). I don't see the row as gene ID and the column as the cell 😅. Thank you so much!

[ATpoint]

A Volcano plot visualizes logFC vs -log10p(pvalue). You just need results from a differential analysis. Can you elaborate more which analysis you're doing and which step you're at.

[Chris]

I follow the instruction on several tutorials, to the end of the analysis and use FeaturePlot() to find genes differently express between 2 conditions and want to make VolcanoPlot().

[Basti]

It would be useful if you provide the code you executed to perform the differential analysis, otherwise it is difficult to know the format of the objects you created

[Chris]

I am sorry that it is pretty long but I think it is easier to help answer the question.

install.packages("harmony")
library(harmony)
library(Seurat)
library(ggplot2)
library(tidyverse)
library(gridExtra)


mtx_obj_1 <- ReadMtx(mtx = 'condition_1/matrix.mtx.gz',
                     features = 'condition_1/features.tsv.gz',
                     cells = 'condition_1/barcodes.tsv.gz')

seurat_mtx_con_1 <- CreateSeuratObject(counts = mtx_obj_1, project = 'scRNA_Seq', min.cells = 5)

seurat_mtx_con_1

seurat_mtx_1[['percent.mt']] <- PercentageFeatureSet(seurat_mtx_1, pattern = 'ˆMT-')


seurat_mtx_con_1 <- subset(seurat_mtx_con_1, subset = nCount_RNA < 40000 &
                          nFeature_RNA > 500 &
                          percent.mt <5)

seurat_mtx_con_1 <- NormalizeData(object = seurat_mtx_con_1)

VlnPlot(seurat_mtx_con_1, features = c('nFeature_RNA','nCount_RNA','percent.mt'), ncol = 3) +
  geom_smooth(method = 'lm')

FeatureScatter(seurat_mtx_con_1, feature1 = 'nCount_RNA',feature2 = 'nFeature_RNA') +
  geom_smooth(method = 'lm')

seurat_mtx_con_1 <- FindVariableFeatures(object = seurat_mtx_con_1)

top10 <- head(VariableFeatures(seurat_mtx_NC),10)

plot1 <- VariableFeaturePlot(seurat_mtx_con_1)
LabelPoints(plot = plot1, points = top10, rebel = TRUE)

mtx_obj_con_2 <- ReadMtx(mtx = 'con_2/matrix.mtx.gz',
                      features = 'con_2/features.tsv.gz',
                      cells = 'con_2/barcodes.tsv.gz')

seurat_mtx_con_2 <- CreateSeuratObject(counts = mtx_obj_con_2, project = 'CV', min.cells = 5)

seurat_mtx_con_2[['percent.mt']] <- PercentageFeatureSet(seurat_mtx_con_2, pattern = 'ˆMT-')

seurat_mtx_con_2 <- subset(seurat_mtx_con_2, subset = nCount_RNA < 40000 &
                          nFeature_RNA > 500 &
                          percent.mt <5)

seurat_mtx_con_2 <- NormalizeData(object = seurat_mtx_con_2)

VlnPlot(seurat_mtx_con_2, features = c('nFeature_RNA','nCount_RNA','percent.mt'), ncol = 3) +
  geom_smooth(method = 'lm')

FeatureScatter(seurat_mtx_con_2, feature1 = 'nCount_RNA',feature2 = 'nFeature_RNA') +
  geom_smooth(method = 'lm')

seurat_mtx_con_2 <- FindVariableFeatures(object = seurat_mtx_con_2)

merged_seurat <- merge(seurat_mtx_con_1, y = c(seurat_mtx_con_2),
                       add.cell.ids = c('con_1','con_2'),
                       project = 'scRNA_Seq')

merged_seurat$sample <- rownames([e-mail])

[e-mail] <- separate([e-mail], col = 'sample', into = c('condition','Barcode'),
                                    sep = '_')

merged_seurat <- FindVariableFeatures(merged_seurat)

merged_seurat <- ScaleData(merged_seurat)

merged_seurat <- RunPCA(merged_seurat)

merged_seurat <- RunUMAP(merged_seurat, reduction = "pca", dims = 1:20)

merged_seurat <- FindNeighbors(merged_seurat, reduction = "pca", dims = 1:20)

merged_seurat <- FindClusters(merged_seurat, resolution = 0.5)

ElbowPlot(merged_seurat)

[e-mail]$condition_1 <- [e-mail]$condition

before <- DimPlot(merged_seurat, reduction = 'umap', group.by = 'condition')
before

CV.harmony <- merged_seurat %>%
  RunHarmony(group.by.vars = 'condition', plot_convergence = F)

CV.harmony$nCount_RNA

[e-mail] <- unite([e-mail], "condition_cluster", condition_1, seurat_clusters, sep = "_")

CV.harmony@reductions

CV.harmony.embed <- Embeddings(CV.harmony,'harmony')
CV.harmony.embed[1:10,1:10]

CV.harmony <- CV.harmony %>%
  RunUMAP(reduction = 'harmony', dims = 1:20) %>%
  FindNeighbors(reduction = 'harmony', dims = 1:20) %>%
  FindClusters(resolution = 0.5)

after <- DimPlot(CV.harmony, reduction = 'umap', group.by = 'condition')
before|after

cluster <- DimPlot(CV.harmony, reduction = 'umap', group.by = 'seurat_clusters', label = T)
condition <- DimPlot(CV.harmony, reduction = 'umap', group.by = 'condition')
condition|cluster

Idents(CV.harmony) <- CV.harmony$condition_cluster

CV.markers <-  FindAllMarkers(CV.harmony,
                              logfc.threshold = 0.25,
                              min.pct = 0.1,
                              only.pos = F)

EnhancedVolcano(?, x = "log2FoldChange", y = "padj", lab = ?,
                pCutoff = 1e-4, FCcutoff = 1)

[Basti]

EnhancedVolcano(CV.markers, x = "avg_log2FC", y = "p_val", lab =   rownames(CV.markers),
                pCutoff = 1e-4, FCcutoff = 1)

"p_val" for nominal p-values but "p_val_adj" could be used for adjusted p-values

[Chris]

Thank you for your help! I got an error that I could not find the solution on the Internet:

library(EnhancedVolcano)

Loading required package: ggrepel Error in value[3L] : Package ‘ggrepel’ version 0.9.2 cannot be unloaded: Error in unloadNamespace(package) : namespace ‘ggrepel’ is imported by ‘Seurat’ so cannot be unloaded

Could you have a suggestion?

Update. I figured it out.