pysumstats 0.4

Package for working with GWAS summary statistics

Patch notes

15-05-2020 (v0.4)

• Fixed export so that saving to a file with '.gz' extension now actaully gzips the output file.
• Fixed data import to enable opening a pysumstats file saved as .pickle by initializing SumStats class. This will also check if the data is still stored when low_ram was set to True.
• Merging pysumstats now computes an overall MAF so meta-analyzed pysumstats can be reimported with pysumstats.
• Added assertions to verify input to all functions a user would use (plot functions, SumStats methods, and cov_matrix_from_phenotype_file)
• Added version number to SumStats objects.
• Added kwargs to SumStats.qc for filtering on custom columns.
• Added custom SumStats warnings. These warnings print always by default you can change this behavior with Pythons default warnings filter.
• Reformatted documentation, now also includes a guide on opening and saving sumstats.
• Moved older patch notes to PATCHNOTES.md to prevent clutter in README

Description

A python package for working with GWAS summary statistics data in Python.
This package is designed to make it easy to read summary statistics, perform QC, merge summary statistics and perform meta-analysis.
Meta-analysis can be performed with .meta() with inverse-variance weighted or samplesize-weighted methods.
GWAMA as described in Baselmans, et al. (2019) can be performed using the .gwama() function in merged summary statistics.
The plotting package uses matplotlib.pyplot for generating figures, so the functions are generally compatible with matplotlib.pyplot colors, and Figure and Axis objects.
Warning: merging with low_memory enabled is still highly experimental.

Reference

Using the pysumstats package for a publication, or something similar? That is awesome!
There is no publication attached to this package, and I am not going to force anyone to reference me or make me a co-author or whatever, I want this to remain easily accessible. But I would greatly appreciate it if you add a link to this github, or a reference to it in the acknowledgements or something like that.
If you have any questions, want to help add methods or want to let me know you are planning a publication with this, you can get in touch via the pypi website of this project.
If you use the .gwama() method, please refer to the original publication: Baselmans, et al. (2019).

Installation

This package was made for Python 3.7. Clone the package directly from this github, or install with

pip3 install --upgrade pysumstats

Usage

import pysumstats as sumstats

Reading files

s1 = sumstats.SumStats("sumstats1.csv.gz", phenotype='GWASsummary1')

Reading data without sample size column: you will manually have to specify gwas sample size

s2 = sumstats.SumStats("sumstats2.txt.gz", phenotype='GWASsummary2', gwas_n=350492)

Reading data with column names not automatically recognized:

s3 = sumstats.SumStats("sumstats3.csv", phenotype='GWASsummary3',
                              column_names={
                                    'rsid': 'weird_name_for_rsid',
                                    'chr': 'weird_name_for_chr',
                                    'bp': 'weird_name_for_bp',
                                    'ea': 'weird_name_for_ea',
                                    'oa': 'weird_name_for_oa',
                                    'maf': 'weird_name_for_maf',
                                    'b': 'weird_name_for_b',
                                    'se': 'weird_name_for_se',
                                    'p': 'weird_name_for_p',
                                    'hwe': 'weird_name_for_p_hwe',
                                    'info': 'weird_name_for_info',
                                    'n': 'weird_name_for_n',
                                    'eaf': 'weird_name_for_eaf',
                                    'oaf': 'weird_name_for_oaf'})

Performing qc

s1.qc(maf=.01)
s2.qc(maf=.01, hwe=1e-6, info=.9)
s3.qc()  # MAF .01 is the default

Merging sumstats, low_memory option is still experimental so be careful with that

merge1 = s1.merge(s2)

Meta analysis

n_weighted_meta = merge1.meta_analyze(name='meta1', method='samplesize')  # N-weighted meta analysis
ivw_meta = merge1.meta_analyze(name='meta1', method='ivw')  # Standard inverse-variance weighted meta analysis
gwama = merge1.gwama(name='meta1', method='ivw')  # GWAMA as described in Baselmans, et al. (2019)

Additionally supports adding SNP heritabilities as weights

exc_meta = exc.gwama(h2_snp={'ntr_exc': .01, 'ukb_ssoe': .02}, name='exc', method='ivw')

And your own covariance matrix (called cov_Z in most R scripts)

# Either read it from a file:
import pandas as pd
cov_z = pd.read_csv('my_cov_z.csv') # Note it should be pandas dataframe with column names and index names equal to your phenotypes

# Or generate it from a phenotype file yourself:
phenotypes = pd.read_csv('my_phenotype_file.csv')
cov_z = sumstats.cov_matrix_from_phenotype_file(phenotypes, phenotypes=['GWASsummary1', 'GWASsummary2'])

gwama = exc.gwama(cov_matrix=cov_z, h2_snp={'GWASsummary1': .01, 'GWASsummary2': .02}, name='meta1', method='ivw')

See a summary of the result

gwama.describe()

See head of the data

gwama.head()

See head of all chromosomes

gwama.head(n_chromosomes=23)

QQ and Manhattan plots of the result

gwama.manhattan(filename='meta_manhattan.png')
gwama.qqplot(filename='meta_qq.png')

Save the result as csv

exc.save('exc_sumstats.csv')

Save the result as a pickle file (way faster to save and load back into Python)

exc.save('exc_sumstats.pickle')

Merge gwama results with another file:

merged = gwama.merge(s3)

Save prepped files for MR analysis in R:

merged.prep_for_mr(exposure='GWASsummary3', outcome='meta1',
                   filename=['GWAS3-Meta.csv', 'Meta-GWAS3.csv'],
                   p_cutoff=5e-8, bidirectional=True, index=False)

The resulting files will have the following column names, per specification of the MendelianRandomization package in R:

rsid chr bp exposure.A1 exposure.A2 outcome.A1 outcome.A2 exposure.se exposure.b outcome.se outcome.b

Some other stuff:

# See column names of the file
gpc_neuro.columns

# SumStats support for standard indexing is growing:
exc[0]  # Get the full output of the first SNP
exc[:10]  # Get the full output of the first 10 SNPs
exc[:10, 'p']  # Get the p value of the first 10 SNPs
exc['p']  # Get the p values of all SNPs
exc['rs78948828']  # Get the full output of 1 specific rsid
exc[['rs78948828', 'rs6057089', 'rs55957973']]  # Get the full output of multiple specific rsids
exc[['rs78948828', 'rs6057089', 'rs55957973'], 'p']  # Get the p-value for specific rsids

# If for whatever reason you want to do stuff with each SNP individually you can also loop over the entire file
for snp_output in exc:
    if exc['p'] < 5e-8:
        print('Yay significant SNP!')
    # do something


# If you only want to loop over some specific columns, you can
for rsid, b, se, p in exc[['rsid', 'b', 'se', 'p']].values:
    if p < 5e-8:
        print('Yay significant SNP!')