-
Notifications
You must be signed in to change notification settings - Fork 2
Bonus (Self‐study): MacParland et al 2018
You may want to download and install the loupe browser. This is a browser designed by 10X to visualize the results of the Cell Ranger analysis that we will conduct below. Alternatively, I will load it on my machine and you can see what it looks like that way.
You'll be working with R to analyze the outputs of the scRNA-seq data analysis. If you'd like to run R on your local machine, you'll need to install R. You'll also likely want an IDE to work with. I would recommend either Rstudio or Visual Studio Code. You'll also need to install the following packages:
- dplyr
- Seurat
- patchwork
- ggplot2
- languageserver
Alternatively you can run R on the computation server (where it is already installed). But you'll likely want to connect to it using Visual Studio Code. So install this on your local machine and install the R extension and the SSH extension.
The processing of single cell RNA-seq data differs from that of bulk RNA-seq data with regards to the programs that are used to process it. However, many of the underlying principles (i.e. mapping and counting) are shared between bulk and single cell RNA-seq.
There are multiple ways to prepare single cell RNA-seq libraries. One popular way is to use the Chromium controller from 10X Genomics. This is what was used in MacParland et al 2018.
We will be working with a program called Cell Ranger which was developed by 10X Genomics.
The Cell Ranger pipeline generally starts from the raw output of the Illumina sequencing platform (raw base call format; BCL). Using the cellranger mkfastq command a set of fastq files are generated from the BCL. The set of fastq files is generated in a particular way where the various index sequences (i.e. barcode and unique molecular identifier (UMI)) are held in separate files to the RNA sequences. In Cell Ranger language, a barcode is unique to a cell (and a GEM well), and a UMI is unique to a transcript (pre-PCR).
Take a look at the data we downloaded from MacParland. You'll notice that it's in bam format, not fastq. This is normal for submission of 10X data to sequencing archives like NCBI. The bam format contains all the information of the fastq files, but with additional information on mapping.
As we touched on earlier, bam format is binary and not human readable. In order to read a bam file you will have to convert it to sam format.
Exercise: Use samtools to look at the first 10 lines of the bam files in human readable format. How will you install it? Would you use conda or install it from source?
To run the Cell Ranger analysis of our data we want to run the cellranger count method. Take a look at the diagram in the below link to see what it is doing.
What is cellranger count doing?
Here is another good resource outlining the Single-Library Analysis that we'll be conducting.
As input, cellranger count takes fastq files. But we have bam files. Fortunately, Cell Ranger also contains a program to recapitulate the fastq files from bam files. It is helpfully named bamtofastq. We can run this now.
Excercise: Install cellranger. Find the bamtofastq executable and make sure you can get it to run.
We won't run it on the actual data because that would take a lot of time and resources that we don't have available to us right now. Instead I've recreated the output files here: /home/humebc/VTK_22/macparland/results/bamtofastq
Now that we have the fastq files we can run the cellranger count program to generate the count tables.
Excercise: Run
cellranger counton one set of the fastq files. You'll need to use a reference transcriptome. These come prebuilt for certain organisms (e.g. Human) and I have already downloaded it for you. It is here:/home/humebc/VTK_22/reference/refdata-gex-GRCh38-2020-A
An explanation of the output files can be found here
Exercise: Examine the output of
cellranger countand compare it to the above documentation.
Exercise: Transfer the .html file that was generated to your local machines hard drive. Open it up.
Exercise: If you installed the Loupe browser, open up the Loupe file on your local machine. If you didn't, watch me!
Great! Now we have one sample analyzed, but what about the other samples?
This is where cellranger aggr is used. Here's the documentation on running aggr. And here's the documentation for the outputs.
Exercise: Do what's necessary to run
cellranger aggr. Is it annoying having to run each sample individually on the command line? Imagine if you had 85 samples? Now how annoying is it? Do you think there are better approaches to running many samples, not just for Cell Ranger, but in general?
The documentation of the output structure for cellranger aggr is here.
Exercise: Again, pull down the .html and cloupe files to visualize the results. Can you see the addition of the additional samples?
We've done the heavy computing on the computational server. If you remember, the sequencing files that we started with were very large and the computation would have taken a very long time on your laptop - if it was possible at all.
However, the files that we've ended up with are not so large - about 100MB.
From here, the computations are less intense and the need for a large parallelization is reduced. We'll doing the next set of computations, and generating figures, in R.
At this point we have a choice. We can either move the files that we generated and that are required for the next analysis onto our local machines and work with them there.
Alternatively, we can again 'connect' to the server and perform our analysis there.
If you choose to run on your local machine, then you'll need to install R and all of the dependencies that we'll be requiring. See the requirements sections at the beginning of this Day's section.
If you choose work on the server then you can either start up an R session on the command line (R) or you can use an IDE such as Visual Studio Code (recommended) to connect over SSH. This will make the work much easier for you.
Exercise: Get your chosen environment setup for the remainder of the analysis in R.
Seurat is a package used to analyze scRNA-seq data. It has gained great popularity in recent years and is widely used by the academic community and industry alike.
Much of the work we'll be doing to recreate the results of the MacParland analysis are well documented by the creators of the Seurat package. For example, much of what we'll be doing is covered in their pbmc3K tutorial.
Exercise: Look through the methods of the MacParland paper to see where we're at with the analysis. Critically appraise how they've written up the methods. Is it easy to follow? Is there enough detail?
The first stage of analysis in R is to create a Seurat object from the features/barcode table that we created using the cellranger aggr command.
The table that we want to import is here: /home/humebc/VTK_22/macparland/results/aggr/macparland/outs/count/filtered_feature_bc_matrix
The commands for creating a Seurat object is as follow:
# First read in the data
mc.data <- Read10X(data.dir = "/home/humebc/VTK_22/macparland/results/aggr/macparland/outs/count/filtered_feature_bc_matrix", min.cells = 3)
Note that we are screening for features that are found in at least 3 cells, the same as they did in MacParland.
Until now we haven't really been able to get a feel for what the count table looks like. Now we can:
# Preview the table
mc.data[1:5, 1:30]
# Then make the Seurat object
mc <- CreateSeuratObject(counts = mc.data, project = "macparland")
You're welcome to call the Seurat variable whatever you like, but its probably easiest if you call it mc like I have.
A common part of processing the scRNA-seq data is to filter the data for 'high quality' cells.
There are several common paramaters by which cells can be filtered. E.g. see here.
One common parameter is the percentage of UMIs that are mitochondrial in origin.
We can calculate this metric for our data as follows:
mc[["percent.mt"]] <- PercentageFeatureSet(mc, pattern = "^MT-")
Other metrics have already been calculated for the dataset as part of creating the Seurat object:
head(pbmc@meta.data, 5)
summary(mc@meta.data)
It can be helpful to visualise these metrics:
# Create a violin plot
VlnPlot(mc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# Visualise the correlation between metrics
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
Finally, we want to perform the filtering out cells from the data in the same way MacParland did. I.e.:
- filter out cell with < 1500 UMIs
- filter out cells with a high percentage of counts of mitochondrial origin
Exercise: Filter out those cells. Have a look at the pbmc3K tutorial for how to do this.
Now that we have completed the pre-processing of the data.
Exercise: Visualize the data again to make sure that the filtering has been applied. Has it?
Exercise: How would we produce the figure from the paper that shows library size plotted against mitochrondrial transcript percent?
You'll hear the term normalization a lot in computational biology and data science.
It can have different meanings depending on the context.
Here, we're normalizing the feature expression measurements for each of the individual cells according to the total expression of the cell. E.g. think about a case where one cell has twice the number of UMIs sequenced than another another. If for a given feature the first cell has a count of 4 and in the second the cell has a count of 8, is that feature upregulated in the second cell?
The authors performed normalization using an R package called scran. However, for simplicity's sake, we will continue using Seurat.
Normalize using Seurat:
mc <- NormalizeData(mc, normalization.method = "LogNormalize", scale.factor = 10000)
From here we are interested in clustering the data. That is assigning each of the cells to a given group.
To do this clustering we need to get an idea for how related each of the cells are to all of the other cells.
Exercise: Discuss, what are we basing this similarity on?
Exercise: Inspect the data, how many cells do we have? How many features?
Generally researchers choose to work on a subset of features for performing similarity analyses and clustering. This is to minimise computational load. All features could be used but do you think that all features provide the same ammount of information with regards to how similar cells are?
The most informative features will be those that vary the most between cells. In other words, they will be those features with the highest count variance across the cells.
We will identify these highly vairable features and use them for our down stream analyses:
mc <- FindVariableFeatures(mc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(mc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(mc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
The next step for clustering is to perform dimensionality reduction.
What is dimensionality reduction?
We will do this by Principal Component Analysis (PCA)
Before performing the PCA we need to standardize. This step gives equal weight to each of the features in downstream analyses, so that highly-expressed genes do not dominate.
Standardization shifts the expression of each gene, so that the mean expression across cells is 0 and scales the expression of each gene, so that the variance across cells is 1.
In Seurat this standardization is referred to as scaling.
Let's scale the data:
# Also referred to as standardization
all.genes <- rownames(mc)
mc <- ScaleData(mc, features = all.genes)
Then we can perform PCA
mc <- RunPCA(mc, features = VariableFeatures(object = mc))
And visualize the resultant Principal Components (PCs)
VizDimLoadings(mc, dims = 1:2, reduction = "pca")
ggsave("visdimloadings.png")
DimPlot(mc, reduction = "pca")
ggsave("dimplot.png")
# NB fast=F must be set in order to return a ggplot object.
DimHeatmap(mc, dims = 1, cells = 500, balanced = T, fast=F, nfeatures=100)
ggsave("dimheatmap.png")
Moving forward, we want to use a certain number of the resultant PCs. There wouldn't be much point in using all of the PCs. If we did that we wouldn't have reduced the dimensionaly of the dataset and the computational complexity would still be very high.
How many PCs do we select though?
One of the simplest ways to assess this is with an elbow plot:
ElbowPlot(mc, ndims=35)
Exercise: How many PCs should we move forwards with? How many did the authors choose?
Seurat implements a graph-based approach to clustering.
Cells are embeded in a graph structure - for example a K-nearest neighbor (KNN) graph, with edges drawn between cells with similar feature expression patterns, and then attempt to partition this graph into highly interconnected ‘quasi-cliques’ or ‘communities’.
The first part of this process is to find the Nearest Neibours of the cells:
mc <- FindNeighbors(mc, dims = 1:10)
Clusters are then called according to the Seurate algorithm from the graph:
mc <- FindClusters(mc, resolution = 0.5, pc.use=1:10)
Take a look at the cluster assignments:
head(Idents(mc), 5)
Finally we can run non-linear dimensional reduction in the form of UMAP or tSNE:
DimPlot(mc, reduction = "umap", pt.size=1)
ggsave("umap.png")
mc <- RunTSNE(mc, dims = 1:10)
DimPlot(mc, reduction = "tsne", pt.size=1)
ggsave("tsne.png")
In MacParland et al they identified the cells in the clusters by looking at cluster biomarkers. That is, genes that are unique markers of the cluster compared to all other clusters. They then used those identified genes with a manually curated set of genes that are kown to be indicative of certain cell types to assign cell types to the clusters.
We will not identify the cell type here, but we will identify the biomarkers (features) of the clusters.
Seurat includes power functionality to this end both to identify the markers and visually display the results for each of the clusters. We will produce an overview figure here as the last part of the practical:
mc.markers <- FindAllMarkers(mc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
mc.markers %>%
group_by(cluster) %>%
slice_max(n = 2, order_by = avg_log2FC)
mc.markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(mc, features = top10$gene) + NoLegend()
ggsave("cluster_diff_expression.png")