Nothing
bmass Introduction -- Real Data
In bmass: Bayesian Multivariate Analysis of Summary Statistics
This is an advanced introductory vignette for running bmass. Here, we will use real data and explore multiple components of the bmass() output. Specifically, we will:
- Download our real dataset
- Quality control (QC) & format our dataset
- Run
bmass
- Access & explore a variety of results
For a more basic introductory vignette using a smaller, simulated dataset, please see here.
Downloading GlobalLipids2013
For this vigenette, we will be download and analyze the GlobalLipids2013 dataset. To download the necessary files from http://csg.sph.umich.edu/willer/public/lipids2013/, run the following:
mkdir my-directory/GlobalLipids2013
cd my-directory/GlobalLipids2013
wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_LDL.txt.gz
wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_HDL.txt.gz
wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_TG.txt.gz
wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_TC.txt.gz
Formatting & QC'ing our dataset
To run bmass, it is recommended that at least a minimal level of QC is performed on the datasets of interest. Additionally, it is important that all the single nucleotide polymorphisms (SNPs) being analyzed are oriented to the same reference allele across all the phenotypes being included. This latter situation is emphasized so that the correct direction of effect is assigned to each phenotype for a given SNP (which is oriented based on reference allele).
Below is example bash code on how to reformat, conduct some basic QC, and orient the reference alleles on the GlobalLipids2013 dataset. First, we will reformat our input datafiles to match the requirements for bmass. bmass expects input files to have the following columns --
- Chr: Chromosome
- BP: Basepair Position
- Marker: Variant ID information, such as a rsID # (should be a unique value for each variant)
- MAF: Minor Allele Frequency
- A1: Reference allele
- A2: AlternativeaAllele
- Direction: Direction of effect size (a column of "+" or "-")
- pValue: p-Value from univariate GWAS
- N: Sample size
To create files that match these specifications, we will conduct the following: extract the columns that are needed, orient SNPs to the minor allele, and extract the direction of effect from the provided effect sizes.
zcat jointGwasMc_HDL.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_HDL.formatted.txt.gz
zcat jointGwasMc_LDL.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_LDL.formatted.txt.gz
zcat jointGwasMc_TG.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_TG.formatted.txt.gz
zcat jointGwasMc_TC.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_TC.formatted.txt.gz
Next, we will conduct some basic QC. These steps include: remove SNPs that have duplicate entries, that do not have entries in every field, do not have MAF information, that are fixed (MAF = 0), are rare (MAF < .01), have sample sample sizes < 50000, and have effect sizes of 0 (ie direction of effect is indeterminable). We will also orient each SNP/phenotype file to its respective minor allele.
join <(zcat jointGwasMc_HDL.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_HDL.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_HDL.formatted.QCed.txt.gz
join <(zcat jointGwasMc_LDL.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_LDL.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | gzip > jointGwasMc_LDL.formatted.QCed.txt.gz
join <(zcat jointGwasMc_TG.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_TG.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | gzip > jointGwasMc_TG.formatted.QCed.txt.gz
join <(zcat jointGwasMc_TC.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_TC.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | gzip > jointGwasMc_TC.formatted.QCed.txt.gz
Lastly, we will orient all phenotype files to the HDL minor allele:
join <(zcat jointGwasMc_HDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $5 }' | sort -k 1,1) <(zcat jointGwasMc_LDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'if ($F[1] eq $F[7]) { ($F[6], $F[7]) = ($F[7], $F[6]); if ($F[8] eq "-") { $F[8] = "+"; } elsif ($F[8] eq "+") { $F[8] = "-"; } else { print STDERR "Error1 -- direction of effect allele matching"; } } print join("\t", @F[2..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_LDL.formatted.QCed.HDLMatch.txt.gz
join <(zcat jointGwasMc_HDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $5 }' | sort -k 1,1) <(zcat jointGwasMc_TG.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'if ($F[1] eq $F[7]) { ($F[6], $F[7]) = ($F[7], $F[6]); if ($F[8] eq "-") { $F[8] = "+"; } elsif ($F[8] eq "+") { $F[8] = "-"; } else { print STDERR "Error1 -- direction of effect allele matching"; } } print join("\t", @F[2..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_TG.formatted.QCed.HDLMatch.txt.gz
join <(zcat jointGwasMc_HDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $5 }' | sort -k 1,1) <(zcat jointGwasMc_TC.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'if ($F[1] eq $F[7]) { ($F[6], $F[7]) = ($F[7], $F[6]); if ($F[8] eq "-") { $F[8] = "+"; } elsif ($F[8] eq "+") { $F[8] = "-"; } else { print STDERR "Error1 -- direction of effect allele matching"; } } print join("\t", @F[2..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_TC.formatted.QCed.HDLMatch.txt.gz
Running bmass
Now with our four GlobalLipids2013 files properly formatted and QCed, we are able to run bmass. To do so, we will also need the list of published univariate GWAS SNPs from the GlobalLipids2013 publication (presented as two columns per SNP, chromosome and basepair position) -- these have been provided in the data subdirectory (GlobalLipids2013.GWASsnps.txt). To run bmass with our prepared input datafiles and our list of univariate GWAS SNPs, we run the following R code:
library("bmass");
HDL <- read.table("jointGwasMc_HDL.formatted.QCed.txt.gz", header=T);
LDL <- read.table("jointGwasMc_LDL.formatted.QCed.HDLMatch.txt.gz", header=T);
TG <- read.table("jointGwasMc_TG.formatted.QCed.HDLMatch.txt.gz", header=T);
TC <- read.table("jointGwasMc_TC.formatted.QCed.HDLMatch.txt.gz", header=T);
load("bmassDirectory/data/GlobalLipids2013.GWASsnps.rda");
Phenotypes <- c("HDL", "LDL", "TG", "TC");
bmassResults <- bmass(Phenotypes, GlobalLipids2013.GWASsnps);
Note once again that to run bmass we supply a vector of phenotype names that also correspond to the variable names referencing each phenotype's respective data file. Also note that you will see a warning that one or more of the GlobalLipids2013 datafiles contains p-values smaller than the threshold at which R can handle converting them into logarthmic scale -- this is to be expected.
Accessing & exploring results
bmass returns a list containing multiple areas of information, including separate sublists pertaining to the previous univariate GWAS SNPs used for training the multivariate model priors and the new multivariate SNPs (if any) bmass has idenfitifed as significant. Run summary(bmassResults) and you should see the following:
> summary(bmassResults)
Length Class Mode
MarginalSNPs 3 -none- list
PreviousSNPs 4 -none- list
NewSNPs 3 -none- list
LogFile 20 -none- character
ZScoresCorMatrix 16 -none- numeric
Models 324 -none- numeric
ModelPriors 1134 -none- numeric
GWASlogBFMinThreshold 1 -none- numeric
We will touch on all the above outputs.
NewSNPs
The NewSNPs sublist contains information pertaining to the new significant multivariate associations found by bmass, if any were identified. Running summary(bmassResults$NewSNPs) you should see the following:
> summary(bmassResults$NewSNPs)
Length Class Mode
SNPs 30 data.frame list
logBFs 5427 -none- numeric
Posteriors 5427 -none- numeric
bmassResults$NewSNPs$SNPs contains the 'main' bmass output of interest -- the list of new multivariate associations found by bmass and those SNPs' related information. The format of this output appears as follows:
> head(bmassResults$NewSNPs$SNPs, n=3)
ChrBP Chr BP Marker MAF A1 HDL_A2 HDL_Direction
1704 10_101902054 10 101902054 rs2862954 0.4631 C T +
72106 10_5839619 10 5839619 rs2275774 0.1781 G A +
118903 11_109521729 11 109521729 rs661171 0.2876 T G +
HDL_pValue HDL_N HDL_ZScore LDL_Direction LDL_pValue LDL_N
1704 1.287e-06 186893 4.841751 + 5.875e-01 172821.0
72106 7.601e-07 179144 4.945343 - 7.773e-05 165198.0
118903 1.705e-06 186946 4.785573 + 1.653e-02 172877.9
LDL_ZScore TG_Direction TG_pValue TG_N TG_ZScore TC_Direction
1704 0.5424624 + 0.013930 177587.1 2.459063 +
72106 -3.9512933 - 0.001035 169853.0 -3.280836 -
118903 2.3969983 - 0.155800 177645.0 -1.419340 +
TC_pValue TC_N TC_ZScore GWASannot mvstat mvstat_log10pVal unistat
1704 2.526e-04 187083 3.659609 0 48.70946 9.173067 23.44255
72106 1.911e-02 179333 -2.343378 0 38.26288 7.004804 24.45642
118903 1.785e-05 187131 4.290215 0 37.84098 6.917765 22.90171
unistat_log10pVal Nmin logBFWeightedAvg
1704 5.890421 172821.0 7.068306
72106 6.119129 165198.0 5.447201
118903 5.768276 172877.9 5.438490
NewSNPs$logBFs contains the log10 Bayes Factor (BF) for each multivariate model tested (in the form of a Model x SNP matrix, with the first columns of the matrix describing each model). NewSNPs$Posteriors is a similar matrix, just with each matrix entry representing the posterior probability for each model/SNP combination instead of the log10 BF.
> dim(bmassResults$NewSNPs$logBFs)
[1] 81 67
> bmassResults$NewSNPs$logBFs[1:5,1:10]
HDL LDL TG TC 10_101902054 10_5839619 11_109521729 11_13313759
[1,] 0 0 0 0 0.000000 0.0000000 0.0000000 0.00000000
[2,] 1 0 0 0 -233.831047 -235.0781367 -234.6670299 -234.15472067
[3,] 2 0 0 0 0.000000 0.0000000 0.0000000 0.00000000
[4,] 0 1 0 0 0.165855 0.3172489 -0.1100418 0.06393596
[5,] 1 1 0 0 -64.774919 -66.2959645 -69.2388829 -69.93309528
11_45696596 11_47251202
[1,] 0.00000000 0.00000000
[2,] -231.68478695 -219.10321932
[3,] 0.00000000 0.00000000
[4,] -0.04838241 0.04997886
[5,] -65.59917325 -45.04269241
PreviousSNPs
PreviousSNPs contains similar information as NewSNPs does, but pertaining to the set of previous univariate GWAS significant SNPs used to train the multivariate model priors and set the empirical significance threshold for weighted average BF. Running summary(bmassResults$PreviousSNPs) you see:
> summary(bmassResults$PreviousSNPs)
Length Class Mode
logBFs 12069 -none- numeric
SNPs 30 data.frame list
DontPassSNPs 30 data.frame list
Posteriors 12069 -none- numeric
PreviousSNPs$SNPs, PreviousSNPs$logBFs, and PreviousSNPs$Posteriors all match the descriptions as above with NewSNPs. However PreviousSNPs$DontPassSNPs refers to any GWAS SNPs which are included in the final publicaiton list but which do meet the original study's univariate GWAS p-value threshold based on the dataset released. This can for instance occur when a final GWAS SNP list is produced from combining a discovery dataset with a replication dataset, and the former data is the only part of the study publicly released. Therefore bmass keeps track of these SNPs so as not to call them a completely novel multivariate SNP, but also to not use them to train the multivariate model priors since they do not technically reach univariate GWAS significance based on the dataset provided alone.
The format of DontPassSNPs is the same seen in NewSNPs$SNPs and PreviousSNPs$SNPs.
MarginalSNPs
MarginalSNPs contains the same information as NewSNPs and PreviousSNPs, but pertaining to the set of marginally significant SNPs extracted from the intermediate MergedDataSources dataset (a data.frame that contains all the input phenotype files combined) based on the SNPMarginalUnivariateThreshold and SNPMarginalMultivariateThreshold values (default of 1e-6 for each). summary(bmassResults$MarginalSNPs) returns:
> summary(bmassResults$MarginalSNPs)
Length Class Mode
SNPs 30 data.frame list
logBFs 20493 -none- numeric
Posteriors 20493 -none- numeric
MarginalSNPs$SNPs, MarginalSNPs$logBFs, and MarginalSNPs$Posteriors are all as described above.
ZScoresCorMatrix
This is the phenotype correlation matrix (V_0 hat in Detailed Methods (Global Lipids Analysis) of Stephens 2013 PLoS ONE) derived from extracting all the 'null' SNPs (abs(ZScore) < 2 for every phenotype) in the dataset.
> bmassResults$ZScoresCorMatrix
HDL_ZScore LDL_ZScore TG_ZScore TC_ZScore
HDL_ZScore 1.0000000 -0.0872789 -0.3655508 0.1523894
LDL_ZScore -0.0872789 1.0000000 0.1607208 0.8223175
TG_ZScore -0.3655508 0.1607208 1.0000000 0.2892982
TC_ZScore 0.1523894 0.8223175 0.2892982 1.0000000
Models
This is a matrix displaying all the model combinations possible given the set of d input phenotypes and 3 model categories of directly associated (1), indirectly associated (2), and unassociated (0). The total number of models displayed should be 3^d.
> dim(bmassResults$Models)
[1] 81 4
> head(bmassResults$Models)
HDL LDL TG TC
[1,] 0 0 0 0
[2,] 1 0 0 0
[3,] 2 0 0 0
[4,] 0 1 0 0
[5,] 1 1 0 0
[6,] 2 1 0 0
ModelPriors
This is the set of trained priors determined for each multivariate model using the previous univariate GWAS SNPs and an EM algorithm (see eq. 36 in Stephens 2013). Note that for any model which does not contain at least one phenotype in the directly associated category, the prior was automatically set to 0 (aside from the global null model of all phenotypes set to unassociated).
Also note that ModelPriors contains a full set of model priors for each value of the hyperparameter sigma_alpha (see "Prior on sigma_alpha" in Stephens 2013); these values are stored in SigmaAlphas, and by default are set to c(0.005,0.0075,0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.15) (hence 81 models * 14 sigma_alphas = 1134 priors). To see a more interpretable form of the model priors, use the GetModelPriorMatrix() function as shown below.
> length(bmassResults$ModelPriors)
[1] 1134
> bmassResults[c("ModelPriorMatrix", "LogFile")] <- GetModelPriorMatrix(Phenotypes, bmassResults$Models, bmassResults$ModelPriors, bmassResults$LogFile)
> head(bmassResults$ModelPriorMatrix)
HDL LDL TG TC Prior Cumm_Prior OrigOrder
1 1 2 1 1 0.32744537 0.3274454 44
2 1 2 2 1 0.13788501 0.4653304 53
3 1 1 1 2 0.11727440 0.5826048 68
4 1 1 1 1 0.07801825 0.6606230 41
5 1 2 1 2 0.06210658 0.7227296 71
6 2 1 2 2 0.05698876 0.7797184 78
ModelPriorMatrix contains the model descriptions as provided in Models along with the combined trained priors across each sigma_alpha value; the ModelPriomatrix is sorted by decreasing model prior value and also shows the cumulative distribution of the priors as well as the original order position of each model.
GWASlogBFMinThreshold
This is the empirical weighted average BF threshold (BFavg) used to determine whether any new SNPs are 'multivariate significant'. As described in Stephens 2013, each SNP has a single summary metric combining all the information across each multivariate model tested; this summary metric is the weighted average of each model's BF, where weights are determined from the previous univariate GWAS SNPs (eq. 1 in Turchin and Stephens 2019 or eq. 8 in Stephens 2013). This summary metric is also calculated for the previous unviariate GWAS SNPs themselves, and the smallest BFavg among them becomes the threshold above which new SNPs are considered multivariate significant.
> bmassResults$GWASlogBFMinThreshold
[1] 4.289906
So every new multivariate SNP bmass identifies in this GlobalLipids2013 analysis has a BFavg value >= 4.289906 (and is also more than a 1Mb window away from a previous univariate GWAS SNP).
[//]: # ##Additional functions
Try the bmass package in your browser
Any scripts or data that you put into this service are public.
bmass documentation built on May 17, 2019, 9:02 a.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.
This is an advanced introductory vignette for running bmass. Here, we will use real data and explore multiple components of the bmass() output. Specifically, we will:
- Download our real dataset
- Quality control (QC) & format our dataset
- Run
bmass - Access & explore a variety of results
For a more basic introductory vignette using a smaller, simulated dataset, please see here.
Downloading GlobalLipids2013
For this vigenette, we will be download and analyze the GlobalLipids2013 dataset. To download the necessary files from http://csg.sph.umich.edu/willer/public/lipids2013/, run the following:
mkdir my-directory/GlobalLipids2013 cd my-directory/GlobalLipids2013 wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_LDL.txt.gz wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_HDL.txt.gz wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_TG.txt.gz wget http://csg.sph.umich.edu/abecasis/public/lipids2013/jointGwasMc_TC.txt.gz
Formatting & QC'ing our dataset
To run bmass, it is recommended that at least a minimal level of QC is performed on the datasets of interest. Additionally, it is important that all the single nucleotide polymorphisms (SNPs) being analyzed are oriented to the same reference allele across all the phenotypes being included. This latter situation is emphasized so that the correct direction of effect is assigned to each phenotype for a given SNP (which is oriented based on reference allele).
Below is example bash code on how to reformat, conduct some basic QC, and orient the reference alleles on the GlobalLipids2013 dataset. First, we will reformat our input datafiles to match the requirements for bmass. bmass expects input files to have the following columns --
- Chr: Chromosome
- BP: Basepair Position
- Marker: Variant ID information, such as a rsID # (should be a unique value for each variant)
- MAF: Minor Allele Frequency
- A1: Reference allele
- A2: AlternativeaAllele
- Direction: Direction of effect size (a column of "+" or "-")
- pValue: p-Value from univariate GWAS
- N: Sample size
To create files that match these specifications, we will conduct the following: extract the columns that are needed, orient SNPs to the minor allele, and extract the direction of effect from the provided effect sizes.
zcat jointGwasMc_HDL.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_HDL.formatted.txt.gz
zcat jointGwasMc_LDL.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_LDL.formatted.txt.gz
zcat jointGwasMc_TG.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_TG.formatted.txt.gz
zcat jointGwasMc_TC.txt.gz | sed 's/:/ /g' | sed 's/chr//g' | awk '{ print $1 "\t" $2 "\t" $5 "\t" $12 "\t" $6 "\t" $7 "\t" $8 "\t" $11 "\t" $10 }' | perl -lane '$F[4] = uc($F[4]); $F[5] = uc($F[5]); if ($F[6] > 0) { $F[6] = "+"; } elsif ($F[6] < 0) { $F[6] = "-"; } else { $F[6] = 0; } print join("\t", @F);' | grep -v SNP_hg18 | gzip > jointGwasMc_TC.formatted.txt.gz
Next, we will conduct some basic QC. These steps include: remove SNPs that have duplicate entries, that do not have entries in every field, do not have MAF information, that are fixed (MAF = 0), are rare (MAF < .01), have sample sample sizes < 50000, and have effect sizes of 0 (ie direction of effect is indeterminable). We will also orient each SNP/phenotype file to its respective minor allele.
join <(zcat jointGwasMc_HDL.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_HDL.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_HDL.formatted.QCed.txt.gz
join <(zcat jointGwasMc_LDL.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_LDL.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | gzip > jointGwasMc_LDL.formatted.QCed.txt.gz
join <(zcat jointGwasMc_TG.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_TG.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | gzip > jointGwasMc_TG.formatted.QCed.txt.gz
join <(zcat jointGwasMc_TC.formatted.txt.gz | awk '{ print $1 "_" $2 }' | grep -v hg18 | sort | uniq -u) <(zcat jointGwasMc_TC.formatted.txt.gz | grep -v hg | grep -v ^rs | grep -v NA | perl -lane 'if ($F[3] > .5) { ($F[4], $F[5]) = ($F[5], $F[4]); if ($F[6] eq "+") { $F[6] = "-"; } elsif ($F[6] eq "-") { $F[6] = "+"; } else { $F[6] = 0; } $F[3] = 1 - $F[3]; } if ($F[6] ne "0") { if ($F[3] > .01) { if ($F[8] >= 50000) { print join("\t", @F); } } }' | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'print join("\t", @F[1..$#F]);' | gzip > jointGwasMc_TC.formatted.QCed.txt.gz
Lastly, we will orient all phenotype files to the HDL minor allele:
join <(zcat jointGwasMc_HDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $5 }' | sort -k 1,1) <(zcat jointGwasMc_LDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'if ($F[1] eq $F[7]) { ($F[6], $F[7]) = ($F[7], $F[6]); if ($F[8] eq "-") { $F[8] = "+"; } elsif ($F[8] eq "+") { $F[8] = "-"; } else { print STDERR "Error1 -- direction of effect allele matching"; } } print join("\t", @F[2..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_LDL.formatted.QCed.HDLMatch.txt.gz
join <(zcat jointGwasMc_HDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $5 }' | sort -k 1,1) <(zcat jointGwasMc_TG.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'if ($F[1] eq $F[7]) { ($F[6], $F[7]) = ($F[7], $F[6]); if ($F[8] eq "-") { $F[8] = "+"; } elsif ($F[8] eq "+") { $F[8] = "-"; } else { print STDERR "Error1 -- direction of effect allele matching"; } } print join("\t", @F[2..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_TG.formatted.QCed.HDLMatch.txt.gz
join <(zcat jointGwasMc_HDL.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $5 }' | sort -k 1,1) <(zcat jointGwasMc_TC.formatted.QCed.txt.gz | awk '{ print $1 "_" $2 "\t" $0 }' | sort -k 1,1) | perl -lane 'if ($F[1] eq $F[7]) { ($F[6], $F[7]) = ($F[7], $F[6]); if ($F[8] eq "-") { $F[8] = "+"; } elsif ($F[8] eq "+") { $F[8] = "-"; } else { print STDERR "Error1 -- direction of effect allele matching"; } } print join("\t", @F[2..$#F]);' | cat <(echo -e "Chr\tBP\tMarker\tMAF\tA1\tA2\tDirection\tpValue\tN") - | gzip > jointGwasMc_TC.formatted.QCed.HDLMatch.txt.gz
Running bmass
Now with our four GlobalLipids2013 files properly formatted and QCed, we are able to run bmass. To do so, we will also need the list of published univariate GWAS SNPs from the GlobalLipids2013 publication (presented as two columns per SNP, chromosome and basepair position) -- these have been provided in the data subdirectory (GlobalLipids2013.GWASsnps.txt). To run bmass with our prepared input datafiles and our list of univariate GWAS SNPs, we run the following R code:
library("bmass"); HDL <- read.table("jointGwasMc_HDL.formatted.QCed.txt.gz", header=T); LDL <- read.table("jointGwasMc_LDL.formatted.QCed.HDLMatch.txt.gz", header=T); TG <- read.table("jointGwasMc_TG.formatted.QCed.HDLMatch.txt.gz", header=T); TC <- read.table("jointGwasMc_TC.formatted.QCed.HDLMatch.txt.gz", header=T); load("bmassDirectory/data/GlobalLipids2013.GWASsnps.rda"); Phenotypes <- c("HDL", "LDL", "TG", "TC"); bmassResults <- bmass(Phenotypes, GlobalLipids2013.GWASsnps);
Note once again that to run bmass we supply a vector of phenotype names that also correspond to the variable names referencing each phenotype's respective data file. Also note that you will see a warning that one or more of the GlobalLipids2013 datafiles contains p-values smaller than the threshold at which R can handle converting them into logarthmic scale -- this is to be expected.
Accessing & exploring results
bmass returns a list containing multiple areas of information, including separate sublists pertaining to the previous univariate GWAS SNPs used for training the multivariate model priors and the new multivariate SNPs (if any) bmass has idenfitifed as significant. Run summary(bmassResults) and you should see the following:
> summary(bmassResults) Length Class Mode MarginalSNPs 3 -none- list PreviousSNPs 4 -none- list NewSNPs 3 -none- list LogFile 20 -none- character ZScoresCorMatrix 16 -none- numeric Models 324 -none- numeric ModelPriors 1134 -none- numeric GWASlogBFMinThreshold 1 -none- numeric
We will touch on all the above outputs.
NewSNPs
The NewSNPs sublist contains information pertaining to the new significant multivariate associations found by bmass, if any were identified. Running summary(bmassResults$NewSNPs) you should see the following:
> summary(bmassResults$NewSNPs) Length Class Mode SNPs 30 data.frame list logBFs 5427 -none- numeric Posteriors 5427 -none- numeric
bmassResults$NewSNPs$SNPs contains the 'main' bmass output of interest -- the list of new multivariate associations found by bmass and those SNPs' related information. The format of this output appears as follows:
> head(bmassResults$NewSNPs$SNPs, n=3) ChrBP Chr BP Marker MAF A1 HDL_A2 HDL_Direction 1704 10_101902054 10 101902054 rs2862954 0.4631 C T + 72106 10_5839619 10 5839619 rs2275774 0.1781 G A + 118903 11_109521729 11 109521729 rs661171 0.2876 T G + HDL_pValue HDL_N HDL_ZScore LDL_Direction LDL_pValue LDL_N 1704 1.287e-06 186893 4.841751 + 5.875e-01 172821.0 72106 7.601e-07 179144 4.945343 - 7.773e-05 165198.0 118903 1.705e-06 186946 4.785573 + 1.653e-02 172877.9 LDL_ZScore TG_Direction TG_pValue TG_N TG_ZScore TC_Direction 1704 0.5424624 + 0.013930 177587.1 2.459063 + 72106 -3.9512933 - 0.001035 169853.0 -3.280836 - 118903 2.3969983 - 0.155800 177645.0 -1.419340 + TC_pValue TC_N TC_ZScore GWASannot mvstat mvstat_log10pVal unistat 1704 2.526e-04 187083 3.659609 0 48.70946 9.173067 23.44255 72106 1.911e-02 179333 -2.343378 0 38.26288 7.004804 24.45642 118903 1.785e-05 187131 4.290215 0 37.84098 6.917765 22.90171 unistat_log10pVal Nmin logBFWeightedAvg 1704 5.890421 172821.0 7.068306 72106 6.119129 165198.0 5.447201 118903 5.768276 172877.9 5.438490
NewSNPs$logBFs contains the log10 Bayes Factor (BF) for each multivariate model tested (in the form of a Model x SNP matrix, with the first columns of the matrix describing each model). NewSNPs$Posteriors is a similar matrix, just with each matrix entry representing the posterior probability for each model/SNP combination instead of the log10 BF.
> dim(bmassResults$NewSNPs$logBFs) [1] 81 67 > bmassResults$NewSNPs$logBFs[1:5,1:10] HDL LDL TG TC 10_101902054 10_5839619 11_109521729 11_13313759 [1,] 0 0 0 0 0.000000 0.0000000 0.0000000 0.00000000 [2,] 1 0 0 0 -233.831047 -235.0781367 -234.6670299 -234.15472067 [3,] 2 0 0 0 0.000000 0.0000000 0.0000000 0.00000000 [4,] 0 1 0 0 0.165855 0.3172489 -0.1100418 0.06393596 [5,] 1 1 0 0 -64.774919 -66.2959645 -69.2388829 -69.93309528 11_45696596 11_47251202 [1,] 0.00000000 0.00000000 [2,] -231.68478695 -219.10321932 [3,] 0.00000000 0.00000000 [4,] -0.04838241 0.04997886 [5,] -65.59917325 -45.04269241
PreviousSNPs
PreviousSNPs contains similar information as NewSNPs does, but pertaining to the set of previous univariate GWAS significant SNPs used to train the multivariate model priors and set the empirical significance threshold for weighted average BF. Running summary(bmassResults$PreviousSNPs) you see:
> summary(bmassResults$PreviousSNPs) Length Class Mode logBFs 12069 -none- numeric SNPs 30 data.frame list DontPassSNPs 30 data.frame list Posteriors 12069 -none- numeric
PreviousSNPs$SNPs, PreviousSNPs$logBFs, and PreviousSNPs$Posteriors all match the descriptions as above with NewSNPs. However PreviousSNPs$DontPassSNPs refers to any GWAS SNPs which are included in the final publicaiton list but which do meet the original study's univariate GWAS p-value threshold based on the dataset released. This can for instance occur when a final GWAS SNP list is produced from combining a discovery dataset with a replication dataset, and the former data is the only part of the study publicly released. Therefore bmass keeps track of these SNPs so as not to call them a completely novel multivariate SNP, but also to not use them to train the multivariate model priors since they do not technically reach univariate GWAS significance based on the dataset provided alone.
The format of DontPassSNPs is the same seen in NewSNPs$SNPs and PreviousSNPs$SNPs.
MarginalSNPs
MarginalSNPs contains the same information as NewSNPs and PreviousSNPs, but pertaining to the set of marginally significant SNPs extracted from the intermediate MergedDataSources dataset (a data.frame that contains all the input phenotype files combined) based on the SNPMarginalUnivariateThreshold and SNPMarginalMultivariateThreshold values (default of 1e-6 for each). summary(bmassResults$MarginalSNPs) returns:
> summary(bmassResults$MarginalSNPs) Length Class Mode SNPs 30 data.frame list logBFs 20493 -none- numeric Posteriors 20493 -none- numeric
MarginalSNPs$SNPs, MarginalSNPs$logBFs, and MarginalSNPs$Posteriors are all as described above.
ZScoresCorMatrix
This is the phenotype correlation matrix (V_0 hat in Detailed Methods (Global Lipids Analysis) of Stephens 2013 PLoS ONE) derived from extracting all the 'null' SNPs (abs(ZScore) < 2 for every phenotype) in the dataset.
> bmassResults$ZScoresCorMatrix HDL_ZScore LDL_ZScore TG_ZScore TC_ZScore HDL_ZScore 1.0000000 -0.0872789 -0.3655508 0.1523894 LDL_ZScore -0.0872789 1.0000000 0.1607208 0.8223175 TG_ZScore -0.3655508 0.1607208 1.0000000 0.2892982 TC_ZScore 0.1523894 0.8223175 0.2892982 1.0000000
Models
This is a matrix displaying all the model combinations possible given the set of d input phenotypes and 3 model categories of directly associated (1), indirectly associated (2), and unassociated (0). The total number of models displayed should be 3^d.
> dim(bmassResults$Models) [1] 81 4 > head(bmassResults$Models) HDL LDL TG TC [1,] 0 0 0 0 [2,] 1 0 0 0 [3,] 2 0 0 0 [4,] 0 1 0 0 [5,] 1 1 0 0 [6,] 2 1 0 0
ModelPriors
This is the set of trained priors determined for each multivariate model using the previous univariate GWAS SNPs and an EM algorithm (see eq. 36 in Stephens 2013). Note that for any model which does not contain at least one phenotype in the directly associated category, the prior was automatically set to 0 (aside from the global null model of all phenotypes set to unassociated).
Also note that ModelPriors contains a full set of model priors for each value of the hyperparameter sigma_alpha (see "Prior on sigma_alpha" in Stephens 2013); these values are stored in SigmaAlphas, and by default are set to c(0.005,0.0075,0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.15) (hence 81 models * 14 sigma_alphas = 1134 priors). To see a more interpretable form of the model priors, use the GetModelPriorMatrix() function as shown below.
> length(bmassResults$ModelPriors) [1] 1134 > bmassResults[c("ModelPriorMatrix", "LogFile")] <- GetModelPriorMatrix(Phenotypes, bmassResults$Models, bmassResults$ModelPriors, bmassResults$LogFile) > head(bmassResults$ModelPriorMatrix) HDL LDL TG TC Prior Cumm_Prior OrigOrder 1 1 2 1 1 0.32744537 0.3274454 44 2 1 2 2 1 0.13788501 0.4653304 53 3 1 1 1 2 0.11727440 0.5826048 68 4 1 1 1 1 0.07801825 0.6606230 41 5 1 2 1 2 0.06210658 0.7227296 71 6 2 1 2 2 0.05698876 0.7797184 78
ModelPriorMatrix contains the model descriptions as provided in Models along with the combined trained priors across each sigma_alpha value; the ModelPriomatrix is sorted by decreasing model prior value and also shows the cumulative distribution of the priors as well as the original order position of each model.
GWASlogBFMinThreshold
This is the empirical weighted average BF threshold (BFavg) used to determine whether any new SNPs are 'multivariate significant'. As described in Stephens 2013, each SNP has a single summary metric combining all the information across each multivariate model tested; this summary metric is the weighted average of each model's BF, where weights are determined from the previous univariate GWAS SNPs (eq. 1 in Turchin and Stephens 2019 or eq. 8 in Stephens 2013). This summary metric is also calculated for the previous unviariate GWAS SNPs themselves, and the smallest BFavg among them becomes the threshold above which new SNPs are considered multivariate significant.
> bmassResults$GWASlogBFMinThreshold [1] 4.289906
So every new multivariate SNP bmass identifies in this GlobalLipids2013 analysis has a BFavg value >= 4.289906 (and is also more than a 1Mb window away from a previous univariate GWAS SNP).
[//]: # ##Additional functions
Try the bmass package in your browser
Any scripts or data that you put into this service are public.
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.