Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.
Table of Contents

...

  1. Open RStudio (you can type it in the Windows search bar)

  2. Create a new R script: ‘File’ → “New File” → “R script”

  3. Save this script where your samples folders are (‘File’ → ‘Save’). These should be on your H or W drive. Save the script file as scrnaseq.R

In the following sections you will be copying and running the R code into your scrnaseq.R script.

Cell Ranger (and nfcore/scrnaseq) generates a default folder and file output structure. There will be a main folder that contains all the sample subfolders (NOTE: this is where you must save your R script). Each sample folder will have an ‘outs’ subfolder. This ‘outs’ folder contains a ‘filtered_feature_bc_matrix’ folder, which contains the files that Seurat uses in its analysis.

...

You can manually set your working directory in RStudio by selecting ‘Session' -> 'Set working directory' -> 'Choose directory'. Choose the same directory as you saved your scrnaseq.R script, previous section. This will output the setwd(...) command with your working directory into the console window (bottom left panel). Copy this command to replace the default setwd(...) line in your R script.

...

Code Block
#### 2k. Visualise gene expression by marker ####

# PCA

p <- FeaturePlot(mat3_filt, features = markers, reduction = "pca", cols = c("red", "lightgrey"), pt.size = 1)
p

# Export as a 300dpi tiff
tiff_exp <- paste0(sample, "_pca_markers.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
# Export as a pdf
pdf_exp <- paste0(sample, "_pca_markers.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")


# UMAP

p <- FeaturePlot(mat3, features = markers, reduction = "umap", cols = c("red", "lightgrey"), pt.size = 1)
p

tiff_exp <- paste0(sample, "_umap_markers.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0(sample, "_umap_markers.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")


# t-SNE

p <- FeaturePlot(mat3, features = markers, reduction = "tsne", cols = c("red", "lightgrey"), pt.size = 1)
p

tiff_exp <- paste0(sample, "_tsne_markers.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0(sample, "_tsne_markers.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")


# Heatmap

p <- DoHeatmap(mat3, features = markers, group.bar = FALSE) + 
scale_fill_gradientn(colors = c("darkorange", "floralwhite", "dodgerblue4")) + 
theme(text = element_text(size = 16))
p

tiff_exp <- paste0(sample, "_hmap_markers.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0(sample, "_hmap_markers.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")


# Clustered heatmap

# You can also cluster the cells and then plot a heatmap based on this cluster data, to see if expression of your chosen markers is related to clusters

# First, generate the clusters for your sample
# Choose the `resolution =` score carefully (between 0 and 1). A lower score means fewer clusters. 
# You can adjust this higher or lower to see how it affects your clusters.
mat3 <- FindNeighbors(mat3, dims = 1:10)
mat3 <- FindClusters(mat3, resolution = 0.5)

# Then generate the heatmap
p <- DoHeatmap(mat3, features = markers, raster = T) + 
scale_fill_gradientn(colors = c("darkorange", "floralwhite", "dodgerblue4")) + 
theme(text = element_text(size = 16)) + labs(color = "Dose (mg)")
p

tiff_exp <- paste0(sample, "_hmap_markers_clustered.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")

pdf_exp <- paste0(sample, "_hmap_markers_clustered.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

2l. Remove non-target cells

In this section we will remove any 'non-target' cells from our dataset. 'Non-target' cells are defined as those with 0 reads (i.e. 0 expression) for our marker(s) of choice.

If you don't want to remove any cells based on expression of specific markers, skip this section

Code Block
#### 2l. Remove non-target cells ####

# View your chosen markers
zm <- mat[row.names(mat) %in% markers, ]
as.matrix(zm)

# Then we can again see a count of cells that have 0 reads (thus 0 expression) for each marker. 
# The 'All markers' row indicates the number of cells that have zero expression in any of the provided markers. 
# These are the cells that will be filtered out from your dataset
a <- length(colnames(zm))
for (i in 1:length(markers)) {a <- c(a, sum(zm[i,] == 0))}
a <- c(a, length(colnames(zm)) - sum(apply(as.matrix(zm) == 0, 2, sum) == 0))
names(a) <- c("Total cells", markers, "All markers")
as.data.frame(a)

# Enter the markers that you want to use in the filtration.
marker_rem <- c("P2ry12", "Tmem119")

# Remove cells from main Seurat object that have zero expression for **any** of these markers.
zm <- mat[row.names(mat) %in% marker_rem, ]
# This line does a sum of every column and then outputs column where this = 0 (if any cell contains reads, this will at least = 1).
zm_1 <- as.matrix(zm)[, apply(as.matrix(zm) == 0, 2, sum) == 0]
# Then we can filter the Seurat object to contain just these cell (i.e. barcode) IDs
mat3_filt <- subset(mat3_filt, cells = colnames(zm_1))

mat3_filt

2m. Clustering by gene expression