5. Differential Gene Expression Analysis
After scaling, normalizing, and clustering the cells from the GEX
libraries we can now investigate which genes are differentially
expressed between either clusters or samples. First, we can investigate
the genes that define each of the clusters by using the GEX output of
the VGM function. Depending on the size of the dataset and the number of
cells this function can be quite slow. The output of this function is a
list in which each element contains the differentially expressed genes
for a given cluster. For example, the first element of the list will
correspond to a data frame describing the genes for cluster0 that we
previously observed on the UMAP. This list object will correspond to the
length of the number of clusters that were previously calculated, and
has the same format as the FindMarkers function from Seurat. The pct.1
will correspond to the percentage of cells expressing the gene in the
cluster of interest, and the pct.2 to the percentage of cells in all
other clusters.
VGM[[2]] <- SeuratObject::JoinLayers(VGM[[2]])
gene_expression_cluster <- GEX_cluster_genes(VGM[[2]], min.pct = 0.2)
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session
length(gene_expression_cluster) # length of 12, corresponding to 12 clusters
## [1] 12
length(unique(VGM[[2]]$seurat_clusters)) # length of 12
## [1] 12
print(sapply(gene_expression_cluster,nrow)) #Nr of differentially expressed genes per cluster
## [1] 1921 1266 1238 1424 2516 1814 1863 1617 1611 1513 1743 1409
We can now look at some of the genes associated with cluster0 in the
previously displayed umap
head(gene_expression_cluster[[1]])
## p_val avg_logFC pct.1 pct.2 p_val_adj SYMBOL cluster
## NKG7 2.631427e-306 3.732392 0.982 0.104 8.336623e-302 NKG7 0
## CCL5 1.070166e-286 3.434706 0.972 0.112 3.390392e-282 CCL5 0
## CTSW 5.876528e-254 3.307826 0.872 0.079 1.861743e-249 CTSW 0
## CD8A 2.962265e-251 3.746899 0.807 0.059 9.384753e-247 CD8A 0
## CD8B1 1.725846e-237 3.285546 0.842 0.080 5.467653e-233 CD8B1 0
## THY1 1.075640e-235 2.854907 0.941 0.145 3.407735e-231 THY1 0
Raw tables are informative, but not visually appealing.
Differentially expressed genes may be plotted in different modes using
Platypus.
It is also possible to create a heatmap displaying differentially
expressed genes for each cluster. This can be customized to sub-sample
cells in case certain clusters are too large for visualization purposes.
Additionally, the user can determine the number of genes to display for
each cluster based on the n.genes.per.cluster argument. The function
GEX_cluster_genes_heatmap can be used to produce a ggplot object, based
on the DoHeatmap from Seurat.
agedCNS_heatmap_clusters <- GEX_cluster_genes_heatmap(GEX = VGM[[2]], GEX_cluster_genes.output = gene_expression_cluster,n.genes.per.cluster = 3,max.cell = 30,metric = "avg_logFC", platypus.version = "v3")
print(agedCNS_heatmap_clusters)
After plotting the ggplot object we can clearly see genes enriched in
various clusters - mainly indicated by the diagonal.
Further, volcano plots may be used to show and annotate a larger
number of differentially expressed genes
agedCNS_heatmap_volcano <- GEX_volcano(DEGs.input = gene_expression_cluster, input.type = "cluster.genes", RP.MT.filter = T, color.p.threshold = 0.01, n.label.up = 10, n.label.down = 10)
## [1] "Top 10 and bottom 10 genes will be labeled"
print(agedCNS_heatmap_volcano[[1]]) #genes specific to cluster 0
We can additionally extract the top N genes per cluster directly
(with filtering) using the following function:
top_10_genes_per_cluster <- GEX_topN_DE_genes_per_cluster(GEX_cluster_genes.output = gene_expression_cluster, n.genes = 10, by_FC = T)
head(top_10_genes_per_cluster)
## p_val avg_logFC pct.1 pct.2 p_val_adj SYMBOL cluster
## CD81 9.129362e-53 -11.663484 0.000 0.357 2.892273e-48 CD81 0
## BLK 5.337741e-39 -10.261625 0.000 0.278 1.691050e-34 BLK 0
## FCRLA 2.255021e-38 -10.195180 0.000 0.274 7.144132e-34 FCRLA 0
## MZB1 4.144279e-36 -10.111562 0.000 0.260 1.312949e-31 MZB1 0
## H2-DMB2 1.198041e-49 -9.893185 0.002 0.342 3.795514e-45 H2-DMB2 0
## SPI1 5.434681e-34 -9.487847 0.000 0.247 1.721761e-29 SPI1 0
We can extract and test for differential expressed genes between the
two samples (or between other subgroups) by using the
GEX_DEgenes_persample function.
This function allows us to also create a heatmap displaying the top
most up- or downregulated genes for each cluster based on log fold
change (avg_logFC) or p value (adj_p_value). Additionally, the user can
determine the number of up- and downreulated genes to be displayed for
each sample. In this case, the output returns a list where the first
element contains the dataframe with the differntial expression
information and the second element contains the heatmap displaying the
most up-/downregulated genes.
If more than two samples are to be examined please refer to
GEX_pairwise_degs two chunks below
DE_genes_samples <- GEX_DEgenes(GEX = VGM[[2]],min.pct = .25, grouping.column = "sample_id",group1 = "s1", group2 = "s2",return.plot = "volcano",up.genes = 5,down.genes = 5,logFC = F)
#This function is flexible and takes any column name as grouping.column to allow easy exploration of differences between custom groups
head(DE_genes_samples[[1]])
## p_val avg_logFC pct.1 pct.2 p_val_adj SYMBOL signif_pval
## RGS10 1.564907e-13 1.649797 0.436 0.138 4.957783e-09 RGS10 TRUE
## CTSD 3.522094e-12 2.038039 0.568 0.297 1.115835e-07 CTSD TRUE
## C1QC 7.536629e-11 1.655614 0.375 0.134 2.387680e-06 C1QC TRUE
## CST3 2.673577e-10 1.862194 0.679 0.455 8.470159e-06 CST3 TRUE
## HPGDS 3.129809e-10 2.328784 0.287 0.077 9.915547e-06 HPGDS TRUE
## C1QA 5.307456e-10 1.574627 0.372 0.138 1.681455e-05 C1QA TRUE
## signif_logfc color_threshold
## RGS10 TRUE TRUE
## CTSD TRUE TRUE
## C1QC TRUE TRUE
## CST3 TRUE TRUE
## HPGDS TRUE TRUE
## C1QA TRUE TRUE
The same function can also return DEGs between seurat clusters or any
2 groups of cells Here DEGs between cluster 1 and 3 are calculated
DE_genes_cl1_vs_3 <- GEX_DEgenes(GEX= VGM[[2]],min.pct = .25, grouping.column = "seurat_clusters",group1 = "1", group2 = "3",return.plot = "heatmap",up.genes = 10,down.genes = 10,logFC = F)
## [1] "MS4A4B" "CD3G" "CD3D" "CD3E" "THY1" "SKAP1" "MS4A6B"
## [8] "LCK" "TMSB10" "GM2682" "H2-AA" "H2-AB1" "H2-EB1" "CD79A"
## [15] "CD74" "CD79B" "LY6D" "EBF1" "BANK1" "H2-DMB2"
#head(DE_genes_cl1_vs_3[[1]])
DE_genes_cl1_vs_3[[2]]
To gain a complete overview of DEGs between groups or clusters the
GEX_pairwise_degs function is used This function calculates DEGs for
every possible pairwise comparison between groups of the selected
column. It is thereby recommended to use this function with a maximum
number of 10 groups/clusters/samples This function automatically saves
plots where top label.n.top.genes genes are labeled by their p value.
Additionally, genes of special interest that should be labeled
irregardless of their p value can be supplied to genes.to.label
DE_clusters_all <- GEX_pairwise_DEGs(GEX = VGM[[2]], group.by = "seurat_clusters", min.pct = 0.25, RP.MT.filter = T, label.n.top.genes = 10, genes.to.label = c("CD74", "EBF1"), save.plot = F)
## Calculating pairwise DEGs 1 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 10
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 2 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 4
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 3 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 9
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 4 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 6
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 5 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 7
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 6 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 2
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 7 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 8 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 9 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 10 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 11 of 66
## [1] 1
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 12 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 4
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 13 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 9
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 14 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 6
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 15 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 7
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 16 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 2
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 17 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 18 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 19 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 20 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 21 of 66
## [1] 10
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 22 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 9
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 23 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 6
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 24 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 7
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 25 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 2
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 26 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 27 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 28 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 29 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 30 of 66
## [1] 4
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 31 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 6
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 32 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 7
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 33 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 2
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 34 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 35 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 36 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 37 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 38 of 66
## [1] 9
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 39 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 7
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 40 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 2
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 41 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 42 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 43 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 44 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 45 of 66
## [1] 6
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 46 of 66
## [1] 7
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 2
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 47 of 66
## [1] 7
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 48 of 66
## [1] 7
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 49 of 66
## [1] 7
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 50 of 66
## [1] 7
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 51 of 66
## [1] 7
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 52 of 66
## [1] 2
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 0
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 53 of 66
## [1] 2
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 54 of 66
## [1] 2
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 55 of 66
## [1] 2
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 56 of 66
## [1] 2
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 57 of 66
## [1] 0
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 3
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 58 of 66
## [1] 0
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 59 of 66
## [1] 0
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 60 of 66
## [1] 0
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 61 of 66
## [1] 3
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 5
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 62 of 66
## [1] 3
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 63 of 66
## [1] 3
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 64 of 66
## [1] 5
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 8
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 65 of 66
## [1] 5
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
## Calculating pairwise DEGs 66 of 66
## [1] 8
## Levels: 11 < 8 < 5 < 3 < 0 < 2 < 7 < 6 < 9 < 4 < 10 < 1
## [1] 11
## Levels: 1 < 10 < 4 < 9 < 6 < 7 < 2 < 0 < 3 < 5 < 8 < 11
Now that cluster defining genes have been examined, a helper plotting
function allows to visualize selected genes.
dottile <- GEX_dottile_plot(GEX = VGM[[2]], genes = c("CD19", "CD74","SDC1", "EBF1","PTPRC","CD93","CD38","CD24A","CD34","CD1D1","CR2","MS4A1","CXCR5","SELL","CD40","CD83","H2-AB1","H2-EB1","CD27","POU2AF1","NT5E","FAS","PDCD1LG2","PRDM1","ITGAM","IL10","IL12A","HAVCR2"), group.by = "seurat_clusters", threshold.to.plot = 5)
#threshold.to.plot specifies how many % of cells have to express a gene to show a dot.
dottile + ggplot2::theme(plot.title = ggplot2::element_blank(), plot.subtitle = ggplot2::element_blank(), legend.position = "bottom")
## Warning: Removed 192 rows containing missing values (`geom_point()`).
#For visualisation purposes in the RMD format
To examine patterns of co-expression two functions allow to generate
both overview and single-cell resolution plots.
#overview
coexpression_dotmap <- GEX_coexpression_coefficient(GEX = VGM[[2]], genes = c("CD19", "CD74","SDC1", "EBF1","PTPRC","CD93","CD38","CD24A","CD34"), plot.dotmap = T)
## Calculating coexpression for 9 genes with 2064 cells
coexpression_dotmap + ggplot2::theme(legend.position = "bottom")
#detail of two selected genes
GEX_scatter_coexpression(GEX = VGM[[2]], gene.1 = "CD19", gene.2 = "EBF1", color.theme = "darkred")
#to save use:
#ggsave(last_plot(), filename = "Coexpression_scatter.png")
5. VDJ Repertoire anaylsis
5.1 Checking VGM output
Now we can analyze the VDJ repertoire data before integrating VDJ and
GEX information. This may be useful if only VDJ libraries have been
sequenced without the accompanying gene expression data. We first start
by examining the structure of VGM[[1]]
print(colnames(VGM[[1]]))
## [1] "barcode" "sample_id"
## [3] "group_id" "clonotype_id_10x"
## [5] "clonotype_id" "clonotype_frequency"
## [7] "celltype" "Nr_of_VDJ_chains"
## [9] "Nr_of_VJ_chains" "VDJ_cdr3s_aa"
## [11] "VJ_cdr3s_aa" "VDJ_cdr3s_nt"
## [13] "VJ_cdr3s_nt" "VDJ_umis"
## [15] "VJ_umis" "VDJ_chain_contig"
## [17] "VJ_chain_contig" "VDJ_chain"
## [19] "VJ_chain" "VDJ_vgene"
## [21] "VJ_vgene" "VDJ_dgene"
## [23] "VDJ_jgene" "VJ_jgene"
## [25] "VDJ_cgene" "VJ_cgene"
## [27] "VDJ_sequence_nt_raw" "VJ_sequence_nt_raw"
## [29] "VDJ_sequence_nt_trimmed" "VJ_sequence_nt_trimmed"
## [31] "VDJ_sequence_aa" "VJ_sequence_aa"
## [33] "VDJ_raw_ref" "VJ_raw_ref"
## [35] "VDJ_trimmed_ref" "VJ_trimmed_ref"
## [37] "VDJ_raw_consensus_id" "VJ_raw_consensus_id"
## [39] "orig_barcode" "specifity"
## [41] "affinity" "GEX_available"
## [43] "orig.ident" "seurat_clusters"
## [45] "PC_1" "PC_2"
## [47] "tSNE_1" "tSNE_2"
## [49] "batch_id"
If GEX data has been integrated, this dataframe also contains
dimensional reduction coordinates and cluster id. If trim.and.align was
set to TRUE, the columns VDJ_sequence_nt_trimmed to VJ_trimmed_ref
contain aligned and reference sequences. The clonotype_frequency column
takes input from the clonotyping output of cellranger, which is saved
int the clonotype_id_10x column.
This object contains both cells with less than 2 V(D)J chains as well
as with more than 2 To show how this is done the VDJ_cgene column is
taken as an example
print(unique(VGM[[1]]$VDJ_cgene))
## [1] "IGHM" "IGHD" "" "IGHM;IGHM" "IGHG2B"
## [6] "IGHM;IGHD" "IGHG1" "IGHA" "IGHM;IGHG2C"
Apart from classical isotypes, we can see an empty (““) element. This
is for cells, for which only a VJ chain is available. Secondly we can
see 3 examples of cells containing more than one VDJ chain. These are
delimited by”;“. Importantly, this format is maintained throughout all
columns.
Both cells with less and more than 2 chains can be cumbersome for
analysis. Platypus V3 functions generally support such cells. To check
their quantity and, if needed, exclude these, the following line can be
used
For advanced filtering and selection options concerning cells with
aberrant number of chains, please refer to the detailed VDJ_GEX_matrix
vignette.
cat(" Cell count by number of VDJ chains")
## Cell count by number of VDJ chains
print(table(VGM[[1]]$Nr_of_VDJ_chains))
##
## 0 1 2
## 92 831 22
cat("\n Cell count by number of VJ chains")
##
## Cell count by number of VJ chains
print(table(VGM[[1]]$Nr_of_VJ_chains))
##
## 0 1 2
## 62 820 63
#Subset the VGM matrix to only include cells with 1 VDJ and 1VJ chain
#VGM[[1]] <- subset(VGM[[1]], Nr_of_VJ_chains == 1 & Nr_of_VJ_chains == 1)
5.2 Changing the clonotype strategy
Often paired nucleotide CDRH3 and CDRL3 clonotyping may not be the
best strategy given somatic hypermutation may occur in the CDR3 region.
Therefore, there could be highly similar clones that likely bind the
same antigen that are officially part of different clonal families. To
address this, we have added a function that allows for various heuristic
clonotyping strategies. This involves clonotyping by identical amino
acid CDRH3 + CDRL3 seuqence, identical germline usage, or seqeunce
homology requirements. This function works both for V3 and V2 of
platypus. The version needs to be specified. For V3 (VGM)
platypus.version a new clonotyping function was added to improve the
existing one. For the full documentation see VDJ_clonotype_v3(). The
deprecated function VDJ_clonotype may still be used for Platypus v2
input Three new parameters have been added therein: 1. global.clonotype
If set to TRUE, clonotypes will be generated across all samples of the
VGM. This may be useful if more than one sample has been taken from the
same animal (e.g. spleen and bone marrow). Moreover this can provide an
alternative approach to investigate clone sharing across donors or mice.
After global clonotyping, the resulting clonotype_id can be used as a
grouping variable and summary statistics can be generated per clone
(including its presence in multiple repertoires.) A key update as part
of the hierarchical clonotyping strategy is the same-sample-criterion.
This is explained below. 2. hierarchical For cells with aberrant numbers
of VDJ and VJ chains: the functions implements a hierarchical
clonotyping strategy for such cells. In brief, clonotypes are determined
based on all cells with exactly 1VDJ and 1VJ chain. Other cells are then
“joined in” in post. This follows three steps 1. hierarchical = “none”
Create clonotypes from 1 VDJ 1 VJ cells. 2. hierarchical =
“single.chains” Cells with a single chain (e.g. 0VDJ and 1VJ) are merged
into existing clonotypes containing the same or a homologous VJ chain,
depending on the clonotyping strategy. Here the same-sample-criterion is
applied when using global clonotyping. If multiple matching clonotypes
have been found for that cell then ones within the same sample_id are
prioritized. Further if multiple matches exist, the cell is merged into
the existing clone with the highest frequency 2. hierarchical =
“double.and.single.chains” A combination of 1VDJ and 2VJ chain will
first be checked for frequency in the dataset (by exact nucleotide CDR3
matching). Only if this count exceeds the triple.chain.count.threshold,
the clone is used as a “hub clone”. This protects from merging
clonotypes on the basis of rare doublets. Existing clonotypes are then
merged into the 1 VDJ 2 VJ clonotypes as they match with the assumption
that e.g. a cell with 1 VDJ 1 VJ is part of that same clonotype, but
missing a VJ chain due to stochastical sampling
In the example that follows, we use VDJ_clonotype to group cells into
clones based on identical CDRH3 + CDRL3 amino acid sequence. We will
compare this to the case in which we group B cells by using the same
germline genes (both heavy chain and light chain).
VGM[[1]] <- VDJ_clonotype(VDJ = VGM[[1]], clone.strategy = "cdr3.aa", global.clonotype = F, VDJ.VJ.1chain = F, hierarchical = "single.chains")
## Clonotyping strategy: cdr3.aa
## Filtered out 22 cells containing more than one VDJ AND VJ chain or two VDJ chains, as these likely correspond to doublets
## Processing sample 1 of 4
## Attempting to merge in 9 aberrant cells
## Processing sample 2 of 4
## Attempting to merge in 28 aberrant cells
## Processing sample 3 of 4
## Attempting to merge in 23 aberrant cells
## Processing sample 4 of 4
## Attempting to merge in 94 aberrant cells
## Backing up 10x default clonotyping in columns clonotype_id_10x and clonotype_frequency_10x before updating clonotype_id and clonotype_frequency columns
cat(" Nr and distribution of clonotypes using exact CDR3.aa matching \n")
## Nr and distribution of clonotypes using exact CDR3.aa matching
print(length(unique(VGM[[1]]$clonotype_id_cdr3.aa)))
## [1] 321
print(table(VGM[[1]]$clonotype_frequency_cdr3.aa)) #Check distribution of clonotypes with identical CDR3 aa sequences
##
## 1 2 3 4 5 6 7 8 9 12 24
## 593 118 69 24 15 36 7 16 9 12 24
VGM[[1]] <- VDJ_clonotype(VDJ = VGM[[1]], clone.strategy = "VDJJ.VJJ", global.clonotype = F, VDJ.VJ.1chain = F, hierarchical = "double.and.single.chains", triple.chain.count.threshold = 4)
## Clonotyping strategy: VDJJ.VJJ
## Filtered out 0 cells containing more than one VDJ AND VJ chain or two VDJ chains, as these likely correspond to doublets
## Found 0 exact matching clones with 3 chains and a frequency of at least 4. These will be used as high confidence clonotypes.
## Processing sample 1 of 4
## Attempting to merge in 9 aberrant cells
## Processing sample 2 of 4
## Attempting to merge in 28 aberrant cells
## Processing sample 3 of 4
## Attempting to merge in 23 aberrant cells
## Processing sample 4 of 4
## Attempting to merge in 94 aberrant cells
cat("\n Nr and distribution of clonotypes using germline gene matching \n ")
##
## Nr and distribution of clonotypes using germline gene matching
##
print(length(unique(VGM[[1]]$clonotype_id_VDJJ.VJJ)))
## [1] 295
print(table(VGM[[1]]$clonotype_frequency_VDJJ.VJJ)) #Check distribution of clonotypes with identical germline genes
##
## 1 2 3 4 5 6 7 8 10 12 24
## 534 148 84 32 30 12 21 16 10 12 24
Columns with new clonotype information are added to the VGM so to
allow for multiple different strategies to be compared easily. For most
Platypus downstream functions, the columns default columns
“clonotype_id” and “clonotype_frequency” are used. To control which
clonotyping strategy is used in these functions, the relevant columns
can be copied over
VGM[[1]]$clonotype_id <- VGM[[1]]$clonotype_id_cdr3.aa
VGM[[1]]$clonotype_frequency <- VGM[[1]]$clonotype_frequency_cdr3.aa
5.4 Visualizing clonal frequencies (V3 only)
For a basic view of clonal expansion the VDJ_clonal_donut produces
circular plots per sample. Label position and size should only be
adjusted after deciding on a plot export format. This function uses the
clonotype_id column as input. If VDJ_clonotyping function was used, its
result are stored in that column. If not, the default 10x clonotyping is
stored there. To retrieve default 10x clonotyping, use the column
clonotyping_id_10x
donuts <- VDJ_clonal_donut(VDJ = VGM[[1]], expanded.colors = c("grey50", "grey65", "grey80"), non.expanded.color = "black", counts.to.use = "freq_column")
## Warning in VDJ_clonal_donut(VDJ = VGM[[1]], expanded.colors = c("grey50", :
## Column name for counts.to.use was not found in VDJ. Defaulting to
## 'clonotype_id_10x'
## Using column clonotype_id for counting clones
## Clones: Expanded: 5 / 6.85%; 1 cell 68 / 93.15%; total: 73
## Cells: Expanded: 10 / 12.82%; 1 cell 68 / 87.18%; Total: 78
## Clones: Expanded: 14 / 10.85%; 1 cell 115 / 89.15%; total: 129
## Cells: Expanded: 45 / 28.12%; 1 cell 115 / 71.88%; Total: 160
## Clones: Expanded: 7 / 3.5%; 1 cell 193 / 96.5%; total: 200
## Cells: Expanded: 15 / 7.21%; 1 cell 193 / 92.79%; Total: 208
## Clones: Expanded: 64 / 18.34%; 1 cell 285 / 81.66%; total: 349
## Cells: Expanded: 214 / 42.89%; 1 cell 285 / 57.11%; Total: 499
#Counts to use = "freq_column" uses the counts in the clonotype_frequency column. Counts to use = "VGM" simply counts the rows of a given clonotype in the VGM table (these counts may differ if cells have been filtered out due to overlapping barcodes or if another clonotyping strategy was used)
donuts[[1]]
The black section shows cells in unexpanded clones. The rest shows
expanded clones by their frequency. The middle label shows the total
number of clones (cells)
5.5 Calculating common repertoire diversity metrics
This requires the Vegan package. We can calculate the diversity for
any column(s) in the VGM. If more than one column is provided, the
content will be pasted together before calculating diversity metrics
#Shannon Evenness for the VDJ chain CDR3
diversity_plot <- VDJ_diversity(VDJ = VGM[[1]],feature.columns = c("VDJ_cdr3s_aa"),grouping.column = "sample_id",metric = c("shannonevenness"), subsample.to.same.n = T)
diversity_plot
#Gini-Simpson index for pasted VDJ and VJ chain CDR3s
diversity_plot <- VDJ_diversity(VDJ = VGM[[1]],feature.columns = c("VDJ_cdr3s_aa", "VJ_cdr3s_aa"),grouping.column = "sample_id", metric = c("ginisimpson"), subsample.to.same.n = T)
diversity_plot
#exact values can be retrived as
#Jaccard index between repertoires of the two samples
diversity_plot <- VDJ_diversity(VDJ = VGM[[1]],feature.columns = c("VDJ_cgene"),grouping.column = "sample_id",metric = c("jaccard"), subsample.to.same.n = T)
diversity_plot
5.6 Sequence similarity networks
Other functions are specifically tailored to repertoire analysis -
such as VDJ_network, which creates a sequence similarity network between
repertoires or within a repertoire by connecting those clones with
sequence similarity. This function relies upon igraph to visually
display and construct the graph - which means that networks with high
number of sequences will not display easily. In the following example we
are using a small dataset. In case of a bigger dataset one may subsample
by using the sample.n function on VGM[[1]]. Setting the per.mouse
argument to false indicates that one network for multiple repertoires
should be produced.
This function is quite specialized and we are working on building a
better networks module. But hopefully the code can be used as a template
for your own analyses.
#subsampled_VGM <- dplyr::sample_n(VGM[[1]], 60)
cl_network <- VDJ_network(VDJ = VGM[[1]],per.sample = T,distance.cutoff = 8, platypus.version = "v3")
for(i in 1:2){
igraph::plot.igraph(cl_network[[1]][[i]],vertex.label=NA,vertex.size=7)
}
#For a plot including clonotype frequencies
cl_network <- VDJ_network(VDJ = VGM[[1]],per.sample = F,distance.cutoff = 8, platypus.version = "v3")
igraph::plot.igraph(cl_network[[4]],vertex.label=NA,vertex.size=3+(0.95*as.numeric(cl_network[[2]]$clonotype_frequency)),vertex.color= as.factor(VGM[[1]]$sample_id))
For more details see the documentation of the VDJ_network function,
but essentially information such as clonal frequency and which sample
(here still indicated by the “sample_id” column) are stored in the
second element of the output list. Here we can see only a few clones
that are showing connections (produced by edges between those with a
distance of 8 amino acids or less between heavy and light chain paired
CDR3 sequence homology)
5.7 Germline gene usage heatmaps
It is also possible to produce heatmaps of the germline gene usage in
the context of heavy chain V gene and light chain V gene. The output of
the VDJ_Vgene_usage function is a matrix for each repertoire
corresponding to the order specified by VDJ_analyze. The outer list
corresponds to the sample and the inner list corresponds to a matrix,
where the rows correspond to the heavy chain V genes and the columns
correspond to the light chains of the V genes. Therefore the
output[[1]][i,j] corresponds to the number of clones using the
combination of IGH-Vgene[i] and IGK/L-Vgene[j].
#First calculate adjacency matrix for V gene usage
Vgene_usage <- VDJ_Vgene_usage(VDJ = VGM[[1]], platypus.version = "v3")
#library(pheatmap)
pheatmap::pheatmap(Vgene_usage[[1]],show_rownames = F,show_colnames = F)
print(class(Vgene_usage[[1]]))
## [1] "matrix" "array"
print(head(rownames(Vgene_usage[[1]])))
## [1] "IGHV8-8" "IGHV1-9" "IGHV1-26" "IGHV14-1" "IGHV5-17" "IGHV5-16"
print(head(colnames(Vgene_usage[[1]])))
## [1] "IGKV10-96" "IGKV6-32" "IGKV19-93" "IGKV8-24" "IGKV1-110" "IGKV4-55"
This can then be easily plotted as a heatmap to observe patterns between
repertoires or can be used to calculate V gene correlation using the
“pheatmap” package.
Platypus also allows a separate analysis of V gene usage for HC and
LC. The VDJ_Vgene_usage_barplot allows the user to plot most frequently
used IgH or IgK/L V genes. By default, this function only provides
visualizations for the HC V genes, but can also provide for the LC if
LC. Vgene is set to TRUE. The User can also select the number of most
used genes to be depicted.
Vgene_usage_barplot <- VDJ_Vgene_usage_barplot(VDJ = VGM[[1]], HC.gene.number = 10, LC.Vgene = T, LC.gene.number = 10, platypus.version = "v3")
## Sample order: s1 ; s2 ; s3 ; s4
Vgene_usage_barplot[[1]]
#VDJ chains
Vgene_usage_stackedbarplot <- VDJ_Vgene_usage_stacked_barplot(VDJ = VGM[[1]], LC.Vgene = F,HC.gene.number = 10, Fraction.HC = 1, platypus.version = "v3")
## Sample order: s1 ; s2 ; s3 ; s4
Vgene_usage_stackedbarplot
Furthermore, we can also produce a circular visualization of how V
and J genes are combined throughout the repertoire. In the example that
follows we use VDJ_VJ_usage_circos to look at the V gene with the
corresponding J gene for each expanded clonotype.
#VDJ and VJ V and J gene pairing
vj_circos_bcells <- VDJ_VJ_usage_circos(VGM = VGM[[1]], c.threshold = 5,label.threshold=5,cell.level = T, platypus.version = "v3")
## Chosen clonotype column: clonotype_id_10x
## Chosen clonotype.frequency column: clonotype_frequency_10x
## Chosen clonotype.frequency column: clonotype_frequency
## Processing sample 1
## Processing sample 2
#VDJ and VJ pairing
HL_circos_bcells <- VDJ_alpha_beta_Vgene_circos(VGM = VGM[[1]], c.threshold = 5,label.threshold=5,cell.level = T, V.or.J= "both", platypus.version = "v3")
## Chosen clonotype.frequency column: clonotype_frequency
## [1] "s2"
## [1] "s4"
## Plotting for sample s2
## Plotting for sample s4
5.8 Assessing CDR3 sequence similarity
Finally, one can also look at any specific HC and LC CDR3 amino acid
patterns arising across the different clones. Using the
VDJ_logoplot_vector function the user can plot a logoplot of the CDR3
region od a certain length, specyfied by the length_cdr3 argument. For
instance, the logoplot below corresponds to all those CDR3 aminoacid
sequences of length 25
pasted_CDR3s <- paste0(VGM[[1]]$VDJ_cdr3s_aa, VGM[[1]]$VJ_cdr3s_aa)
agedCNS_CDR3_logoplot <- VDJ_logoplot_vector(cdr3.vector = pasted_CDR3s, seq_type = "aa", length_cdr3 = "auto")
## Returning logoplot based on 106 sequences of a total 945 input sequences
print(agedCNS_CDR3_logoplot)
5.9 Assessing Nr of variants per clone
Depending on the clonotype strategy, within a clone several different
variants of CDR3s or full length sequences may be contained. To count
variants, the function VDJ_variants_per_clone is used. It returns a
table per sample_id with stats about variants
variants_agedCNS <- VDJ_variants_per_clone(VDJ = VGM[[1]], variants.of = c("VDJ_cdr3s_aa", "VJ_cdr3s_aa"), clonotypes.col = "clonotype_id_10x", split.by = "sample_id", stringDist.method = "levenshtein")
head(variants_agedCNS[[1]])
## split_group clonotype_id n sequences n variants mean variants distance
## 1 s1 clonotype7 1 1 NA
## 2 s1 clonotype58 1 1 NA
## 3 s1 clonotype65 1 1 NA
## 4 s1 clonotype43 1 1 NA
## 5 s1 clonotype33 1 1 NA
## 6 s1 clonotype49 1 1 NA
## max variants distance min variants distance most frequent variant
## 1 NA NA CARIGYAMDYW/CQQGNTLPPTF
## 2 NA NA CARDFTTVVARGYFDVW/CQQDYSSPWTF
## 3 NA NA CARGITTVVAYYYAMDYW/CLQYDNLYTF
## 4 NA NA CTTGPYDYYAMDYW/CQQHYSTPYTF
## 5 NA NA CARPHDYDGVDYW/CSQSTHVPPWTF
## 6 NA NA CARRYYSNYAWFAYW/CQQWSSYPYTF
## n VDJ cgene variants most frequent VDJ cgene
## 1 1 IGHM
## 2 1 IGHD
## 3 1 IGHD
## 4 1 IGHD
## 5 1 IGHM
## 6 1 IGHM
#set split.by to "none" if clonotyping was conducted across all samples
5.10 Searching for overlap between repertoires
Public clones or sequences are sequences that are shared between
individuals. The VDJ_overlap_heatmap function can quantify this overlap
and return a list of overlapping items
#overlap of VDJ V genes
VDJv_overlap <- VDJ_overlap_heatmap(VDJ = VGM[[1]],feature.columns = c("VDJ_vgene") ,grouping.column = "sample_id", axis.label.size = 20, pvalues.label.size = 12, add.barcode.table = T, plot.type = "ggplot")
head(VDJv_overlap[[2]]) #summary of overlap
## V1 V2 overlap
## 1 s1 s2 30
## 2 s1 s3 33
## 3 s1 s4 38
## 4 s2 s3 40
## 5 s2 s4 47
## 6 s3 s4 53
## items.overlapping
## 1 IGHV8-8;IGHV1-9;IGHV1-26;IGHV5-17;IGHV5-16;IGHV1-54;IGHV1-72;IGHV1-64;IGHV5-9-1;IGHV1-19;IGHV1-18;IGHV8-5;IGHV3-6;IGHV1-52;IGHV1-80;IGHV1-74;IGHV1-39;IGHV1-22;IGHV10-1;IGHV1-62-2;IGHV9-3;IGHV1-53;IGHV1-69;IGHV9-1;IGHV1-76;IGHV1-50;IGHV1-81;IGHV1-12;IGHV1-82;IGHV1-85
## 2 IGHV8-8;IGHV1-9;IGHV1-26;IGHV5-17;IGHV5-16;IGHV1-54;IGHV1-72;IGHV1-64;IGHV5-9-1;IGHV1-71;IGHV1-19;IGHV1-18;IGHV1-59;IGHV8-5;IGHV3-6;IGHV1-52;IGHV1-80;IGHV1-74;IGHV1-39;IGHV10-1;IGHV1-62-2;IGHV9-3;IGHV2-5;IGHV1-53;IGHV1-69;IGHV1-76;IGHV1-81;IGHV1-12;IGHV1-66;IGHV14-2;IGHV1-4;IGHV1-82;IGHV1-85
## 3 IGHV8-8;IGHV1-9;IGHV1-26;IGHV14-1;IGHV5-17;IGHV5-16;IGHV1-54;IGHV1-72;IGHV1-64;IGHV5-9-1;IGHV2-3;IGHV1-71;IGHV1-19;IGHV1-18;IGHV1-75;IGHV1-59;IGHV4-1;IGHV3-6;IGHV1-52;IGHV1-80;IGHV1-74;IGHV1-39;IGHV1-22;IGHV10-1;IGHV1-62-2;IGHV9-3;IGHV2-5;IGHV1-53;IGHV1-69;IGHV9-1;IGHV1-76;IGHV1-50;IGHV1-81;IGHV1-12;IGHV1-66;IGHV14-2;IGHV1-82;IGHV1-85
## 4 IGHV5-9-1;IGHV5-16;IGHV9-3;IGHV1-12;IGHV1-39;IGHV14-3;IGHV1-9;IGHV1-82;IGHV1-26;IGHV10-3;IGHV7-3;IGHV1-19;IGHV1-74;IGHV1-76;IGHV3-6;IGHV1-53;IGHV1-69;IGHV1-85;IGHV1-81;IGHV1-18;IGHV1-54;IGHV8-8;IGHV1-15;IGHV2-9-1;IGHV1-64;IGHV1-52;IGHV5-12;IGHV10-1;IGHV1-78;IGHV5-4;IGHV6-3;IGHV6-6;IGHV5-17;IGHV1-72;IGHV1-55;IGHV1-80;IGHV8-5;IGHV5-9;IGHV2-2;IGHV1-62-2
## 5 IGHV5-9-1;IGHV5-16;IGHV9-3;IGHV1-50;IGHV1-12;IGHV1-39;IGHV14-3;IGHV1-2;IGHV8-12;IGHV1-9;IGHV1-82;IGHV1-58;IGHV1-36;IGHV1-26;IGHV10-3;IGHV7-3;IGHV1-19;IGHV1-74;IGHV1-76;IGHV3-6;IGHV1-53;IGHV1-69;IGHV1-85;IGHV1-81;IGHV1-18;IGHV1-54;IGHV8-8;IGHV1-15;IGHV2-9-1;IGHV1-22;IGHV1-64;IGHV1-52;IGHV1-42;IGHV5-12;IGHV10-1;IGHV1-78;IGHV5-4;IGHV6-3;IGHV6-6;IGHV5-17;IGHV9-1;IGHV1-72;IGHV1-55;IGHV1-80;IGHV2-2;IGHV2-6;IGHV1-62-2
## 6 IGHV1-55;IGHV1-9;IGHV1-71;IGHV1-64;IGHV1-53;IGHV1-82;IGHV7-1;IGHV1-18;IGHV5-4;IGHV1-52;IGHV1-26;IGHV1-7;IGHV5-9-1;IGHV2-9-1;IGHV5-12;IGHV13-2;IGHV1-72;IGHV1-69;IGHV1-19;IGHV10-3;IGHV1-81;IGHV1-80;IGHV6-3;IGHV1-5;IGHV1-15;IGHV6-6;IGHV1-84;IGHV5-16;IGHV9-3;IGHV1-78;IGHV1-66;IGHV2-5;IGHV3-6;IGHV1-12;IGHV1-39;IGHV7-3;IGHV1-62-2;IGHV1-74;IGHV10-1;IGHV1-59;IGHV8-8;IGHV1-77;IGHV15-2;IGHV14-3;IGHV14-2;IGHV5-17;IGHV1-61;IGHV3-1;IGHV1-54;IGHV2-2;IGHV1-63;IGHV1-76;IGHV1-85
## overlap_lab
## 1 30
## 2 33
## 3 38
## 4 40
## 5 47
## 6 53
#overlap of clones
clone_overlap <- VDJ_overlap_heatmap(VDJ = VGM[[1]],feature.columns = c("VDJ_cdr3s_aa","VJ_cdr3s_aa") ,grouping.column = "sample_id", axis.label.size = 20, pvalues.label.size = 12, add.barcode.table = T, plot.type = "ggplot")
print(clone_overlap[[1]])
#Pheatmap plots will function only with length(unique(grouping.column)) > 4
Given the small size of the two repertoires, it is unsurprising that
no clones are shared
