Nathan Skene & Julien Bryois 2018-11-08
- Introduction
- Installation
- Using the package (basic usage)
- Set parameters to be used for the analysis
- Install and load all the required R packages and data
- Prepare quantile groups for celltypes
- Download summary statistics file & check it is properly formatted
- Map SNPs to Genes
- Run the main cell type association analysis
- Run the conditional cell type association analysis
- Controlling for a second GWAS
- Who do I talk to?
- Citation
This R package contains code used for testing which cell types can explain the heritability signal from GWAS summary statistics. The method was described in our 2018 Nature Genetics paper, "Genetic identification of brain cell types underlying schizophrenia". This package takes GWAS summary statistics & Single Cell Transcriptome specificity data (in EWCE's CTD format) as input. As output it calculates the associations between the GWAS trait and the celltypes.
Before installing this package it is neccesary to install the magma software package. Please download it from https://ctg.cncr.nl/software/magma. The executable should then be copied to /usr/local/bin so that R can find it. Then install this package as follows:
install.packages("devtools")
library(devtools)
install_github("nathanskene/ewce")
install_github("nathanskene/MAGMA_Celltyping")
# Set path the 1000 genomes reference data.
genome_ref_dir = "~/Downloads/g1000_eur"
if(!file.exists(sprintf("%s/g1000_eur.bed",genome_ref_dir))){
download.file("https://ctg.cncr.nl/software/MAGMA/ref_data/g1000_eur.zip",destfile=sprintf("%s.zip",genome_ref_dir))
unzip(sprintf("%s.zip",genome_ref_dir),exdir=genome_ref_dir)
}
genome_ref_path = sprintf("%s/g1000_eur",genome_ref_dir)The EWCE package comes with a celltype specificity dataset which we use as an example. The One2One package is used to obtain 1:1 orthologs.
library(MAGMA.Celltyping)
# Load the celltype data
data(ctd)
# Load the mouse to human 1:1 orthologs
data(ortholog_data_Mouse_Human)Note that the cell type dataset loaded in the code above is the Karolinksa cortex/hippocampus data only. For the full Karolinska dataset with hypothalamus and midbrain instead use the following:
data(ctd_allKI)
Or for the DRONC seq or AIBS datasets use
data(ctd_AIBS)
data(ctd_DRONC_human)
data(ctd_DRONC_human)
First we need to calculate the quantile groups for each celltype within the single cell dataset. This is done using the prepare.quantile.groups function. If your single cell dataset is not from mouse, then change the specificity_species argument. If you wish to use a smaller number of bins then
ctd = prepare.quantile.groups(ctd,specificity_species="mouse",numberOfBins=40)
# Examine how the quantile groups look
print(ctd[[1]]$quantiles[c("Gfap","Dlg4","Aif1"),])
print(table(ctd[[1]]$quantiles[,1]))We need to have a summary statistics file to analyse as input. As an example download summary statistics for Fluid Intelligence, based on the UK Biobank, generated by Ben Neale's group.
The function check.sumstats.formatted.for.magma does some basic processing to get it into the right format.
# Download and unzip the summary statistics file
library(R.utils)
#gwas_sumstats_path = "~/Downloads/20016.assoc.tsv"
gwas_sumstats_path = "/Users/natske/GWAS_Summary_Statistics/20016.assoc.tsv"
if(!file.exists(gwas_sumstats_path)){
download.file("https://www.dropbox.com/s/shsiq0brkax886j/20016.assoc.tsv.gz?raw=1",destfile=sprintf("%s.gz",gwas_sumstats_path))
gunzip(sprintf("%s.gz",gwas_sumstats_path),gwas_sumstats_path)
}
# Format it (i.e. column headers etc)
col_headers = format_sumstats_for_magma(gwas_sumstats_path)genesOutPath = map.snps.to.genes(gwas_sumstats_path,genome_ref_path=genome_ref_path)The analyses can be run in either linear or top10% enrichment modes. Let's start with linear:
ctAssocsLinear = calculate_celltype_associations(ctd,gwas_sumstats_path,genome_ref_path=genome_ref_path,specificity_species = "mouse")
FigsLinear = plot_celltype_associations(ctAssocsLinear,ctd=ctd)Now let's add the top 10% mode
ctAssocsTop = calculate_celltype_associations(ctd,gwas_sumstats_path,genome_ref_path=genome_ref_path,EnrichmentMode = "Top 10%")
FigsTopDecile = plot_celltype_associations(ctAssocsTop,ctd=ctd)Then plot linear together with the top decile mode
ctAssocMerged = merge_magma_results(ctAssocsLinear,ctAssocsTop)
FigsMerged = plot_celltype_associations(ctAssocMerged,ctd=ctd)# Conditional analysis
ctCondAssocs = calculate_conditional_celltype_associations(ctd,gwas_sumstats_path,genome_ref_path=genome_ref_path,analysis_name = "Conditional")
plot_celltype_associations(ctCondAssocs,ctd=ctd)We now want to test enrichments that remain in a GWAS after we control for a second GWAS. So let's download a second GWAS sumstats file and prepare it for analysis.
20018.assoc.tsv is the sumstats file for 'Prospective memory result' from the UK Biobank.
20016.assoc.tsv is the sumstats file for 'Fluid Intelligence Score' from the UK Biobank.
So let's subtract genes associated with prospective memory from those involved in fluid intelligence.
# Download and unzip the summary statistics file
library(R.utils)
gwas_sumstats_path = "~/GWAS_Summary_Statistics/20018.assoc.tsv"
if(!file.exists(gwas_sumstats_path)){
download.file("https://www.dropbox.com/s/shsiq0brkax886j/20016.assoc.tsv.gz?raw=1",destfile=sprintf("%s.gz",gwas_sumstats_path))
gunzip(sprintf("%s.gz",gwas_sumstats_path),gwas_sumstats_path)
}
# Format & map SNPs to genes
col_headers = format_sumstats_for_magma(gwas_sumstats_path)
genesOutPath = map.snps.to.genes(gwas_sumstats_path,genome_ref_path=genome_ref_path)gwas_sumstats_path_Memory = "~/GWAS_Summary_Statistics/20018.assoc.tsv"
gwas_sumstats_path_Intelligence = "~/GWAS_Summary_Statistics/20016.assoc.tsv"
ctAssocsLinearMemory = calculate_celltype_associations(ctd,gwas_sumstats_path_Memory,genome_ref_path=genome_ref_path,specificity_species = "mouse")
ctAssocsLinearIntelligence = calculate_celltype_associations(ctd,gwas_sumstats_path_Intelligence,genome_ref_path=genome_ref_path,specificity_species = "mouse")
plot_celltype_associations(ctAssocsLinearMemory,ctd=ctd)ctAssocMerged_MemInt = merge_magma_results(ctAssocsLinearMemory,ctAssocsLinearIntelligence)
FigsMerged_MemInt = magma.tileplot(ctd=ctd,results=ctAssocMerged_MemInt[[1]]$results,annotLevel=1,fileTag="Merged_MemInt_lvl1",output_path = "~/Desktop")
FigsMerged_MemInt = magma.tileplot(ctd=ctd,results=ctAssocMerged_MemInt[[2]]$results,annotLevel=2,fileTag="Merged_MemInt_lvl2",output_path = "~/Desktop")Check which cell types 'Fluid Intelligence' is associated with after controlling for 'Prospective memory'
gwas_sumstats_path = "/Users/natske/GWAS_Summary_Statistics/20016.assoc.tsv"
memoryGenesOut = "~/GWAS_Summary_Statistics/MAGMA_Files/20018.assoc.tsv.10UP.1.5DOWN/20018.assoc.tsv.10UP.1.5DOWN.genes.out"
ctAssocsLinear = calculate_celltype_associations(ctd,gwas_sumstats_path,genome_ref_path=genome_ref_path,specificity_species = "mouse",genesOutCOND=memoryGenesOut,analysis_name = "ControllingForPropMemory")
FigsLinear = plot_celltype_associations(ctAssocsLinear,ctd=ctd,fileTag = "ControllingForPropMemory")We find that after controlling for prospective memory, there is no significant enrichment left associated with fluid intelligence.
If you have any issues using the package then please get in touch with Nathan Skene (nathan.skene at ki.se). Bug reports etc are all most welcome, we want the package to be easy to use for everyone!
If you use the software then please cite
The package utilises the MAGMA package developed in the Complex Trait Genetics lab at VU university (not us!) so please also cite their work:
de Leeuw, et al. MAGMA: Generalized gene-set analysis of GWAS data. PLoS Comput Biol, 2015.
If you use the EWCE package as well then please cite
If you use the cortex/hippocampus single cell data associated with this package then please cite the following papers:
If you use the midbrain and hypothalamus single cell datasets associated with the 2018 paper then please cite the following papers: