README.md
In PavlidisLab/markerGeneProfile: Cell type marker gene profiles in whole tissue samples
markerGeneProfile
This package includes functions responsible for marker gene selection
and marker gene profile estimation estimation as described in Mancarci
et al. 2017. It
also includes a copy of mouse brain cell type markers from the
neuroExpressoAnalysis
package for convenience along with mock data for easy testing.
Installation
Use devtools to install. Additional github packages needs to be
installed for it work.
devtools::install_github('oganm/markerGeneProfile')
In this document additional packages are used that are not package
dependencies
install.packages('ggplot2')
install.packages('gplots')
install.packages('viridis')
install.packages('dplyr')
install.packages('knitr')
Usage
Marker genes
A list of marker genes specific to or enriched in a cell type is
required in order to estimate cell type profiles. In this package, a
copy of mouse brain cell type-specific markers from
neuroExpressoAnalysis,
the package that summarizes the entire analysis performed in Mancarci et
al. 2017 is included (mouseMarkerGenes). If an different marker gene
set is needed, steps outlined below can be followed to create one from a
new dataset.
Sample data for marker gene selection
This package includes a sample cell type-specific transcriptomic dataset
representing expression profiles from multiple purified cell types aimed
to demonstrate the minimal information required for the selection of
marker genes.
mgp_sampleProfilesMeta includes the basic metadata required for the
cell type specific expression dataset.
data(mgp_sampleProfilesMeta)
knitr::kable(head(mgp_sampleProfilesMeta))
| sampleName | replicate | PMID | CellType | region | RegionToParent | RegionToChildren |
| :--------- | --------: | ---: | :------- | :------- | :------------- | :--------------- |
| Sample01 | 1 | 1 | Cell A | Region 1 | TRUE | TRUE |
| Sample02 | 1 | 1 | Cell A | Region 1 | TRUE | TRUE |
| Sample03 | 1 | 1 | Cell A | Region 1 | TRUE | TRUE |
| Sample04 | 2 | 2 | Cell A | Region 1 | TRUE | TRUE |
| Sample05 | 2 | 2 | Cell A | Region 1 | TRUE | TRUE |
| Sample06 | 2 | 2 | Cell A | Region 1 | TRUE | TRUE |
sampleName: name of the samples. This needs to correspond to column
names in the expression file.
replicate: A vector marking which samples are replicates of each
other.
PMID: A vector marking which samples come from the same study.
Normally taking PMIDs of the papers is a good idea.
CellType: A vector marking the cell types that the samples
represent.
region: The regions samples are extracted from. Only needed if
region specific genes are to be selected.
RegionToParent: If region specific genes are to be selected and a
region hierarchy is to be used, this column controls whether or not the
sample should be included in the parent regions of the indicated region.
If not provided it will default to TRUE. The name of this column is
hard coded and should not be changed.
RegionToChildren: Same as above except it controls if the sample
should be included in the children regions. If not provided it will
default to TRUE. The name of this column is hard coded and should not
be changed.
mgp_sampleProfiles is a sample expression data. Gene.Symbol column
is the gene identifier that should be composed of unique IDs while the
rest are sample names that corresponds to the relevant column in the
metadata file. Other columns can be present before the sample data but
they should not be of class double.
data(mgp_sampleProfiles)
knitr::kable(mgp_sampleProfiles)
| Gene.Symbol | Sample01 | Sample02 | Sample03 | Sample04 | Sample05 | Sample06 | Sample07 | Sample08 | Sample09 | Sample10 | Sample11 | Sample12 | Sample13 | Sample14 | Sample15 | Sample16 | Sample17 | Sample18 |
| :---------- | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: |
| Gene1 | 16 | 16 | 16 | 16 | 16 | 16 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| Gene2 | 1 | 1 | 1 | 1 | 1 | 1 | 16 | 16 | 16 | 16 | 16 | 16 | 1 | 1 | 1 | 1 | 1 | 1 |
| Gene3 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 16 | 16 | 16 | 16 | 16 | 16 |
| Gene4 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 13 | 13 | 13 | 13 | 13 | 13 |
| Gene5 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 9 | 9 | 9 | 9 | 9 | 9 |
| Gene6 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 7 | 7 | 7 | 7 | 7 | 7 |
mpg_sampleRegionHiearchy is a sample region hiearchy. It is a nested
named list.
data(mpg_sampleRegionHiearchy)
mpg_sampleRegionHiearchy
## $All
## $All$`Region 1`
## [1] ""
##
## $All$`Region 2`
## [1] ""
nametree(mpg_sampleRegionHiearchy)
## All
## ├──Region 1
## └──Region 2
In this example Region 1 and Region 2 are subsets of All region
Selection of marker genes
Marker gene selection is performed using three functions:
markerCandidates, pickMarkers and rotateSelect. By default
markerCandidates will return files for each cell type in a region that
lists the gene that are above a given silhouette and fold change
thresholds. Other variables are briefly explained below but see the
package documentation for in depth explanations
markerCandidates(design = mgp_sampleProfilesMeta, # the design file
expression = mgp_sampleProfiles, # expression file
outLoc = 'README_files/quickSelection', # output directory
groupNames = 'CellType', # name of the column with cell types. can be a vector
regionNames = 'region', # name of the column with brain regions. leave NULL if no region seperation is desired
PMID = 'PMID', # name of the column with study identifiers
sampleName = 'sampleName', # name of the column with sample names
replicates = 'replicate', # name of the column with replicates
foldChangeThresh = 10, # threshold of fold change for gene selection (default is 10)
minimumExpression = 8, # minimum expression level that a gene can be considered a marker gene (default is 8)
background = 6, # background level of expression (default is 6)
regionHierarchy = mpg_sampleRegionHiearchy, # hierarchy of brain regions to be used
geneID = 'Gene.Symbol', # column name with with gene idenditifers
cores = 8 # number of cores to use in parallelization
)
This creates 3 directories in the output directory
list.files('README_files/quickSelection')
## [1] "All_CellType" "CellType" "Region 2_CellType"
The CellType directory is a list of marker genes that disregards all
region specifications (redundant with All_CellType in this case) while
Region 2_CellType and All_CellType directories inlcude cell types
from the relevant region. Note the absence of Region 1_CellType since
that region only has a single cell
type.
read.table('README_files/quickSelection/All_CellType/Cell C') %>% knitr::kable()
| V1 | V2 | V3 | V4 |
| :---- | -: | -: | -: |
| Gene3 | 15 | 1 | 1 |
| Gene4 | 12 | 1 | 1 |
| Gene5 | 8 | 1 | 1 |
This file shows the candidate genes for cell type Cell C in region
All. The first column is the gene identifier, the second is change in
expression in log_2 scale and the third one is the silhouette
coefficient. Note that Gene6 is absent since its expression level was
below the minimum level allowed. markerCandidates function does not
apply a threshold for silhouette coefficient it also doesn’t check to
see if a gene satisfies fold change threshold for multiple genes.
pickMarkers function does that.
pickMarkers('README_files/quickSelection/All_CellType/',
foldChange = 1, # this is a fixed fold change threshold that ignores some leniency that comes from markerCandidates. setting it to 1 makes it irrelevant
silhouette = 0.5)
## $`Cell A`
## [1] "Gene1"
##
## $`Cell B`
## [1] "Gene2"
##
## $`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
If all genes for all regions needs to be seen
pickMarkersAll('README_files/quickSelection',
foldChange = 1,
silhouette = 0.5)
## $All_CellType
## $All_CellType$`Cell A`
## [1] "Gene1"
##
## $All_CellType$`Cell B`
## [1] "Gene2"
##
## $All_CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $CellType
## $CellType$`Cell A`
## [1] "Gene1"
##
## $CellType$`Cell B`
## [1] "Gene2"
##
## $CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $`Region 2_CellType`
## $`Region 2_CellType`$`Cell B`
## [1] "Gene2"
##
## $`Region 2_CellType`$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
Better selection of marker genes
The above method is a quick way to pick markers but it does not handle
bimodality in expression distribution well. To ensure robustness of the
results it is better to perform multiple selections with permutations.
markerCandidates function has variables to handle permutations for
you. rotate controls what is the percentage of samples that should be
removed every time. seed controls the random seed and is there to ensure
reproducibility.
for (i in 1:10){
markerCandidates(design = mgp_sampleProfilesMeta, # the design file
expression = mgp_sampleProfiles, # expression file
outLoc = file.path('README_files/Rotation',i), # output directory
groupNames = 'CellType', # name of the column with cell types. can be a vector
regionNames = 'region', # name of the column with brain regions. leave NULL if no region seperation is desired
PMID = 'PMID', # name of the column with study identifiers
sampleName = 'sampleName', # name of the column with sample names
replicates = 'replicate', # name of the column with replicates
foldChangeThresh = 10, # threshold of fold change for gene selection (default is 10)
minimumExpression = 8, # minimum expression level that a gene can be considered a marker gene (default is 8)
background = 6, # background level of expression (default is 6)
regionHierarchy = mpg_sampleRegionHiearchy, # hierarchy of brain regions to be used
geneID = 'Gene.Symbol', # column name with with gene idenditifers
cores = 8, # number of cores to use in parallelization
rotate = 0.33,
seed = i
)
}
This creates multiple selection directories. rotateSelect can be used
to count the number of times a gene is selected for each cell type in
each region. This creates another directory similar to the output of
markerCandidates. Again, valid markers can be acquired using
pickMarkers
rotateSelect(rotationOut='README_files/Rotation',
rotSelOut='README_files/RotSel',
cores = 8,
foldChange = 1 # this is a fixed fold change threshold that ignores some leniency that comes from markerCandidates. setting it to 1 makes it irrelevant
)
pickMarkers('README_files/RotSel/All_CellType/',rotationThresh = 0.95)
## $`Cell A`
## [1] "Gene1"
##
## $`Cell B`
## [1] "Gene2"
##
## $`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
pickMarkersAll('README_files/RotSel',rotationThresh = 0.95)
## $All_CellType
## $All_CellType$`Cell A`
## [1] "Gene1"
##
## $All_CellType$`Cell B`
## [1] "Gene2"
##
## $All_CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $CellType
## $CellType$`Cell A`
## [1] "Gene1"
##
## $CellType$`Cell B`
## [1] "Gene2"
##
## $CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $`Region 2_CellType`
## $`Region 2_CellType`$`Cell B`
## [1] "Gene2"
##
## $`Region 2_CellType`$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
Marker gene profiles (MGP)
###Sample data for MGP estimation
The package includes mouse brain cell type markers published in Mancarci
et al. 2017 and gene expression data from substantia nigra samples from
healthy donors and Parkinson’s disease patients by Lesnick et
al. 2007.
Mouse marker genes is available in mouseMarkerGenes object as a nested
list.
data(mouseMarkerGenes)
names(mouseMarkerGenes)
## [1] "All" "Amygdala" "BasalForebrain"
## [4] "Brainstem" "Cerebellum" "Cerebrum"
## [7] "Cortex" "Hippocampus" "LocusCoeruleus"
## [10] "Midbrain" "SpinalCord" "Striatum"
## [13] "Subependymal" "SubstantiaNigra" "Thalamus"
lapply(mouseMarkerGenes$Midbrain[1:3],head, 14)
## $Astrocyte
## [1] "Aass" "Acsbg1" "Acsl6" "Acss1" "Add3" "Adhfe1"
## [7] "AI464131" "Aldh1l1" "Aldh6a1" "Antxr1" "Aox1" "Apoe"
## [13] "Aqp4" "Axl"
##
## $BrainstemCholin
## [1] "2310030G06Rik" "Anxa2" "Cabp1" "Calca"
## [5] "Calcb" "Cd24a" "Cd55" "Cda"
## [9] "Chodl" "Ecel1" "Fxyd7" "Hebp2"
## [13] "Hspb1" "Hspb8"
##
## $Dopaminergic
## [1] "Cacna2d2" "Cadps2" "Chrna6" "Mapk8ip2" "Nr4a2" "Ntn1"
## [7] "Prkcg" "Rian" "Scn2a" "Slc6a3" "Snhg11" "Tenm1"
## [13] "Th" "Zim3"
Available Lesnick et al. data is stored in mgp_LesnickParkinsonsExp
and mgp_LesnickParkinsonsMeta objects
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:testthat':
##
## matches
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
data(mgp_LesnickParkinsonsExp)
mgp_LesnickParkinsonsExp %>%
dplyr::select(-GeneNames) %>%
head %>% {.[,1:6]}
## Probe Gene.Symbol NCBIids GSM184354.cel GSM184355.cel
## 3366 1007_s_at DDR1 780 10.236880 9.891552
## 44099 1053_at RFC2 5982 5.421790 5.280541
## 14880 117_at HSPA7|HSPA6 3311|3310 5.164445 4.651754
## 11594 121_at PAX8 7849 7.076004 7.035090
## 11907 1255_g_at GUCA1A 2978 3.107388 3.418976
## 2069 1294_at UBA7 7318 6.644858 6.182664
## GSM184356.cel
## 3366 10.498371
## 44099 5.852467
## 14880 4.729189
## 11594 6.698765
## 11907 3.491832
## 2069 5.982642
data(mgp_LesnickParkinsonsMeta)
mgp_LesnickParkinsonsMeta %>% head
## GSM disease
## 1 GSM184354 Control
## 2 GSM184355 Control
## 3 GSM184356 Control
## 4 GSM184357 Control
## 5 GSM184358 Control
## 6 GSM184359 Control
Before MGP estimation, it is important to filter expression data of low
expressed genes and make sure all genes are represented only once in the
dataset. While there are many probeset summarization and methods, for
this work we chose the most variable probeset and remove all probes with
a median expression below the median expression of the dataset.
unfilteredParkinsonsExp = mgp_LesnickParkinsonsExp # keep this for later
medExp = mgp_LesnickParkinsonsExp %>%
sepExpr() %>% {.[[2]]} %>%
unlist %>% median
# mostVariable function is part of this package that does probe selection and filtering for you
mgp_LesnickParkinsonsExp = mostVariable(mgp_LesnickParkinsonsExp,
threshold = medExp,
threshFun= median)
MGP estimation
MGPs are simply the first principal component of marker gene expression.
This method of marker gene profile estimation is similar to the
methodology of multiple previous works that aim to estimate relative
abundance of cell types in a whole tissue sample (Chikina et al., 2015;
Westra et al., 2015; Xu et al., 2013).
Primary function that deals with MGP estimation is mgpEstimate. This
function will take in expression data and marker gene lists to return
MGPs. exprData is the expression matrix which should include gene
identifiers as a column. Other columns can be present at the beginning
of the data frame but should not be of class double. genes is the
list of markers. It is assumed to be a list of character vectors, each
vector containing gene names for a cell type. geneColName is the name
of the column where gene names are found. geneTransform is a function
that will be applied to the gene names provided in genes. Note that
this by default tranforms mouse gene names to human gene names. Set it
to NULL if this is not desired.
Below a basic estimation performed with other important variables
briefly explained.
estimations = mgpEstimate(exprData=mgp_LesnickParkinsonsExp,
genes=mouseMarkerGenes$Midbrain,
geneColName='Gene.Symbol',
outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
geneTransform =function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
groups=mgp_LesnickParkinsonsMeta$disease, #if there are experimental groups provide them here. if not desired set to NULL
seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
removeMinority = TRUE) # removes genes if they are the minority in terms of rotation sign from estimation process
Dopamine rgic cell loss is a known effect of Parkinson’s Disease. To see
if this effect can be observed we can look at dopaminergic MGPs in
healthy donors vs Parkinson’s Disease patients
library(ggplot2)
ls(estimations)
## [1] "allMarkerExpression" "estimates"
## [3] "fullPCAs" "groups"
## [5] "meanUsedMarkerExpression" "removedMarkerRatios"
## [7] "rotations" "simpleScaledEstimation"
## [9] "trimmedPCAs" "usedMarkerExpression"
ls(estimations$estimates)
## [1] "Astrocyte" "BrainstemCholin"
## [3] "Dopaminergic" "Microglia"
## [5] "Microglia_activation" "Microglia_deactivation"
## [7] "Oligo" "Serotonergic"
dopaminergicFrame =
data.frame(`Dopaminergic MGP` = estimations$estimates$Dopaminergic,
state = estimations$groups$Dopaminergic, # note that unless outlierSampleRemove is TRUE this will be always the same as the groups input
check.names=FALSE)
ggplot2::ggplot(dopaminergicFrame,
aes(x = state, y = `Dopaminergic MGP`)) +
geom_boxvio() + geom_jitter(width = .05) # this is just a convenience function that outputs a list of ggplot elements.
To see if the difference between the two groups is statistically
significant we use wilcoxon test (Mann-Whitney test). Wilcoxon-test is a
non parametric tests that does not make assumptions about the
distribution of the data. We chose to use it because marker gene
profiles are not normally
distributed.
group1 = estimations$estimates$Dopaminergic[estimations$groups$Dopaminergic %in% "Control"]
group2 = estimations$estimates$Dopaminergic[estimations$groups$Dopaminergic %in% "PD"]
wilcox.test(group1,group2)
##
## Wilcoxon rank sum test
##
## data: group1 and group2
## W = 119, p-value = 0.006547
## alternative hypothesis: true location shift is not equal to 0
Based on these results we can say that there is indeed a significant
difference between doparminergic marker gene profiles between control
and parkinson’s disease patients.
This difference can be also observed by looking at gene expression of
markers
estimations$usedMarkerExpression$Dopaminergic%>%
as.matrix %>%
gplots::heatmap.2(trace = 'none',
scale='row',Rowv = FALSE,Colv = FALSE, dendrogram = 'none',
col= viridis::viridis(10),cexRow = 1.5, cexCol = 1,
ColSideColors = estimations$groups$Dopaminergic %>%
toColor(palette = c('Control' = 'blue',
"PD" = "red")) %$% cols ,
margins = c(7,8))
Indeed in general most marker genes have a high expression is control
samples which are shown in green.
It is important to note that not all marker genes are used in the
estimation of marker gene profiles. Below we see all genes that are
identified as dopaminergic cell markers with human orthologues. Yet not
all these genes appear in the heatmap
above
mouseHumanGeneTable = mouseMarkerGenes$Midbrain$Dopaminergic %>% homologene::mouse2human()
allHumanDopaGenes = mouseHumanGeneTable %$% humanGene
mouseHumanGeneTable
## mouseGene humanGene mouseID humanID
## 1 Cacna2d2 CACNA2D2 56808 9254
## 2 Cadps2 CADPS2 320405 93664
## 3 Chrna6 CHRNA6 11440 8973
## 4 Mapk8ip2 MAPK8IP2 60597 23542
## 5 Nr4a2 NR4A2 18227 4929
## 6 Ntn1 NTN1 18208 9423
## 7 Prkcg PRKCG 18752 5582
## 8 Slc6a3 SLC6A3 13162 6531
## 9 Tenm1 TENM1 23963 10178
## 10 Th TH 21823 7054
Some of these genes are removed because they are not included in our
dataset, either because they are not in the platform, or because the
gene was filtered in the pre-processing stage due to low
expression.
allHumanDopaGenes[!allHumanDopaGenes %in% mgp_LesnickParkinsonsExp$Gene.Symbol]
## [1] "CHRNA6"
allGenes = allHumanDopaGenes[allHumanDopaGenes %in% mgp_LesnickParkinsonsExp$Gene.Symbol]
Other genes can be removed because mgpEstimate thinks they don’t
correlate well with the rest of the genes. Details of this process can
be found in Mancarci et al. 2017
manuscript
allGenes[!allGenes %in% rownames(estimations$usedMarkerExpression$Dopaminergic)]
## [1] "PRKCG"
genesUsed = rownames(estimations$usedMarkerExpression$Dopaminergic)
By looking at the unfiltered expression data, we can see how these two
genes were unsuitable for MGP estimation
toPlot =
unfilteredParkinsonsExp[unfilteredParkinsonsExp$Gene.Symbol %in%
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic)$humanGene,] %>%
mostVariable(threshold = 0)
toPlot %>%
mostVariable(threshold = 0) %>%
sepExpr() %>%
{.[[2]]} %>% as.matrix() %>%
apply(1,function(x){scale01(x)}) %>%
t %>%
{rownames(.) =
toPlot$Gene.Symbol[toPlot$Gene.Symbol %in%
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic)$humanGene];.} %>%
reshape2::melt() %>%
{colnames(.) = c('Gene','Sample','Expression');.} %>%
dplyr::mutate(`Is used?` = rep('used',length(Gene)) %>%
{.[Gene %in% 'CHRNA6'] = '(CHRNA6) not expressed';.[Gene %in% 'PRKCG'] = '(PRKCG) not correlated';.}) %>%
ggplot(aes(y = Expression, x = Sample, group = Gene, color = `Is used?`)) +
geom_line() +
cowplot::theme_cowplot() +
theme( axis.text.x= element_blank(),
axis.title = element_text(size = 20),
legend.text = element_text(size = 17),
legend.title = element_text(size=22),
plot.title = element_text(size = 20)
) +
ggtitle('Scaled expression of markers')
toPlot %>%
mostVariable(threshold = 0) %>%
sepExpr() %>%
{.[[2]]} %>% as.matrix()%>%
{rownames(.) =
toPlot$Gene.Symbol[toPlot$Gene.Symbol %in%
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic)$humanGene];.} %>%
reshape2::melt() %>%
{colnames(.) = c('Gene','Sample','Expression');.} %>%
dplyr::mutate(`Is used?` = rep('used',length(Gene)) %>%
{.[Gene %in% 'CHRNA6'] = 'CHRNA6 - not expressed';.[Gene %in% 'PRKCG'] = 'PRKCG - not correlated';.}) %>%
ggplot(aes(y = Expression, x = Sample, group = Gene, color = `Is used?`)) +
geom_line() +
cowplot::theme_cowplot() +
theme( axis.text.x= element_blank(),
axis.title = element_text(size = 20),
legend.text = element_text(size = 17),
legend.title = element_text(size=22),
plot.title = element_text(size = 20)
)+
ggtitle('Nonscaled expression of markers')
Both of these genes seem to be non-expressed though PRKCG managed to be
just above the removal threshold. Luckily lack of correlation reveals
that it is not behaving as other marker genes.
The ratio of genesUsed and allGenes (all markers available in the
study) can be used as a confidence metric. If a significant portion of
the genes do not correlate well with each other that may point to
presence of non cell type specific signal (regulation or noise). For all
cell types this ratio is outputted. You’ll see a warning if this ratio
ever exceeds
0.4
estimations$removedMarkerRatios
## Astrocyte BrainstemCholin Dopaminergic
## 0.12264151 0.35294118 0.11111111
## Microglia Microglia_activation Microglia_deactivation
## 0.09734513 0.08333333 0.12149533
## Oligo Serotonergic
## 0.05882353 0.00000000
The ratio for dopaminergic cells is fairly low. We can also look at the
amount of varience explained in the first PC of marker gene expression.
If this value is low, it can point to higher amounts of confounding
noise or regulation.
estimations$trimmedPCAs$Dopaminergic %>% summary()
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.3056 1.0829 0.87149 0.53287 0.45796 0.38075
## Proportion of Variance 0.6645 0.1466 0.09494 0.03549 0.02622 0.01812
## Cumulative Proportion 0.6645 0.8110 0.90597 0.94147 0.96768 0.98580
## PC7 PC8
## Standard deviation 0.28376 0.18177
## Proportion of Variance 0.01006 0.00413
## Cumulative Proportion 0.99587 1.00000
For instance unlike dopaminergic cells, cholinergic cells seem to have a
higher proportion of their marker genes removed. Looking at the
variation explained by first PC reveals that first PC only explains 28%
of all variance.
estimations$trimmedPCAs$BrainstemCholin %>% summary()
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7194 1.4430 1.2543 1.1992 0.8893 0.82035 0.77928
## Proportion of Variance 0.2688 0.1893 0.1430 0.1307 0.0719 0.06118 0.05521
## Cumulative Proportion 0.2688 0.4581 0.6011 0.7318 0.8037 0.86491 0.92012
## PC8 PC9 PC10 PC11
## Standard deviation 0.59534 0.48891 0.44505 0.29527
## Proportion of Variance 0.03222 0.02173 0.01801 0.00793
## Cumulative Proportion 0.95234 0.97407 0.99207 1.00000
Such a result can occur if
- Marker genes are highly regulated.
- There is little difference in cell type proportions or wholescale
regulation of marker genes and MGP is driven by noise and other
confounds.
Note about homologene
In this tutorial homologene database is used through
homologene
package. Note that the default version of this database is somewhat old
which means some annotations may have been changed. The homologene
package also includes an updated version of the homologene package to
make sure you get all the correct matches.
If we look at dopaminergic cell markers for instance we could have
gotten an extra gene (SCN2A) if we used the updated
version.
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic) %$% humanGene
## [1] "CACNA2D2" "CADPS2" "CHRNA6" "MAPK8IP2" "NR4A2" "NTN1"
## [7] "PRKCG" "SLC6A3" "TENM1" "TH"
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic,
db = homologene::homologeneData2) %$% humanGene
## [1] "CACNA2D2" "CADPS2" "CHRNA6" "MAPK8IP2" "NR4A2" "NTN1"
## [7] "PRKCG" "SCN2A" "SLC6A3" "TENM1" "TH"
PavlidisLab/markerGeneProfile documentation built on July 11, 2019, 11:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
markerGeneProfile
This package includes functions responsible for marker gene selection and marker gene profile estimation estimation as described in Mancarci et al. 2017. It also includes a copy of mouse brain cell type markers from the neuroExpressoAnalysis package for convenience along with mock data for easy testing.
Installation
Use devtools to install. Additional github packages needs to be installed for it work.
devtools::install_github('oganm/markerGeneProfile')
In this document additional packages are used that are not package dependencies
install.packages('ggplot2')
install.packages('gplots')
install.packages('viridis')
install.packages('dplyr')
install.packages('knitr')
Usage
Marker genes
A list of marker genes specific to or enriched in a cell type is
required in order to estimate cell type profiles. In this package, a
copy of mouse brain cell type-specific markers from
neuroExpressoAnalysis,
the package that summarizes the entire analysis performed in Mancarci et
al. 2017 is included (mouseMarkerGenes). If an different marker gene
set is needed, steps outlined below can be followed to create one from a
new dataset.
Sample data for marker gene selection
This package includes a sample cell type-specific transcriptomic dataset representing expression profiles from multiple purified cell types aimed to demonstrate the minimal information required for the selection of marker genes.
mgp_sampleProfilesMeta includes the basic metadata required for the
cell type specific expression dataset.
data(mgp_sampleProfilesMeta)
knitr::kable(head(mgp_sampleProfilesMeta))
| sampleName | replicate | PMID | CellType | region | RegionToParent | RegionToChildren | | :--------- | --------: | ---: | :------- | :------- | :------------- | :--------------- | | Sample01 | 1 | 1 | Cell A | Region 1 | TRUE | TRUE | | Sample02 | 1 | 1 | Cell A | Region 1 | TRUE | TRUE | | Sample03 | 1 | 1 | Cell A | Region 1 | TRUE | TRUE | | Sample04 | 2 | 2 | Cell A | Region 1 | TRUE | TRUE | | Sample05 | 2 | 2 | Cell A | Region 1 | TRUE | TRUE | | Sample06 | 2 | 2 | Cell A | Region 1 | TRUE | TRUE |
sampleName: name of the samples. This needs to correspond to column names in the expression file.
replicate: A vector marking which samples are replicates of each other.
PMID: A vector marking which samples come from the same study. Normally taking PMIDs of the papers is a good idea.
CellType: A vector marking the cell types that the samples represent.
region: The regions samples are extracted from. Only needed if region specific genes are to be selected.
RegionToParent: If region specific genes are to be selected and a
region hierarchy is to be used, this column controls whether or not the
sample should be included in the parent regions of the indicated region.
If not provided it will default to TRUE. The name of this column is
hard coded and should not be changed.
RegionToChildren: Same as above except it controls if the sample
should be included in the children regions. If not provided it will
default to TRUE. The name of this column is hard coded and should not
be changed.
mgp_sampleProfiles is a sample expression data. Gene.Symbol column
is the gene identifier that should be composed of unique IDs while the
rest are sample names that corresponds to the relevant column in the
metadata file. Other columns can be present before the sample data but
they should not be of class double.
data(mgp_sampleProfiles)
knitr::kable(mgp_sampleProfiles)
| Gene.Symbol | Sample01 | Sample02 | Sample03 | Sample04 | Sample05 | Sample06 | Sample07 | Sample08 | Sample09 | Sample10 | Sample11 | Sample12 | Sample13 | Sample14 | Sample15 | Sample16 | Sample17 | Sample18 | | :---------- | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | -------: | | Gene1 | 16 | 16 | 16 | 16 | 16 | 16 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | | Gene2 | 1 | 1 | 1 | 1 | 1 | 1 | 16 | 16 | 16 | 16 | 16 | 16 | 1 | 1 | 1 | 1 | 1 | 1 | | Gene3 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 16 | 16 | 16 | 16 | 16 | 16 | | Gene4 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 13 | 13 | 13 | 13 | 13 | 13 | | Gene5 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 9 | 9 | 9 | 9 | 9 | 9 | | Gene6 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 7 | 7 | 7 | 7 | 7 | 7 |
mpg_sampleRegionHiearchy is a sample region hiearchy. It is a nested
named list.
data(mpg_sampleRegionHiearchy)
mpg_sampleRegionHiearchy
## $All
## $All$`Region 1`
## [1] ""
##
## $All$`Region 2`
## [1] ""
nametree(mpg_sampleRegionHiearchy)
## All
## ├──Region 1
## └──Region 2
In this example Region 1 and Region 2 are subsets of All region
Selection of marker genes
Marker gene selection is performed using three functions:
markerCandidates, pickMarkers and rotateSelect. By default
markerCandidates will return files for each cell type in a region that
lists the gene that are above a given silhouette and fold change
thresholds. Other variables are briefly explained below but see the
package documentation for in depth explanations
markerCandidates(design = mgp_sampleProfilesMeta, # the design file
expression = mgp_sampleProfiles, # expression file
outLoc = 'README_files/quickSelection', # output directory
groupNames = 'CellType', # name of the column with cell types. can be a vector
regionNames = 'region', # name of the column with brain regions. leave NULL if no region seperation is desired
PMID = 'PMID', # name of the column with study identifiers
sampleName = 'sampleName', # name of the column with sample names
replicates = 'replicate', # name of the column with replicates
foldChangeThresh = 10, # threshold of fold change for gene selection (default is 10)
minimumExpression = 8, # minimum expression level that a gene can be considered a marker gene (default is 8)
background = 6, # background level of expression (default is 6)
regionHierarchy = mpg_sampleRegionHiearchy, # hierarchy of brain regions to be used
geneID = 'Gene.Symbol', # column name with with gene idenditifers
cores = 8 # number of cores to use in parallelization
)
This creates 3 directories in the output directory
list.files('README_files/quickSelection')
## [1] "All_CellType" "CellType" "Region 2_CellType"
The CellType directory is a list of marker genes that disregards all
region specifications (redundant with All_CellType in this case) while
Region 2_CellType and All_CellType directories inlcude cell types
from the relevant region. Note the absence of Region 1_CellType since
that region only has a single cell
type.
read.table('README_files/quickSelection/All_CellType/Cell C') %>% knitr::kable()
| V1 | V2 | V3 | V4 | | :---- | -: | -: | -: | | Gene3 | 15 | 1 | 1 | | Gene4 | 12 | 1 | 1 | | Gene5 | 8 | 1 | 1 |
This file shows the candidate genes for cell type Cell C in region
All. The first column is the gene identifier, the second is change in
expression in log_2 scale and the third one is the silhouette
coefficient. Note that Gene6 is absent since its expression level was
below the minimum level allowed. markerCandidates function does not
apply a threshold for silhouette coefficient it also doesn’t check to
see if a gene satisfies fold change threshold for multiple genes.
pickMarkers function does that.
pickMarkers('README_files/quickSelection/All_CellType/',
foldChange = 1, # this is a fixed fold change threshold that ignores some leniency that comes from markerCandidates. setting it to 1 makes it irrelevant
silhouette = 0.5)
## $`Cell A`
## [1] "Gene1"
##
## $`Cell B`
## [1] "Gene2"
##
## $`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
If all genes for all regions needs to be seen
pickMarkersAll('README_files/quickSelection',
foldChange = 1,
silhouette = 0.5)
## $All_CellType
## $All_CellType$`Cell A`
## [1] "Gene1"
##
## $All_CellType$`Cell B`
## [1] "Gene2"
##
## $All_CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $CellType
## $CellType$`Cell A`
## [1] "Gene1"
##
## $CellType$`Cell B`
## [1] "Gene2"
##
## $CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $`Region 2_CellType`
## $`Region 2_CellType`$`Cell B`
## [1] "Gene2"
##
## $`Region 2_CellType`$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
Better selection of marker genes
The above method is a quick way to pick markers but it does not handle
bimodality in expression distribution well. To ensure robustness of the
results it is better to perform multiple selections with permutations.
markerCandidates function has variables to handle permutations for
you. rotate controls what is the percentage of samples that should be
removed every time. seed controls the random seed and is there to ensure
reproducibility.
for (i in 1:10){
markerCandidates(design = mgp_sampleProfilesMeta, # the design file
expression = mgp_sampleProfiles, # expression file
outLoc = file.path('README_files/Rotation',i), # output directory
groupNames = 'CellType', # name of the column with cell types. can be a vector
regionNames = 'region', # name of the column with brain regions. leave NULL if no region seperation is desired
PMID = 'PMID', # name of the column with study identifiers
sampleName = 'sampleName', # name of the column with sample names
replicates = 'replicate', # name of the column with replicates
foldChangeThresh = 10, # threshold of fold change for gene selection (default is 10)
minimumExpression = 8, # minimum expression level that a gene can be considered a marker gene (default is 8)
background = 6, # background level of expression (default is 6)
regionHierarchy = mpg_sampleRegionHiearchy, # hierarchy of brain regions to be used
geneID = 'Gene.Symbol', # column name with with gene idenditifers
cores = 8, # number of cores to use in parallelization
rotate = 0.33,
seed = i
)
}
This creates multiple selection directories. rotateSelect can be used
to count the number of times a gene is selected for each cell type in
each region. This creates another directory similar to the output of
markerCandidates. Again, valid markers can be acquired using
pickMarkers
rotateSelect(rotationOut='README_files/Rotation',
rotSelOut='README_files/RotSel',
cores = 8,
foldChange = 1 # this is a fixed fold change threshold that ignores some leniency that comes from markerCandidates. setting it to 1 makes it irrelevant
)
pickMarkers('README_files/RotSel/All_CellType/',rotationThresh = 0.95)
## $`Cell A`
## [1] "Gene1"
##
## $`Cell B`
## [1] "Gene2"
##
## $`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
pickMarkersAll('README_files/RotSel',rotationThresh = 0.95)
## $All_CellType
## $All_CellType$`Cell A`
## [1] "Gene1"
##
## $All_CellType$`Cell B`
## [1] "Gene2"
##
## $All_CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $CellType
## $CellType$`Cell A`
## [1] "Gene1"
##
## $CellType$`Cell B`
## [1] "Gene2"
##
## $CellType$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
##
##
## $`Region 2_CellType`
## $`Region 2_CellType`$`Cell B`
## [1] "Gene2"
##
## $`Region 2_CellType`$`Cell C`
## [1] "Gene3" "Gene4" "Gene5"
Marker gene profiles (MGP)
###Sample data for MGP estimation
The package includes mouse brain cell type markers published in Mancarci et al. 2017 and gene expression data from substantia nigra samples from healthy donors and Parkinson’s disease patients by Lesnick et al. 2007.
Mouse marker genes is available in mouseMarkerGenes object as a nested
list.
data(mouseMarkerGenes)
names(mouseMarkerGenes)
## [1] "All" "Amygdala" "BasalForebrain"
## [4] "Brainstem" "Cerebellum" "Cerebrum"
## [7] "Cortex" "Hippocampus" "LocusCoeruleus"
## [10] "Midbrain" "SpinalCord" "Striatum"
## [13] "Subependymal" "SubstantiaNigra" "Thalamus"
lapply(mouseMarkerGenes$Midbrain[1:3],head, 14)
## $Astrocyte
## [1] "Aass" "Acsbg1" "Acsl6" "Acss1" "Add3" "Adhfe1"
## [7] "AI464131" "Aldh1l1" "Aldh6a1" "Antxr1" "Aox1" "Apoe"
## [13] "Aqp4" "Axl"
##
## $BrainstemCholin
## [1] "2310030G06Rik" "Anxa2" "Cabp1" "Calca"
## [5] "Calcb" "Cd24a" "Cd55" "Cda"
## [9] "Chodl" "Ecel1" "Fxyd7" "Hebp2"
## [13] "Hspb1" "Hspb8"
##
## $Dopaminergic
## [1] "Cacna2d2" "Cadps2" "Chrna6" "Mapk8ip2" "Nr4a2" "Ntn1"
## [7] "Prkcg" "Rian" "Scn2a" "Slc6a3" "Snhg11" "Tenm1"
## [13] "Th" "Zim3"
Available Lesnick et al. data is stored in mgp_LesnickParkinsonsExp
and mgp_LesnickParkinsonsMeta objects
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:testthat':
##
## matches
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
data(mgp_LesnickParkinsonsExp)
mgp_LesnickParkinsonsExp %>%
dplyr::select(-GeneNames) %>%
head %>% {.[,1:6]}
## Probe Gene.Symbol NCBIids GSM184354.cel GSM184355.cel
## 3366 1007_s_at DDR1 780 10.236880 9.891552
## 44099 1053_at RFC2 5982 5.421790 5.280541
## 14880 117_at HSPA7|HSPA6 3311|3310 5.164445 4.651754
## 11594 121_at PAX8 7849 7.076004 7.035090
## 11907 1255_g_at GUCA1A 2978 3.107388 3.418976
## 2069 1294_at UBA7 7318 6.644858 6.182664
## GSM184356.cel
## 3366 10.498371
## 44099 5.852467
## 14880 4.729189
## 11594 6.698765
## 11907 3.491832
## 2069 5.982642
data(mgp_LesnickParkinsonsMeta)
mgp_LesnickParkinsonsMeta %>% head
## GSM disease
## 1 GSM184354 Control
## 2 GSM184355 Control
## 3 GSM184356 Control
## 4 GSM184357 Control
## 5 GSM184358 Control
## 6 GSM184359 Control
Before MGP estimation, it is important to filter expression data of low expressed genes and make sure all genes are represented only once in the dataset. While there are many probeset summarization and methods, for this work we chose the most variable probeset and remove all probes with a median expression below the median expression of the dataset.
unfilteredParkinsonsExp = mgp_LesnickParkinsonsExp # keep this for later
medExp = mgp_LesnickParkinsonsExp %>%
sepExpr() %>% {.[[2]]} %>%
unlist %>% median
# mostVariable function is part of this package that does probe selection and filtering for you
mgp_LesnickParkinsonsExp = mostVariable(mgp_LesnickParkinsonsExp,
threshold = medExp,
threshFun= median)
MGP estimation
MGPs are simply the first principal component of marker gene expression. This method of marker gene profile estimation is similar to the methodology of multiple previous works that aim to estimate relative abundance of cell types in a whole tissue sample (Chikina et al., 2015; Westra et al., 2015; Xu et al., 2013).
Primary function that deals with MGP estimation is mgpEstimate. This
function will take in expression data and marker gene lists to return
MGPs. exprData is the expression matrix which should include gene
identifiers as a column. Other columns can be present at the beginning
of the data frame but should not be of class double. genes is the
list of markers. It is assumed to be a list of character vectors, each
vector containing gene names for a cell type. geneColName is the name
of the column where gene names are found. geneTransform is a function
that will be applied to the gene names provided in genes. Note that
this by default tranforms mouse gene names to human gene names. Set it
to NULL if this is not desired.
Below a basic estimation performed with other important variables briefly explained.
estimations = mgpEstimate(exprData=mgp_LesnickParkinsonsExp,
genes=mouseMarkerGenes$Midbrain,
geneColName='Gene.Symbol',
outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
geneTransform =function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
groups=mgp_LesnickParkinsonsMeta$disease, #if there are experimental groups provide them here. if not desired set to NULL
seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
removeMinority = TRUE) # removes genes if they are the minority in terms of rotation sign from estimation process
Dopamine rgic cell loss is a known effect of Parkinson’s Disease. To see if this effect can be observed we can look at dopaminergic MGPs in healthy donors vs Parkinson’s Disease patients
library(ggplot2)
ls(estimations)
## [1] "allMarkerExpression" "estimates"
## [3] "fullPCAs" "groups"
## [5] "meanUsedMarkerExpression" "removedMarkerRatios"
## [7] "rotations" "simpleScaledEstimation"
## [9] "trimmedPCAs" "usedMarkerExpression"
ls(estimations$estimates)
## [1] "Astrocyte" "BrainstemCholin"
## [3] "Dopaminergic" "Microglia"
## [5] "Microglia_activation" "Microglia_deactivation"
## [7] "Oligo" "Serotonergic"
dopaminergicFrame =
data.frame(`Dopaminergic MGP` = estimations$estimates$Dopaminergic,
state = estimations$groups$Dopaminergic, # note that unless outlierSampleRemove is TRUE this will be always the same as the groups input
check.names=FALSE)
ggplot2::ggplot(dopaminergicFrame,
aes(x = state, y = `Dopaminergic MGP`)) +
geom_boxvio() + geom_jitter(width = .05) # this is just a convenience function that outputs a list of ggplot elements.
To see if the difference between the two groups is statistically significant we use wilcoxon test (Mann-Whitney test). Wilcoxon-test is a non parametric tests that does not make assumptions about the distribution of the data. We chose to use it because marker gene profiles are not normally distributed.
group1 = estimations$estimates$Dopaminergic[estimations$groups$Dopaminergic %in% "Control"]
group2 = estimations$estimates$Dopaminergic[estimations$groups$Dopaminergic %in% "PD"]
wilcox.test(group1,group2)
##
## Wilcoxon rank sum test
##
## data: group1 and group2
## W = 119, p-value = 0.006547
## alternative hypothesis: true location shift is not equal to 0
Based on these results we can say that there is indeed a significant difference between doparminergic marker gene profiles between control and parkinson’s disease patients.
This difference can be also observed by looking at gene expression of markers
estimations$usedMarkerExpression$Dopaminergic%>%
as.matrix %>%
gplots::heatmap.2(trace = 'none',
scale='row',Rowv = FALSE,Colv = FALSE, dendrogram = 'none',
col= viridis::viridis(10),cexRow = 1.5, cexCol = 1,
ColSideColors = estimations$groups$Dopaminergic %>%
toColor(palette = c('Control' = 'blue',
"PD" = "red")) %$% cols ,
margins = c(7,8))
Indeed in general most marker genes have a high expression is control samples which are shown in green.
It is important to note that not all marker genes are used in the estimation of marker gene profiles. Below we see all genes that are identified as dopaminergic cell markers with human orthologues. Yet not all these genes appear in the heatmap above
mouseHumanGeneTable = mouseMarkerGenes$Midbrain$Dopaminergic %>% homologene::mouse2human()
allHumanDopaGenes = mouseHumanGeneTable %$% humanGene
mouseHumanGeneTable
## mouseGene humanGene mouseID humanID
## 1 Cacna2d2 CACNA2D2 56808 9254
## 2 Cadps2 CADPS2 320405 93664
## 3 Chrna6 CHRNA6 11440 8973
## 4 Mapk8ip2 MAPK8IP2 60597 23542
## 5 Nr4a2 NR4A2 18227 4929
## 6 Ntn1 NTN1 18208 9423
## 7 Prkcg PRKCG 18752 5582
## 8 Slc6a3 SLC6A3 13162 6531
## 9 Tenm1 TENM1 23963 10178
## 10 Th TH 21823 7054
Some of these genes are removed because they are not included in our dataset, either because they are not in the platform, or because the gene was filtered in the pre-processing stage due to low expression.
allHumanDopaGenes[!allHumanDopaGenes %in% mgp_LesnickParkinsonsExp$Gene.Symbol]
## [1] "CHRNA6"
allGenes = allHumanDopaGenes[allHumanDopaGenes %in% mgp_LesnickParkinsonsExp$Gene.Symbol]
Other genes can be removed because mgpEstimate thinks they don’t correlate well with the rest of the genes. Details of this process can be found in Mancarci et al. 2017 manuscript
allGenes[!allGenes %in% rownames(estimations$usedMarkerExpression$Dopaminergic)]
## [1] "PRKCG"
genesUsed = rownames(estimations$usedMarkerExpression$Dopaminergic)
By looking at the unfiltered expression data, we can see how these two genes were unsuitable for MGP estimation
toPlot =
unfilteredParkinsonsExp[unfilteredParkinsonsExp$Gene.Symbol %in%
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic)$humanGene,] %>%
mostVariable(threshold = 0)
toPlot %>%
mostVariable(threshold = 0) %>%
sepExpr() %>%
{.[[2]]} %>% as.matrix() %>%
apply(1,function(x){scale01(x)}) %>%
t %>%
{rownames(.) =
toPlot$Gene.Symbol[toPlot$Gene.Symbol %in%
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic)$humanGene];.} %>%
reshape2::melt() %>%
{colnames(.) = c('Gene','Sample','Expression');.} %>%
dplyr::mutate(`Is used?` = rep('used',length(Gene)) %>%
{.[Gene %in% 'CHRNA6'] = '(CHRNA6) not expressed';.[Gene %in% 'PRKCG'] = '(PRKCG) not correlated';.}) %>%
ggplot(aes(y = Expression, x = Sample, group = Gene, color = `Is used?`)) +
geom_line() +
cowplot::theme_cowplot() +
theme( axis.text.x= element_blank(),
axis.title = element_text(size = 20),
legend.text = element_text(size = 17),
legend.title = element_text(size=22),
plot.title = element_text(size = 20)
) +
ggtitle('Scaled expression of markers')
toPlot %>%
mostVariable(threshold = 0) %>%
sepExpr() %>%
{.[[2]]} %>% as.matrix()%>%
{rownames(.) =
toPlot$Gene.Symbol[toPlot$Gene.Symbol %in%
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic)$humanGene];.} %>%
reshape2::melt() %>%
{colnames(.) = c('Gene','Sample','Expression');.} %>%
dplyr::mutate(`Is used?` = rep('used',length(Gene)) %>%
{.[Gene %in% 'CHRNA6'] = 'CHRNA6 - not expressed';.[Gene %in% 'PRKCG'] = 'PRKCG - not correlated';.}) %>%
ggplot(aes(y = Expression, x = Sample, group = Gene, color = `Is used?`)) +
geom_line() +
cowplot::theme_cowplot() +
theme( axis.text.x= element_blank(),
axis.title = element_text(size = 20),
legend.text = element_text(size = 17),
legend.title = element_text(size=22),
plot.title = element_text(size = 20)
)+
ggtitle('Nonscaled expression of markers')
Both of these genes seem to be non-expressed though PRKCG managed to be just above the removal threshold. Luckily lack of correlation reveals that it is not behaving as other marker genes.
The ratio of genesUsed and allGenes (all markers available in the
study) can be used as a confidence metric. If a significant portion of
the genes do not correlate well with each other that may point to
presence of non cell type specific signal (regulation or noise). For all
cell types this ratio is outputted. You’ll see a warning if this ratio
ever exceeds
0.4
estimations$removedMarkerRatios
## Astrocyte BrainstemCholin Dopaminergic
## 0.12264151 0.35294118 0.11111111
## Microglia Microglia_activation Microglia_deactivation
## 0.09734513 0.08333333 0.12149533
## Oligo Serotonergic
## 0.05882353 0.00000000
The ratio for dopaminergic cells is fairly low. We can also look at the amount of varience explained in the first PC of marker gene expression. If this value is low, it can point to higher amounts of confounding noise or regulation.
estimations$trimmedPCAs$Dopaminergic %>% summary()
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.3056 1.0829 0.87149 0.53287 0.45796 0.38075
## Proportion of Variance 0.6645 0.1466 0.09494 0.03549 0.02622 0.01812
## Cumulative Proportion 0.6645 0.8110 0.90597 0.94147 0.96768 0.98580
## PC7 PC8
## Standard deviation 0.28376 0.18177
## Proportion of Variance 0.01006 0.00413
## Cumulative Proportion 0.99587 1.00000
For instance unlike dopaminergic cells, cholinergic cells seem to have a higher proportion of their marker genes removed. Looking at the variation explained by first PC reveals that first PC only explains 28% of all variance.
estimations$trimmedPCAs$BrainstemCholin %>% summary()
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7194 1.4430 1.2543 1.1992 0.8893 0.82035 0.77928
## Proportion of Variance 0.2688 0.1893 0.1430 0.1307 0.0719 0.06118 0.05521
## Cumulative Proportion 0.2688 0.4581 0.6011 0.7318 0.8037 0.86491 0.92012
## PC8 PC9 PC10 PC11
## Standard deviation 0.59534 0.48891 0.44505 0.29527
## Proportion of Variance 0.03222 0.02173 0.01801 0.00793
## Cumulative Proportion 0.95234 0.97407 0.99207 1.00000
Such a result can occur if
- Marker genes are highly regulated.
- There is little difference in cell type proportions or wholescale regulation of marker genes and MGP is driven by noise and other confounds.
Note about homologene
In this tutorial homologene database is used through
homologene
package. Note that the default version of this database is somewhat old
which means some annotations may have been changed. The homologene
package also includes an updated version of the homologene package to
make sure you get all the correct matches.
If we look at dopaminergic cell markers for instance we could have gotten an extra gene (SCN2A) if we used the updated version.
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic) %$% humanGene
## [1] "CACNA2D2" "CADPS2" "CHRNA6" "MAPK8IP2" "NR4A2" "NTN1"
## [7] "PRKCG" "SLC6A3" "TENM1" "TH"
homologene::mouse2human(mouseMarkerGenes$Midbrain$Dopaminergic,
db = homologene::homologeneData2) %$% humanGene
## [1] "CACNA2D2" "CADPS2" "CHRNA6" "MAPK8IP2" "NR4A2" "NTN1"
## [7] "PRKCG" "SCN2A" "SLC6A3" "TENM1" "TH"
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.