Section 1: running seismic on the included sample data
+seismicGWAS includes a preloaded small single cell RNAseq dataset in the form of a SingleCellExperiment object tmfacs_sce_small
tmfacs_sce_small
+#> class: SingleCellExperiment
+#> dim: 23341 4000
+#> metadata(0):
+#> assays(1): logcounts
+#> rownames(23341): 0610005C13Rik 0610007C21Rik ... l7Rn6
+#> zsGreen-transgene
+#> rowData names(1): symbol
+#> colnames(4000): M4.MAA000844.3_10_M.1.1 A18.MAA000388.3_11_M.1.1 ...
+#> M5.MAA000910.3_10_M.1.1 L2.MAA000563.3_10_M.1.1
+#> colData names(17): nReads orig.ident ... sizeFactor cluster_name
+#> reducedDimNames(0):
+#> mainExpName: RNA
+#> altExpNames(0):and two sample MAGMA files, one for Rheumatoid Arthritis (Okada et al. 2013)
+head(ra_magma)
+#> GENE CHR START STOP NSNPS NPARAM N ZSTAT P
+#> 1 148398 1 824993 889961 1 1 58284 -0.025069 0.51000
+#> 2 339451 1 860967 911099 1 1 58284 -0.025069 0.51000
+#> 3 84069 1 866872 920488 12 2 58284 -0.671780 0.74914
+#> 4 26155 1 869583 929679 36 4 58284 -1.717100 0.95702
+#> 5 84808 1 900579 952473 80 6 58284 -1.124900 0.86968
+#> 6 9636 1 913847 959920 93 7 58284 -0.991290 0.83923and one for Type 2 Diabetes (Xue et al. 2018)
+head(t2d_magma)
+#> GENE CHR START STOP NSNPS NPARAM N ZSTAT P
+#> 1 254173 1 1074286 1143315 21 3 572082 1.75840 0.039337
+#> 2 8784 1 1128888 1177163 6 1 572870 1.93210 0.026671
+#> 3 55210 1 1412523 1480067 9 1 569925 0.24436 0.403470
+#> 4 339453 1 1460158 1510740 37 2 571069 0.14879 0.440860
+#> 5 29101 1 1467053 1545262 39 2 571016 0.15608 0.437990
+#> 6 9906 1 1646277 1712438 16 2 576212 -0.48124 0.684830Calculating cell type specificity scores
+The first step of seismic is to calculate cell type specificity scores which will be used for subsequent cell type-trait association calculations and influential gene analyses. Here we provide the sample SingleCellExperiment object and column in the data containing cell type labels cluster_name.
tmfacs_sscore <- calc_specificity(tmfacs_sce_small, ct_label_col='cluster_name')
+head(tmfacs_sscore[, 1:3]) # show a truncated view of the matrix
+#> Bladder.Bladder mesenchymal cell
+#> 0610005C13Rik 0.00000000
+#> 0610007C21Rik 0.05171198
+#> 0610007L01Rik 0.05053099
+#> 0610007N19Rik 0.10618537
+#> 0610007P08Rik 0.04881961
+#> 0610007P14Rik 0.04976406
+#> Bladder.Luminal bladder epithelial cell
+#> 0610005C13Rik 0.000000e+00
+#> 0610007C21Rik 5.416472e-02
+#> 0610007L01Rik 5.907819e-02
+#> 0610007N19Rik 4.726793e-05
+#> 0610007P08Rik 6.890189e-02
+#> 0610007P14Rik 6.196749e-02
+#> Brain_Myeloid.microglial cell
+#> 0610005C13Rik 3.639054e-19
+#> 0610007C21Rik 3.768650e-02
+#> 0610007L01Rik 3.409995e-07
+#> 0610007N19Rik 4.422689e-185
+#> 0610007P08Rik 2.758088e-08
+#> 0610007P14Rik 2.225707e-17The original single cell experiment in the tmfacs_sce_small was derived from Tabula Muris. To associate with human GWAS traits we need to convert the gene identifiers. For convenience seismicGWAS has a function for translating between gene ids for a specificity score matrix. We change here from mouse gene symbols to human entrez ids to match the MAGMA data.
tmfacs_sscore_hsa <- translate_gene_ids(tmfacs_sscore, from='mmu_symbol')
+head(tmfacs_sscore_hsa[,1:3]) # show a truncated view of the matrix
+#> 6 x 3 Matrix of class "dgeMatrix"
+#> Bladder.Bladder mesenchymal cell Bladder.Luminal bladder epithelial cell
+#> 81932 2.542912e-04 0.01671191
+#> 81577 5.003187e-04 0.01487892
+#> 8614 9.234388e-02 0.02736780
+#> 256281 3.213228e-07 0.08150739
+#> 221687 0.000000e+00 0.00000000
+#> 346389 0.000000e+00 0.17272120
+#> Brain_Myeloid.microglial cell
+#> 81932 7.503793e-12
+#> 81577 3.645070e-10
+#> 8614 0.000000e+00
+#> 256281 7.970981e-07
+#> 221687 0.000000e+00
+#> 346389 0.000000e+00Calculating cell type trait associations
+We can now get cell type-trait associations for the RA and Type 2 Diabetes data and plot the top 20 associations per trait. Default plots show FDR corrected p-values plotted on a -log scale so that stronger associations appear at the top of the boxplot.
+# calculate and plot the cell type-trait associations for RA
+ra <- get_ct_trait_associations(tmfacs_sscore_hsa, ra_magma)
+plot_top_associations(ra, limit = 20)
# calculate and plot the cell type-trait associations for T2D
+t2d <- get_ct_trait_associations(tmfacs_sscore_hsa, t2d_magma)
+plot_top_associations(t2d, limit = 20)
Influential gene analysis
+With seismicGWAS we can also look at the genes driving a particular cell type-trait association, in a process we call influential gene analysis. Here, we need to select a particular cell type (in our example we choose pancreas beta cell) and a trait (we use Type 2 Diabetes again). It is important to note that influential gene calculations are only meaningful for cell type-trait associations that are significant. So we have chosen a cell type and trait that were found to be significant in the cell type-trait associations above.
+# calculate influential genes
+ct <- "Pancreas.beta cell"
+t2d_beta_inf_genes <- find_inf_genes(ct, tmfacs_sscore_hsa, t2d_magma)
+head(t2d_beta_inf_genes)
+#> gene specificity zstat dfbetas is_influential
+#> <char> <num> <num> <num> <lgcl>
+#> 1: 3767 0.4410499 6.1866 0.1941221 TRUE
+#> 2: 169026 0.4932408 5.5197 0.1886886 TRUE
+#> 3: 2645 0.4184336 5.2990 0.1513507 TRUE
+#> 4: 3651 0.8982126 3.2276 0.1386102 TRUE
+#> 5: 6833 0.4991599 4.1322 0.1293991 TRUE
+#> 6: 7466 0.2439394 7.5000 0.1256267 TRUENote: Only genes with positive contributions to the cell type-trait association are marked influential by default. While negative contributions may be informative, they are not as directly interpretable and are driven by genes from other cell types or are somewhat cell type specific but not associated with the given trait.
+Next, to aid in result interpretation we make use of the biomaRt package to convert entrez ids to gene symbols.
# convert ids to symbols using biomaRt
+library(biomaRt)
+
+ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
+hsa.map <- getBM(attributes = c('entrezgene_id', 'hgnc_symbol'),
+ filters = 'entrezgene_id',
+ values = t2d_beta_inf_genes$gene,
+ mart = ensembl)
+
+hsa.map$entrezgene_id <- as.character(hsa.map$entrezgene_id)
+t2d_beta_inf_genes <- merge(t2d_beta_inf_genes, hsa.map,
+ by.x='gene', by.y='entrezgene_id')We can also use the gene symbols to label the influential points in the influential gene plot.
+# plot influential genes
+plot_inf_genes(t2d_beta_inf_genes, gene_col = 'hgnc_symbol', num_labels = 5)
