The genetic architecture of fibromyalgia across 2.5 million individuals.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Methods › Harmonization and meta-analysis ↔ fibromyalgia_utils.py, lines 3650–3794 · score 0.94 · single nucleotide variants, multiple rs, allele flips, alt allele, minimal representation, alternate alleles
[2] § Methods › Harmonization and meta-analysis ↔ fibromyalgia_utils.py, lines 2453–2519 · score 0.55 · reference allele, alternate allele, harmonization, rs, plink, matched
[3] § Methods › Genotyping, imputation and quality control › Estonian Biobank ↔ fibromyalgia_utils.py, lines 2453–2519 · score 0.55 · PLINK format, allele frequency, match, filtered, chromosome, variants
[4] § Methods › Genotyping, imputation and quality control › All of Us ↔ fibromyalgia_utils.py, lines 2839–2905 · score 0.52 · minor allele frequency, MAF, sex, GWAS, chromosome, genome

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 3,893 lines · 195 KB · no license · 4 matches

from __future__ import annotations
import os
import polars as pl
import signal
import subprocess
import sys
from functools import cache, reduce
from typing import Any, Iterable
pl.enable_string_cache()
def raise_error_if_on_compute_node(message=None):
"""
Raises an error if the user is on a compute node
Args:
message: Message to print when user is not on a login node; if None,
prints a default message
"""
import re
import socket
if re.search('(nia|nc|nl|ng)[0-9]', socket.gethostname()):
# nia matches the Niagara compute nodes; nc/nl/ng match the Narval ones
import inspect
calling_function = inspect.currentframe().f_back.f_code.co_name
raise RuntimeError(message if message is not None else
f'{calling_function}() needs internet access! Run '
f'once on the login node to download required '
f'files, then re-run here on this compute node')
def check_cluster(cluster):
"""
Check the cluster the user is on.
Returns:
The cluster: "narval" or "niagara". Raises an error if $CLUSTER is not
set to one of those two.
"""
if cluster is None:
raise RuntimeError('The environment variable $CLUSTER is not set; it '
'must be set to "narval" or "niagara"')
if cluster != 'narval' and cluster != 'niagara':
raise RuntimeError(f"The environment variable $CLUSTER is set to "
f"{cluster!r}, but must be set to 'narval' or "
f"'niagara'")
@cache
def get_base_data_directory():
"""
Get the location of the "base" data directory, where data will be stored.
Returns:
The path of the base data directory
"""
cluster = os.environ.get('CLUSTER')
return '/home/wainberg/projects/def-wainberg' if cluster == 'narval' else \
'/scratch/w/wainberg/wainberg' if cluster == 'niagara' else '.'
def plural(string: str, count: int) -> str:
"""
Adds an s to the end of string, unless `count` is 1 or -1.
Args:
string: a string
count: a count
Returns:
`string`, with an s at the end if `count` is 1 or -1
"""
return string if abs(count) == 1 else f'{string}s'
def check_type(variable: Any, variable_name: str,
expected_types: type | tuple[type, ...],
expected_type_name: str) -> None:
"""
Check whether `variable` has the expected type.
Args:
variable: the variable to be checked
variable_name: the name of the variable, used in the error message
expected_types: the expected type or types (specifying int, float, or
bool also implicitly includes their NumPy equivalents)
expected_type_name: the name of the expected type, used in the error
message (e.g. 'a polars DataFrame')
"""
if isinstance(variable, expected_types):
return
if not isinstance(expected_types, tuple):
expected_types = expected_types,
for t in expected_types:
if t is int:
import numpy as np
if isinstance(variable, np.integer):
return
elif t is float:
import numpy as np
if isinstance(variable, np.floating):
return
elif t is bool:
import numpy as np
if isinstance(variable, np.bool_):
return
error_message = (
f'{variable_name} must be {expected_type_name}, but has type '
f'{type(variable).__name__!r}')
raise TypeError(error_message)
def check_types(variable: Iterable[Any],
variable_name: str,
expected_types: type | tuple[type, ...],
expected_type_name: str):
"""
Check whether all elements of `variable` are of the expected type(s).
Args:
variable: the variable to be checked
variable_name: the name of the variable, used in the error message
expected_types: the expected type or types
expected_type_name: the name of the expected type, used in the error
message (e.g. 'polars DataFrames')
"""
if not isinstance(expected_types, tuple):
expected_types = expected_types,
for element in variable:
if not isinstance(element, expected_types):
for t in expected_types:
if t is int:
import numpy as np
if isinstance(variable, np.integer):
break
elif t is float:
import numpy as np
if isinstance(variable, np.floating):
break
elif t is bool:
import numpy as np
if isinstance(variable, np.bool_):
break
else:
error_message = (
f'all elements of {variable_name} must be '
f'{expected_type_name}, but it contains an element of '
f'type {type(element).__name__!r}')
raise TypeError(error_message)
def check_dtype(series: pl.Series,
series_name: str,
expected_dtypes: pl.datatypes.classes.DataTypeClass | str |
tuple[pl.datatypes.classes.DataTypeClass |
str, ...]) -> None:
"""
Check whether `series` has the expected polars dtype.
Args:
series: the polars Series to be checked
series_name: the name of the variable, used in the error message
expected_dtypes: the expected dtype or dtypes. Specify the string
`'integer'` to include all integer dtypes, and
`'floating-point'` to include all floating-point
dtypes.
"""
base_type = series.dtype.base_type()
if not isinstance(expected_dtypes, tuple):
expected_dtypes = expected_dtypes,
for expected_type in expected_dtypes:
if base_type == expected_type or expected_type == 'integer' and \
base_type in pl.INTEGER_DTYPES or \
expected_type == 'floating-point' and \
base_type in pl.FLOAT_DTYPES:
return
if len(expected_dtypes) == 1:
expected_dtypes = str(expected_dtypes[0])
elif len(expected_dtypes) == 2:
expected_dtypes = ' or '.join(map(str, expected_dtypes))
else:
expected_dtypes = ', '.join(map(str, expected_dtypes[:-1])) + \
', or ' + str(expected_dtypes[-1])
error_message = (
f'{series_name} must be {expected_dtypes}, but has data type '
f'{base_type!r}')
raise TypeError(error_message)
class ProcessPool(object):
"""
Like multiprocessing.Pool but 1) child processes ignore KeyboardInterrupts
and 2) apply_async() is called submit() and takes actual *args and **kwargs
rather than a list of args and a dictionary of kwargs (like
ProcessPoolExecutor.submit() from the concurrent.futures module)
Attributes:
pool (multiprocessing.Pool): the underlying multiprocessing.Pool object
"""
def __init__(self, max_concurrent, start_method='forkserver'):
"""
Sets up the multiprocessing.Pool object underlying this ProcessPool.
Args:
max_concurrent: The number of worker processes to be spawned by the
Pool object, i.e. the maximum number of concurrent
processes that will be allowed to run at once.
start_method: How worker processes will be started. Possible
values are 'fork', 'spawn', 'forkserver'. For
details, see docs.python.org/3/library/
multiprocessing.html#contexts-and-start-methods.
"""
import multiprocessing
self.pool = multiprocessing.get_context(start_method) \
.Pool(max_concurrent, initializer=self.ignore_keyboard_interrupts)
@staticmethod
def ignore_keyboard_interrupts():
"""
When provided as an initializer to a multiprocessing.Pool object,
this function tells the Pool to ignore KeyboardInterrupts (i.e. SIGINT)
so that pressing Ctrl + C doesn't kill all your background processes.
"""
signal.signal(signal.SIGINT, signal.SIG_IGN)
def submit(self, func, *args, **kwargs):
"""
Submits a function to the process pool. A thin wrapper over
Pool.apply_async() that allows the user to pass actual *args and
**kwargs rather than a list of args and a dictionary of kwargs.
Args:
func: the function to be submitted
*args: positional arguments to be passed to function
**kwargs: keyword arguments to be passed to function
Returns:
The multiprocessing.pool.AsyncResult object returned by
Pool.apply_async().
"""
return self.pool.apply_async(func, args, kwargs)
@cache
def get_process_pool(max_concurrent, start_method='forkserver'):
"""
Creates a ProcessPool for a given value of max_concurrent and start_method,
which due to the @cache decorator will persist across multiple calls to
functions like run_background() or run_function_background(), so long as
they call this function with the same max_concurrent and start_method.
Args:
max_concurrent: The number of worker processes to be spawned by the
ProcessPool, i.e. the maximum number of concurrent
processes that will be allowed to run at once.
start_method: How worker processes will be started. Possible
values are 'fork', 'spawn', 'forkserver'. For
details, see docs.python.org/3/library/
multiprocessing.html#contexts-and-start-methods.
Returns:
A ProcessPool for the given values of max_concurrent and start_method.
A new ProcessPool will be created the first time this function is run
for a given value of max_concurrent and start_method, which will then
be cached for all subsequent calls with the same max_concurrent and
start_method.
"""
return ProcessPool(max_concurrent, start_method=start_method)
@cache
def cython_inline(code, debug=False, boundscheck=None, cdivision=True,
initializedcheck=None, wraparound=False,
warn_undeclared=True, include_dirs=None, libraries=None,
clang=False, extra_compiler_flags=None,
extra_linker_flags=None, verbose=False,
**other_cython_settings):
"""
A drop-in replacement for `cython.inline()` that supports cimports. It
turns on the major Cython optimizations (`boundscheck=False`,
`cdivision=True`, `initializedcheck=False`, `wraparound=False`) and sets
`language_level=3` for full Python 3 compatibility.
Args:
code: a string of Cython code to compile
debug: whether to turn off most compiler optimizations (`-Ofast`,
`-funroll-loops`) and turn on debug compilation (`-g -Og`) and
Cython's boundscheck and initializedcheck
boundscheck: whether to perform array bounds checking when indexing;
always affects array/memoryview indexing, but also affects
list, tuple, and string indexing when wraparound=False
cdivision: whether to use C-style rather than Python-style division and
remainder operations; disabling leads to a ~35% speed
penalty for these operations
initializedcheck: whether to check whether memoryviews and C++ classes
are initialized before using them
wraparound: whether to support Python-style negative indexing
warn_undeclared: whether to warn about undeclared variables (i.e. those
without a cdef type declaration)
include_dirs: an optional tuple of include directories of libraries to
link against; `np.get_include()` will always be included
libraries: an optional tuple of libraries to link against,
e.g. ('hdf5',)
clang: whether to compile with clang instead of GCC
extra_compiler_flags: a tuple of extra compiler flags to use on top of
the defaults
extra_linker_flags: a tuple of extra linker flags to use on top of the
defaults
verbose: if True, print Cython's compilation logs
**other_cython_settings: other Cython settings, which will be written
into the source code as #cython compiler
directives
Returns:
The {function_name: function} dictionary of compiled functions that
would be returned by cython.inline().
"""
import numpy as np
from hashlib import md5
from inspect import getmembers
from textwrap import dedent
# ~ is read-only on Niagara compute nodes, so build in CYTHON_CACHE_DIR in
# scratch instead
cython_cache_dir = os.path.abspath(os.environ.get(
'CYTHON_CACHE_DIR', os.path.expanduser('~/.cython')))
os.makedirs(cython_cache_dir, exist_ok=True)
# If `boundscheck` and/or `initializedcheck` are `None`, turn them on if
# `debug` is `True`, and turn them off otherwise
if boundscheck is None:
boundscheck = debug
if initializedcheck is None:
initializedcheck = debug
# Remove extra levels of indentation from the code string (since it's
# usually defined inside a function, so there's at least one extra level of
# indentation that would cause a syntax error if not removed) and remove
# any leading newlines if present (since when users define code strings,
# they usually use triple-quoted strings, and the code usually doesn't
# start until the line after the three opening quotes, leading to a single
# leading newline)
settings = dict(language_level=3, boundscheck=boundscheck,
cdivision=cdivision, initializedcheck=initializedcheck,
wraparound=wraparound,
**{'warn.undeclared': warn_undeclared})
settings.update(other_cython_settings)
code = ''.join(f'#cython: {setting_name}={setting}\n'
for setting_name, setting in settings.items()) + \
f'#distutils: define_macros=NPY_NO_DEPRECATED_API=' \
f'NPY_1_7_API_VERSION\n' + dedent(code)
# Define build options and environment variables
if include_dirs is not None:
include_dirs = (np.get_include(),) + include_dirs
else:
include_dirs = np.get_include(),
include_dirs = \
'[' + ', '.join(f'{include_dir!r}'
for include_dir in include_dirs) + ']'
if libraries is not None:
libraries = \
f'[' + ', '.join(f'{library!r}' for library in libraries) + ']'
narval = os.environ.get('CLUSTER') == 'narval'
niagara = os.environ.get('CLUSTER') == 'niagara'
march = 'znver2' if narval else 'skylake' if niagara else 'native'
extra_compiler_flags = \
''.join(f', {flag!r}' for flag in extra_compiler_flags) \
if extra_compiler_flags is not None else ''
extra_linker_flags = \
''.join(f', {flag!r}' for flag in extra_linker_flags) \
if extra_linker_flags is not None else ''
build_options = f'''
language='c++',
include_dirs={include_dirs},
libraries={libraries},
extra_compile_args=['-g', '-Og', '-march={march}', '-fopenmp',
'-Wall', '-Wextra', '-Wpedantic',
'-Werror', '-Wno-maybe-uninitialized',
'-Wno-ignored-qualifiers'
{extra_compiler_flags}]
if {debug} else ['-Ofast', '-march={march}',
'-funroll-loops', '-fopenmp', '-Wall',
'-Wextra', '-Wpedantic', '-Werror',
'-Wno-maybe-uninitialized',
'-Wno-ignored-qualifiers'
{extra_compiler_flags}],
extra_link_args=['-fopenmp'{extra_linker_flags}] if {debug}
else
['-Ofast', '-fopenmp'{extra_linker_flags}]'''
build_environment_variables = \
f'CXX={"clang++" if clang else "g++"}{" CFLAGS=-g " if debug else ""}'
# Make a short alphabetic module name by taking the MD5 hash of the
# concatenated code, build options, and build environment variables, and
# converting the hexadecimal digits to letters (0 -> a, 1 -> b, ...,
# 9 --> j, a --> k, ..., f --> p)
module_name = ''.join(chr(ord(c) + (49 if c <= '9' else 10))
for c in md5((code + build_options +
build_environment_variables)
.encode('utf-8')).hexdigest())
code_file = os.path.join(cython_cache_dir, f'{module_name}.pyx')
# Try to import the module; build it if it does not exist
sys.path.append(cython_cache_dir)
try:
module = __import__(module_name)
except ModuleNotFoundError:
# Create the code file
with open(code_file, 'w') as f:
# noinspection PyTypeChecker
print(code, file=f)
# Write the build script to a temp file based on the module name
build_file = os.path.join(cython_cache_dir, f'{module_name}_build.py')
with open(build_file, 'w') as f:
build_script = dedent(f'''
from setuptools import Extension, setup
from Cython.Build import cythonize
setup(name='{module_name}', ext_modules=cythonize([
Extension('{module_name}', ['{code_file}'],
{build_options})],
build_dir='{cython_cache_dir}'))''')
# noinspection PyTypeChecker
print(build_script, file=f)
# Build the code (note: `sys.executable` is the location of Python)
run(f'cd {cython_cache_dir} && {build_environment_variables} '
f'{sys.executable} {build_file} build_ext --inplace'
f'{"" if verbose else " > /dev/null"}')
# Remove the temp file
os.unlink(build_file)
# Try again
module = __import__(module_name)
finally:
# noinspection PyInconsistentReturns
sys.path = sys.path[:-1]
# Create a dict of all the Cython functions defined in the module
function_dict = {function_name: function
for function_name, function in getmembers(module)
if repr(function).startswith('<cyfunction')}
# Return the dict of Cython functions
return function_dict
def print_df(df, num_rows=-1, num_columns=-1):
"""
Prints the entirety of a polars DataFrame without truncating.
Args:
df: the DataFrame to print
num_rows: the number of rows to print (-1 to print all rows)
num_columns: the number of columns to print (-1 to print all columns)
"""
with pl.Config(tbl_rows=num_rows, tbl_cols=num_columns):
print(df)
def map_df(df, map_col, other_df, key_col, value_col, *,
retain_missing=False):
"""
Maps df[map_col] based on the mapping other_df[key_col] ->
other_df[value_col].
In other words, for each element of df[map_col], check if it's in
other_df[key_col], and if so, replace it with the corresponding entry of
other_df[value_col].
Equivalent to df.with_columns(pl.col(map_col).replace_strict(dict(zip(
other_df[key_col], other_df[value_col])), default=pl.first() if
retain_missing else None)).
Implementation detail: uses a join, but prefixes other_df's key_col and
value_col with "__MAP_DF_" to handle the possibility that map_col might
have the same name as key_col or value_col.
Args:
df: a polars DataFrame
map_col: a column in df
other_df: another polars DataFrame
key_col: a column in other_df with the mapping keys; all values must be
unique, although this is not checked, for speed
value_col: a column in other_df with the mapping values
retain_missing: if False, sets elements of map_col that don't appear in
key_col to null; if True, leaves them unchanged
Returns:
df with map_col transformed so that each of its values that are in
key_col are transformed to the corresponding value in value_col.
"""
if key_col == value_col:
raise ValueError(f'Both key_col and value_col are set to the column '
f'name "{key_col}"')
if isinstance(df, pl.LazyFrame):
other_df = other_df.lazy()
if isinstance(other_df, pl.LazyFrame):
df = df.lazy()
prefix = '__MAP_DF_'
df = df.join(other_df.select(pl.col(key_col, value_col)
.name.prefix(prefix)),
left_on=map_col, right_on=prefix + key_col, how='left')
if retain_missing:
df = df.with_columns(pl.col(prefix + value_col)
.fill_null(pl.col(map_col)))
df = df.with_columns(pl.col(prefix + value_col).alias(map_col)) \
.drop(prefix + value_col)
return df
def thread_string(num_threads):
"""
Generates a string for setting relevant environment variables to limit the
number of threads for a called process.
Args:
num_threads: The number of threads.
Returns:
A string for setting the relevant environment variables to num_threads.
"""
return f'export MKL_NUM_THREADS={num_threads}; ' \
f'export OMP_NUM_THREADS={num_threads}; ' \
f'export OPENBLAS_NUM_THREADS={num_threads}; ' \
f'export NUMEXPR_MAX_THREADS={num_threads}; ' \
if num_threads is not None else ''
def run(cmd, *, log_file=None, unbuffered=False, pipefail=True,
num_threads=None, **kwargs):
"""
Runs a bash code segment interactively.
Args:
cmd: the command to be run
log_file: a filename to log stdout/stderr to, in addition to printing
unbuffered: set to True for unbuffered I/O (currently broken when cmd
contains multiple commands)
pipefail: set to False when piping commands to head to avoid errors
num_threads: set to a positive integer to limit how many threads cmd
uses, or to None to not limit the number of threads
**kwargs: passed on to subprocess.run()
Returns:
The CompletedProcess object returned by subprocess.run().
"""
run_kwargs = dict(check=True, shell=True, executable='/bin/bash')
run_kwargs.update(**kwargs)
return subprocess.run(
f'{thread_string(num_threads)}'
f'set -eu{"o pipefail" if pipefail else ""}; '
f'{"stdbuf -i0 -o0 -e0 " if unbuffered else ""}{cmd}'
f'{f" 2>&1 | tee {log_file}" if log_file is not None else ""}',
**run_kwargs)
def read_csv_from_command(cmd, *, run_kwargs={}, **kwargs):
"""
Read a columnar file with polars from the output of a bash command.
Args:
cmd: the bash command
run_kwargs: keyword arguments to utils.run()
**kwargs: keyword arguments to pl.read_csv()
Returns:
A polars DataFrame with the contents of the columnar file produced by
the bash command.
"""
from io import BytesIO
# noinspection PyTypeChecker,PyArgumentList
return pl.read_csv(BytesIO(run(cmd, stdout=subprocess.PIPE, **run_kwargs)
.stdout), **kwargs)
def read_csv_delim_whitespace(whitespace_delimited_file, **kwargs):
"""
Reads a whitespace-delimited file with polars, mimicking the behavior of
delim_whitespace=True in pandas.read_csv().
Explanation of the sed command:
s/^[[:space:]]\\+// removes leading whitespace
s/[[:space:]]\\+$// removes trailing whitespace
s/[[:space:]]\\+/\t/g replaces internal whitespace with tabs
Args:
whitespace_delimited_file: the whitespace-delimited file
**kwargs: keyword arguments to pl.read_csv()
Returns:
A polars DataFrame with the contents of the whitespace-delimted file.
"""
if not os.path.exists(whitespace_delimited_file):
error_message = \
f'No such file or directory: {whitespace_delimited_file}'
raise FileNotFoundError(error_message)
if whitespace_delimited_file.endswith('.gz'):
whitespace_delimited_file = f'<(zcat {whitespace_delimited_file})'
return read_csv_from_command(
f"sed 's/^[[:space:]]\\+//; s/[[:space:]]\\+$//; "
f"s/[[:space:]]\\+/\t/g' {whitespace_delimited_file}",
separator='\t', **kwargs)
def savefig(filename, *, despine=False, **kwargs):
"""
A drop-in replacement for plt.savefig(). Saves a matplotlib figure to
filename, but additionally:
- removes the right and top spines for a cleaner-looking plot (if
despine=True)
- saves with dpi=450 for higher-resolution plots
- saves with bbox_inches='tight' and pad_inches=0 to avoid extra whitespace
around plots
- saves PDF files with transparent=True so transparency info isn't
discarded
- closes the plot with plt.close() to avoid memory leaks
Args:
filename: the filename to save the matplotlib figure to
despine: whether to remove the right and top spines using seaborn's
despine() function
**kwargs:
"""
import matplotlib.pyplot as plt
if despine:
spines = plt.gca().spines
spines['top'].set_visible(False)
spines['right'].set_visible(False)
all_kwargs = dict(dpi=450, bbox_inches='tight', pad_inches=0,
transparent=filename.endswith('pdf'))
all_kwargs.update(kwargs)
plt.savefig(filename, **all_kwargs)
plt.close()
def manhattan_plot(sumstats, *, clumped_variants=None, genome_build=None,
SNP_col='SNP', chrom_col='CHROM', bp_col='BP',
ref_col='REF', alt_col='ALT', p_col='P',
max_p=1e-3, max_p_to_label=5e-8,
p_thresholds={5e-8: ('#D62728', 'dashed')},
label_overrides={}, label_kwargs={},
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations'):
"""
Make a Manhattan plot of the variants in sumstats; highlight lead variants.
If clumped_variants is not None, highlight the LD clumps as well.
If genome_build is not None, label each lead variant with its nearest
coding gene(s) from that genome build.
After running, call savefig() to save the plot, or customize further first.
Style inspired by ncbi.nlm.nih.gov/pmc/articles/PMC6481311/figure/F1.
Requires the textalloc package (github.com/ckjellson/textalloc) to make
sure the gene labels don't overlap. Install it with:
pip install --no-deps --no-build-isolation textalloc
Args:
sumstats: the summary statistics to plot
clumped_variants: a DataFrame of clumped variants from ld_clump() to
highlight lead variants on the Manhattan plot, or
None to skip
genome_build: if not None, a genome build used to label each lead
variant with its nearest coding gene(s)
SNP_col: the name of the variant ID column in sumstats
chrom_col: the name of the chromosome column in sumstats
bp_col: the name of the base-pair position column in sumstats
ref_col: the name of the reference allele column in sumstats
alt_col: the name of the alternate allele column in sumstats
p_col: the name of the p-value column in sumstats
max_p: defines the bottom y limit of the Manhattan plot; variants with
p-values >= this value will not be plotted
max_p_to_label: variants with p-values >= this value will not be
labeled; only has an effect if genome_build is not None
p_thresholds: a dictionary of p-value thresholds to draw horizontal
lines at; values are tuples of (line color, line style)
label_overrides: a dictionary of {gene_name: label} used to override
specific gene labels.
label_kwargs: a dictionary of keyword arguments to be passed to
`textalloc.allocate()` when adding gene labels to
control the text properties, such as:
- `textcolor`/`textsize`: the text color and size
- `x_scatter`/`y_scatter`: the x/y coordinates of
points in the scatter plot, to repel labels away
from. Defaults to all points in the plot.
- `min_distance`/`max_distance`: the minimum and
maximum distances from each point to its label, as
a proportion of the width of the x-axis. Defaults
to 0 and 0.02, instead of textalloc's defaults of
0.015 and 0.2
- `draw_lines`: whether to draw lines between each
label and its corresponding point. Defaults to
`False`, instead of textalloc's default of
`True`.
See github.com/ckjellson/textalloc#parameters for the
full list of possible arguments.
gene_annotation_dir: the directory where the coding gene locations
returned by get_coding_genes() will be cached.
Must be run on the login node to generate this
cache, if it doesn't exist. Because generating the
cache can take a long time, you will probably want
to leave this argument at its default value.
"""
import matplotlib.pyplot as plt
import numpy as np
# noinspection PyUnresolvedReferences
import textalloc
if sumstats.is_empty():
raise ValueError(f'sumstats is empty!')
if SNP_col not in sumstats:
raise ValueError(f'{SNP_col!r} not in sumstats; specify SNP_col')
if chrom_col not in sumstats:
raise ValueError(f'{chrom_col!r} not in sumstats; specify chrom_col')
if bp_col not in sumstats:
raise ValueError(f'{bp_col!r} not in sumstats; specify bp_col')
if p_col not in sumstats:
raise ValueError(f'{p_col!r} not in sumstats; specify p_col')
if genome_build is not None:
check_valid_genome_build(genome_build=genome_build)
# Subset sumstats to p < max_p; standardize chromosome names; alternate
# each chromosome in a different color (dark blue then light blue); offset
# each chromosome by the cumulative number of base pairs from the start of
# chr1 to the start of that chromosome
sumstats = sumstats \
.filter(pl.col(p_col) < max_p) \
.with_columns(pl.col(chrom_col)
.pipe(standardize_chromosomes, omit_chr_prefix=True)) \
.with_columns(color=pl.col(chrom_col).cast(pl.Categorical)
.to_physical().mod(2)
.replace_strict({0: '#0A6FA5', 1: '#008FCD'}))
# Get the total number of bps to the start and end of each chromosome
cumulative_bp = sumstats \
.group_by(chrom_col, maintain_order=True) \
.agg(pl.max(bp_col)) \
.with_columns(end=pl.col(bp_col).cum_sum()) \
.with_columns(start=pl.col.end.shift().fill_null(0)) \
.drop(bp_col)
# Get x and y coordinates to plot
sumstats = sumstats \
.join(cumulative_bp, on=chrom_col, how='left') \
.with_columns(x=pl.col.start + pl.col(bp_col),
y=-pl.col(p_col).log10())
# Plot horizontal lines indicating significance thresholds
for p_threshold, (color, linestyle) in p_thresholds.items():
plt.axhline(y=-np.log10(p_threshold), color=color, linestyle=linestyle,
zorder=-1)
# Overplot three times: first non-clumped variants, then clumped
# variants, then lead variants. Rasterize non-clumped and (optionally)
# clumped variants for quick plotting and rendering.
w, h = plt.gcf().get_size_inches()
plt.gcf().set_size_inches(1.2 * w, h)
if clumped_variants is None:
plt.scatter(sumstats['x'], sumstats['y'], c=sumstats['color'], s=1,
rasterized=True)
else:
non_clump_sumstats = sumstats.filter(
~pl.col(SNP_col).is_in(clumped_variants[SNP_col]))
clump_sumstats = sumstats.filter(
pl.col(SNP_col).is_in(clumped_variants[SNP_col]),
~pl.col(SNP_col).is_in(clumped_variants[f'{SNP_col}_lead']))
join_columns = SNP_col, chrom_col, bp_col, ref_col, alt_col
lead_sumstats = clumped_variants \
.filter('is_lead') \
.select(join_columns) \
.join(sumstats, on=join_columns, how='left')
plt.scatter(non_clump_sumstats['x'], non_clump_sumstats['y'],
c=non_clump_sumstats['color'], s=1, rasterized=True)
plt.scatter(clump_sumstats['x'], clump_sumstats['y'], c='#DBA756', s=1)
plt.scatter(lead_sumstats['x'], lead_sumstats['y'],
c='#DBA756', marker='D', edgecolors='k', s=10)
if genome_build is not None:
lead_sumstats = lead_sumstats \
.pipe(lambda df: df.join(
df.filter(pl.col(p_col) < max_p_to_label)
.pipe(get_nearest_gene, genome_build=genome_build,
SNP_col=SNP_col, chrom_col=chrom_col, bp_col=bp_col,
gene_annotation_dir=gene_annotation_dir),
on=(SNP_col, chrom_col, bp_col), how='left')) \
.drop_nulls('gene')
labels = lead_sumstats['gene']
if label_overrides:
labels = \
labels.list.eval(pl.element().replace(label_overrides))
labels = labels.list.join(', ')
default_label_kwargs = dict(
ax=plt.gca(), x=lead_sumstats['x'], y=lead_sumstats['y'],
text_list=labels, x_scatter=sumstats['x'],
y_scatter=sumstats['y'], min_distance=0, max_distance=0.02,
draw_lines=False)
label_kwargs = default_label_kwargs | label_kwargs \
if label_kwargs is not None else default_label_kwargs
textalloc.allocate(**label_kwargs)
# Put each chrom's label in the chrom's center; hide chr17/19/21 labels
plt.xticks((cumulative_bp['start'] + cumulative_bp['end']) / 2,
cumulative_bp[chrom_col]
.replace({'17': '', '19': '', '21': ''}))
# Set x and y labels and limits, resize
plt.xlabel('Chromosome')
plt.ylabel('-log$_{10}$(p)')
padding = 10_000_000 # avoid clipping SNPs at the far left or right
plt.xlim(sumstats['x'].min() - padding, sumstats['x'].max() + padding)
plt.ylim(bottom=-np.log10(max_p))
def qqplot(ps, *, equal_aspect=True, **kwargs):
"""
Generates a quantile-quantile (Q-Q) plot of p-values.
Inspired by qqplot() from the qmplot package.
Args:
ps: a polars Series or 1D NumPy array of p-values to plot
equal_aspect: if True, calls ax.set_aspect('equal') so the x and y axes
use the same number of inches per unit increase in x and
y
**kwargs: passed to ax.scatter()
"""
#
import matplotlib.pyplot as plt
import numpy as np
if ps.is_empty():
raise ValueError(f'ps is empty!')
ax = kwargs.pop('ax') if 'ax' in kwargs else plt.gca()
rasterized = kwargs.pop('rasterized') if 'rasterized' in kwargs else True
ppoints = lambda n, a=0.5: (np.arange(n) + 1 - a) / (n + 1 - 2 * a)
ax.scatter(-np.log10(ppoints(len(ps))),
-np.log10(np.sort(ps).clip(5e-324)),
edgecolors='none', rasterized=rasterized, **kwargs)
xmax = ax.get_xlim()[1]
ax.plot([0, xmax], [0, xmax], c='k', zorder=-1)
if equal_aspect:
ax.set_aspect('equal')
ax.set_xlim(0)
ax.set_ylim(0)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
def standardize(array):
"""
Standardize a polars Series, DataFrame or expression or a NumPy array to
zero mean and unit variance (using 1 delta degree of freedom), columnwise.
Args:
array: the Series, DataFrame, expression or array to standardize
Returns:
The standardized Series, DataFrame, expression or array.
"""
if isinstance(array, (pl.Series, pl.Expr)):
return (array - array.mean()) / array.std()
if isinstance(array, pl.DataFrame):
return array.with_columns((pl.selectors.numeric() -
pl.selectors.numeric().mean()) /
pl.selectors.numeric().std())
import numpy as np
if not isinstance(array, np.ndarray):
raise ValueError(f'array must be a polars Series, DataFrame or '
f'expression or a NumPy array!')
return (array - np.mean(array, axis=0)) / np.std(array, ddof=1, axis=0)
def inflation_factor(pvalues):
"""
Calculates the genomic inflation factor from a vector of p-values.
Args:
pvalues: a polars Series or expression or NumPy array of p-values
Returns:
A single floating-point number with the genomic inflation factor.
"""
from scipy.special import chdtri # chdtri(df, p) == chi2.isf(p, df)
if isinstance(pvalues, (pl.Series, pl.Expr)):
return chdtri(1, pvalues.median()) / chdtri(1, 0.5)
import numpy as np
if not isinstance(pvalues, np.ndarray):
raise ValueError('pvalues must be a polars Series or expression or a '
'NumPy array!')
return chdtri(1, np.median(pvalues)) / chdtri(1, 0.5)
def inverse_normal_transform(values, *, c=3 / 8):
"""
Calculates the rank-based inverse normal transform of a polars Series or
expression or 1D NumPy array.
Args:
values: a polars Series or expression or 1D NumPy array.
c: a parameter of the transformation: a fractional shift to the ranks
that's applied before transforming. By default, we use the Blom
transform (c = 3/8). For more details on c, see:
ncbi.nlm.nih.gov/pmc/articles/PMC2921808/#S2title
Returns:
The rank-based inverse normal transform of values.
"""
from scipy.special import ndtri # ntdri(x) == norm.ppf(x)
if isinstance(values, (pl.Series, pl.Expr)):
rank = values.rank()
transformed_rank = (rank - c) / \
(rank.len() - rank.null_count() - 2 * c + 1)
else:
import numpy as np
from scipy.stats import rankdata
if not isinstance(values, np.ndarray) or values.ndim != 1:
raise ValueError('values must be a polars Series or expression or '
'a 1D NumPy array!')
rank = rankdata(values)
transformed_rank = (rank - c) / ((~np.isnan(rank)).sum() - 2 * c + 1)
return ndtri(transformed_rank)
@polars_numpy_autoconvert()
def z_to_p(z_scores, *, high_precision=False):
"""
Converts a polars Series or NumPy array of z-scores to p-values
Args:
z_scores: the polars Series or NumPy array of z-scores
high_precision: if True, uses R's pnorm function for high-precision
output - important for very small p-values
Returns:
The corresponding p-values; same type as z_scores.
"""
import numpy as np
if high_precision:
from ryp import to_py, to_r
to_r(np.abs(z_scores), 'abs.z')
pvalues = 2 * np.exp(np.float128(to_py(
'pnorm(abs.z, lower.tail=False, log.p=True)')))
else:
from scipy.special import ndtr
pvalues = 2 * ndtr(-np.abs(z_scores)) # ndtr(-x) == norm.sf(x)
return pvalues
def p_to_abs_z(pvalues):
"""
Converts a polars Series or expression or NumPy array of p-values to
abs(z-scores)
Args:
pvalues: the polars Series or expression or NumPy array of p-values
Returns:
The corresponding abs(z-scores); same type as pvalues.
"""
from scipy.special import ndtri
abs_z_scores = -ndtri(pvalues / 2) # -ndtri(x) == norm.isf(x)
return abs_z_scores
def standardize_chromosomes(chromosome, *, return_numeric=False,
omit_chr_prefix=False):
"""
Standardizes the chromosome name(s) in chromosome, which can be an integer,
string, polars Series or polars expression.
Args:
chromosome: the chromosome name(s) to standardize
return_numeric: whether to return chromosomes as numbers
omit_chr_prefix: whether to leave out the 'chr' prefix. Mutually
exclusive with return_numeric.
Returns:
The corresponding standardized chromosome name(s):
- 'chr1' to 'chr22': return as-is
- 1 to 22 or '1' to '22': convert to 'chr1' to 'chr22'
- 23, '23', 'chr23', 'X': convert to 'chrX'
- 24, '24', 'chr24', 'Y': convert to 'chrY'
- 'M' or 'MT': convert to 'chrM'
- 25, '25', 'chr25': disallow; could refer to either chrM or chrXY (the
pseudoautosomal region of the X and Y chromosomes)
Or, if return_numeric=True, disallow chrM and its aliases and convert
everything to a number between 1 and 24, inclusive.
Or, if omit_chr_prefix=True, leave out the 'chr' prefix.
"""
import numpy as np
if return_numeric and omit_chr_prefix:
raise ValueError(f'Only one of return_numeric and omit_chr_prefix can '
f'be True!')
check_type(chromosome, 'chromosome', (int, str, pl.Series, pl.Expr),
'an int, str, or polars Series or expression')
if isinstance(chromosome, pl.Expr):
return chromosome.map_batches(lambda col: standardize_chromosomes(
col, return_numeric=return_numeric,
omit_chr_prefix=omit_chr_prefix))
elif isinstance(chromosome, str) or \
isinstance(chromosome, (int, np.integer)):
if isinstance(chromosome, str):
original_chromosome = chromosome
chromosome = chromosome.removeprefix('chr').replace('X', '23') \
.replace('Y', '24').replace('MT', 'M')
if not ((chromosome.isdigit() and 1 <= int(chromosome) <= 24) or
(chromosome == 'M' and not return_numeric)):
raise ValueError(f'Invalid chromosome '
f'{original_chromosome!r}!')
else:
if chromosome not in range(1, 25):
raise ValueError(f'chromosome == "{chromosome}" but should be '
f'between 1 and 24 (chrY) inclusive!')
if return_numeric:
# noinspection PyTypeChecker
return int(chromosome)
else:
chromosome = str(chromosome).replace('23', 'X').replace('24', 'Y')
return chromosome if omit_chr_prefix else f'chr{chromosome}'
else:
check_dtype(chromosome, 'chromosome',
(pl.String, pl.Categorical, pl.Enum, 'integer'))
if chromosome.dtype in pl.INTEGER_DTYPES:
mask = chromosome.is_in(range(1, 25))
if not mask.all():
error_message = (
'chroms were specified as integers, but some were not '
'between 1 and 24 (chrY) inclusive!')
raise ValueError(error_message)
chromosome = chromosome.cast(pl.String)
else:
# Remove this if statement once polars allows Categorical
# .str.len_bytes(): github.com/pola-rs/polars/issues/9773
if chromosome.dtype == pl.Categorical or \
chromosome.dtype == pl.Enum:
chromosome = chromosome.cast(pl.String)
chromosome = chromosome.str.replace('^chr', '') \
.str.replace('X', '23').str.replace('Y', '24') \
.str.replace('MT', 'M')
valid = chromosome.cast(pl.Int8, strict=False).is_in(range(1, 25))
if not return_numeric:
valid |= chromosome == 'M'
if not valid.all():
raise ValueError('Some chromosomes are invalid!')
if return_numeric:
return chromosome.cast(pl.Int8)
else:
chromosome = chromosome.replace({'23': 'X', '24': 'Y'})
return chromosome if omit_chr_prefix else 'chr' + chromosome
def load_alias_to_gene_map(gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations'):
"""
Load a map from old gene names ("aliases") to their current gene names.
Remove "ambiguous" aliases that map to multiple current gene names (like
ACSM2, which maps to both ACSM2A and ACSM2B). This is done with
pl.col.alias.is_unique() below.
Also remove aliases that are gene names themselves, e.g. TTLL5 is an alias
of TTLL10, but TTLL5 is also a gene name. This can happen when a gene
"splits" into two genes. This is done with the two ~pl.col.alias.is_in
conditions below.
Args:
gene_annotation_dir: the directory where the alias-to-gene map will be
cached. Must be run on the login node to generate
this cache, if it doesn't exist. Because
generating the cache can take a long time, you
will probably want to leave this argument at its
default value.
Returns:
A two-column DataFrame with old gene names in the "alias" column and
their current gene names in the "gene" column.
"""
gene_alias_file = f'{gene_annotation_dir}/gene_aliases.tsv'
if not os.path.exists(gene_alias_file):
raise_error_if_on_compute_node()
os.makedirs(gene_annotation_dir, exist_ok=True)
run(f'curl -fsSL https://ftp.ebi.ac.uk/pub/databases/genenames/new/'
f'tsv/locus_groups/protein-coding_gene.txt | cut -f2,9,11 > '
f'{gene_alias_file}')
Ensembl_genes = get_Ensembl_map(gene_annotation_dir=gene_annotation_dir,
no_unaliasing=True) \
.filter(pl.col.most_recent_Ensembl_version ==
pl.col.most_recent_Ensembl_version.first()) \
['gene_name']
alias_to_gene = pl.scan_csv(gene_alias_file, separator='\t') \
.select(alias=pl.col.alias_symbol.str.split('|')
.list.concat(pl.col.prev_symbol), gene='symbol') \
.drop_nulls() \
.explode('alias') \
.filter(pl.col.alias.is_unique(), ~pl.col.alias.is_in(pl.col.gene),
~pl.col.alias.is_in(Ensembl_genes)) \
.collect()
return alias_to_gene
def unalias(df, gene_col, *,
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations'):
"""
"Unaliases" the genes in df's gene_col by mapping old gene names to their
current gene names, according to the map returned by
load_alias_to_gene_map().
Before matching a gene list from a third-party dataset
Before matching gene lists from two third-party datasets to each other, run
them both through this function. You don't have to run this on the result
of Ensembl_to_gene(), though.
Args:
df: a polars DataFrame
gene_col: the name of the column with the gene names to be unaliased
gene_annotation_dir: the directory where the alias-to-gene map returned
by load_alias_to_gene_map() will be cached. Must
be run on the login node to generate this cache,
if it doesn't exist. Because generating the cache
can take a long time, you will probably want to
leave this argument at its default value.
Returns:
df with each matching gene name in gene_col mapped to its alias. gene
names not matching any of the old gene names in
load_alias_to_gene_map() are left as-is.
"""
return map_df(df, gene_col, load_alias_to_gene_map(
gene_annotation_dir=gene_annotation_dir),
key_col='alias', value_col='gene', retain_missing=True)
def get_Ensembl_map(ENSP=False,
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations',
no_unaliasing=False):
"""
Gets a map from ENSGs (or ENSPs, if ENSP=True) to their most recent gene
symbols in the Ensembl database, according to the map returned by
get_Ensembl_map(). The returned gene names are unaliased, so you don't have
to run unalias() on them.
Ensembl IDs retired on or before release-42 (in 2006) are not included
since these releases lack GTF files on the Ensembl website, and these IDs
would be unlikely to be encountered in modern genomics data anyway.
Args:
ENSP: whether to convert ENSPs (protein IDs) instead of ENSGs (gene
IDs)
gene_annotation_dir: the directory where the Ensembl to gene symbol map
will be cached. Must be run on the login node to
generate this cache, if it doesn't exist. Because
generating the cache can take a long time, you
will probably want to leave this argument at its
default value.
no_unaliasing: if True, gene names will not be unalised
Returns:
A two-column polars DataFrame with an "Ensembl_ID" column containing
each of the Ensembl IDs in the Ensembl database and a "gene_name"
column containing their gene names.
"""
Ensembl_ID_type = 'ENSP' if ENSP else 'ENSG'
mapping_file = os.path.join(gene_annotation_dir,
f'{Ensembl_ID_type}_to_gene_name.tsv')
if not os.path.exists(mapping_file):
print(f'Generating "{mapping_file}"...')
raise_error_if_on_compute_node()
os.makedirs(gene_annotation_dir, exist_ok=True)
latest_Ensembl_version = max(map(int, run(
'curl -fsSL https://ftp.ensembl.org/pub/current_gtf/'
'homo_sapiens/ | grep -oP "GRCh38\\.\\K[0-9]{3}" | uniq',
stdout=subprocess.PIPE).stdout.split()))
get_build = lambda i: "GRCh38" if i in range(76, 82) else "GRCh37" \
if i in range(55, 76) else "NCBI36"
GTF_basenames = [
basename
for i in range(latest_Ensembl_version, 81, -1) for
basename in (
[f'release-{i}/gtf/homo_sapiens/'
f'Homo_sapiens.GRCh38.{i}.chr_patch_hapl_scaff',
f'grch37/release-{i}/gtf/homo_sapiens/'
f'Homo_sapiens.GRCh37.{i}.chr_patch_hapl_scaff']
if i in (87, 85, 82) else
[f'release-{i}/gtf/homo_sapiens/'
f'Homo_sapiens.GRCh38.{i}.chr_patch_hapl_scaff'])
] + [
f'release-{i}/gtf/homo_sapiens/Homo_sapiens.{get_build(i)}.{i}'
for i in range(81, 47, -1)
] + [
'release-47/gtf/Homo_sapiens.NCBI36.47',
'release-46/homo_sapiens_46_36h/data/gtf/Homo_sapiens.NCBI36.46',
'release-44/homo_sapiens_44_36f/data/gtf/Homo_sapiens.NCBI36.44',
'release-43/homo_sapiens_43_36e/data/gtf/Homo_sapiens.NCBI36.43']
sed_command = \
's/.*gene_name "(\\S*)".*protein_id "(\\S*)".*/\\2\\t\\1/p' \
if ENSP else \
's/.*gene_id "(\\S*)".*gene_name "(\\S*)".*/\\1\\t\\2/p'
cache_dir = os.path.join(gene_annotation_dir,
f'{Ensembl_ID_type}_by_version')
os.makedirs(cache_dir, exist_ok=True)
version_mapping_files = [os.path.join(
cache_dir, f'{basename.split("/")[1].replace("-", "_")}_grch37.tsv'
if basename.startswith('grch37/') else
f'{basename.split("/")[0].replace("-", "_")}.tsv')
for basename in GTF_basenames]
for basename, version_mapping_file in zip(GTF_basenames,
version_mapping_files):
if not os.path.exists(version_mapping_file):
print(f'Generating "{version_mapping_file}"...')
run(f'curl -fsSL https://ftp.ensembl.org/pub/'
f'{basename}.gtf.gz | zcat | sed -nr \'{sed_command}\' | '
f'awk \'!seen[$1]++\' > {version_mapping_file}')
run(f'awk \'!seen[$1]++ {{print $1, $2, gensub(/.*\\/([^/]+)\\..*/, '
f'"\\\\1", "", FILENAME)}}\' OFS="\t" '
f'{" ".join(version_mapping_files)} > {mapping_file}')
return pl.read_csv(mapping_file, separator='\t', has_header=False,
new_columns=['Ensembl_ID', 'gene_name',
'most_recent_Ensembl_version'],
comment_prefix='#') \
.pipe(lambda df: df if no_unaliasing else unalias(
df, 'gene_name', gene_annotation_dir=gene_annotation_dir))
def check_valid_genome_build(genome_build):
"""
Checks whether a genome build is valid (for the functions here)
Args:
genome_build: the genome build; must be hg19 or hg38
"""
assert genome_build == 'hg19' or genome_build == 'hg38', genome_build
def get_coding_genes(genome_build='hg38', gencode_version=46,
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations',
return_file=False):
"""
Get the coordinates, gene symbols and Ensembl IDs of coding genes on
autosomes and sex chromosomes.
Args:
genome_build: the genome build (hg19 or hg38) to get coordinates for
gencode_version: the Gencode version to take coding genes from
gene_annotation_dir: the directory where a BED file of the coding genes
will be cached. Must be run on the login node to
generate this cache, if it doesn't exist. Because
generating the cache can take a long time, you
will probably want to leave this argument at its
default value.
return_file: If True, return a BED file path instead of a DataFrame.
Note: the start coordinates in the BED file are one less
than in the returned DataFrame, because BED files are
0-based while the returned DataFrame (and most sumstats)
are 1-based.
Returns:
A DataFrame with chrom, start, end, strand, gene, and Ensembl_IDs as
columns, or a BED file of the same if return_file=True. The gene names
are unaliased with unalias(), so you don't have to run unalias() again.
Genes in the pseudoautosomal regions have chrX as their chromosome, and
their chrY positions are not reported. Ensembl_IDs is a list column
since gene symbols very occasionally have multiple Ensembl IDs.
"""
check_valid_genome_build(genome_build)
coding_genes_file = f'{gene_annotation_dir}/coding_genes_{genome_build}_' \
f'gencode_v{gencode_version}.bed'
if os.path.exists(coding_genes_file):
return coding_genes_file if return_file else pl.read_csv(
coding_genes_file, separator='\t', has_header=False,
new_columns=['chrom', 'start', 'end', 'strand', 'gene',
'Ensembl_IDs']) \
.with_columns(pl.col.start + 1, # convert BED to one-based
pl.col.Ensembl_IDs.str.split(','))
os.makedirs(gene_annotation_dir, exist_ok=True)
coding_genes_intermediate_file = \
coding_genes_file.removesuffix('.bed') + '.intermediate.bed'
if not os.path.exists(coding_genes_intermediate_file):
raise_error_if_on_compute_node()
gencode_URL = (
f'https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/'
f'release_{gencode_version}'
f'{"/GRCh37_mapping" if genome_build == "hg19" else ""}/gencode.v'
f'{gencode_version}{"lift37" if genome_build == "hg19" else ""}.'
f'annotation.gtf.gz')
# Subtract 1 from the start coordinate because BED is 0-based, whereas
# GTF is 1-based
run(f'curl -fsSL {gencode_URL} | zcat | tr -d ";\\"" | awk \'$3 == '
f'"gene" && $0 ~ /protein_coding|IG_.*_gene|TR_.*_gene/ {{for '
f'(x = 1; x <= NF; x++) {{ if ($x == "gene_name") gene_name = '
f'$(x + 1); else if ($x == "gene_id") gene_id = $(x + 1); }} '
f'print $1, $4 - 1, $5, $7, gene_name, gene_id}}\' OFS="\t" | '
f'sort -k1,1V -k2,3n > {coding_genes_intermediate_file}')
coding_genes = pl.scan_csv(
coding_genes_intermediate_file, separator='\t', has_header=False,
new_columns=['chrom', 'start', 'end', 'strand', 'gene',
'Ensembl_IDs']) \
.filter(pl.col.chrom != 'chrM') \
.with_columns(Ensembl_IDs=pl.col.Ensembl_IDs.str.split_exact('.', 1)
.struct.field('field_0')) \
.collect()
# For hg19, 7 coding genes were subsequently merged into another gene
# (PRAMEF21 -> PRAMEF20, NBPF16 -> NBPF15, FOXD4L2 -> FOXD4L4,
# MRC1L1 -> MRC1, ANXA8L2 -> ANXA8L1, ASAH2C -> ASAH2B, CT45A4 -> CT45A3);
# to avoid the same gene being present in two different locations, do NOT
# unalias. Also, for both hg19 and hg38, ZNF475 is an alias of ZFP1, but
# also a gene in its own right; we can similarly ignore the alias.
assert len(coding_genes
.with_columns(unaliased_gene='gene')
.pipe(unalias, 'unaliased_gene',
gene_annotation_dir=gene_annotation_dir)
.filter(pl.col.unaliased_gene != pl.col.gene)) == \
(1 if genome_build == 'hg38' else 8)
# No genes should have more than one chromosome, except those in the
# pseudoautosomal regions
chrX_PAR1_end = 2_781_479 if genome_build == 'hg38' else 2_699_520
chrX_PAR2_start = 155_701_383 if genome_build == 'hg38' else 154_931_044
chrY_PAR1_end = 2_781_479 if genome_build == 'hg38' else 2_649_520
chrY_PAR2_start = 56_887_903 if genome_build == 'hg38' else 59_034_050
non_PAR_coding_genes = coding_genes.filter(
~(pl.col.chrom.eq('chrX') & (pl.col.end.lt(chrX_PAR1_end) |
pl.col.start.ge(chrX_PAR2_start))),
~(pl.col.chrom.eq('chrY') & (pl.col.end.lt(chrY_PAR1_end) |
pl.col.start.ge(chrY_PAR2_start))))
assert len(non_PAR_coding_genes.filter(
pl.col.chrom.n_unique().over('gene') != 1)) == 0
# No genes should have more than one strand
assert len(coding_genes.filter(
pl.col.strand.n_unique().over('gene') != 1)) == 0
# Occasionally, non-pseudoautosomal genes may appear more than once, always
# under different Ensembl IDs. Confirm that these genes are always
# overlapping, then merge them, and make a list of their Ensembl IDs.
duplicated_non_PAR_genes = \
non_PAR_coding_genes.filter(pl.col.gene.is_duplicated())
assert len(duplicated_non_PAR_genes.filter(
pl.col.start.max().over('gene') >= pl.col.end.min().over('gene'))) == 0
assert not duplicated_non_PAR_genes['Ensembl_IDs'].is_duplicated().any()
# Merge isoforms and genes with multiple Ensembl IDs into a single entry
# per gene symbol. The start and end are the min and max across isoforms,
# and the Ensembl IDs are aggregated into a list column. Only report the
# chrX position for pseudoautosomal genes by filtering out the chrY
# pseudoautosomal regions ahead of time.
coding_genes = coding_genes \
.lazy() \
.filter(~(pl.col.chrom.eq('chrY') & (
pl.col.end.lt(chrY_PAR1_end) | pl.col.start.ge(chrY_PAR2_start)))) \
.group_by('gene', maintain_order=True) \
.agg(pl.first('chrom'), pl.min('start'), pl.max('end'),
pl.first('strand'), 'Ensembl_IDs') \
.with_columns(pl.col.Ensembl_IDs.list.unique().list.sort()) \
.select('chrom', 'start', 'end', 'strand', 'gene', 'Ensembl_IDs') \
.collect()
assert not coding_genes['gene'].is_duplicated().any()
# Write to a file, and return
coding_genes \
.with_columns(pl.col.Ensembl_IDs.list.join(',')) \
.write_csv(coding_genes_file, separator='\t', include_header=False)
run(f"rm '{coding_genes_intermediate_file}'")
return coding_genes_file if return_file else \
coding_genes.with_columns(pl.col.start + 1)
def get_coding_exons(genome_build='hg38', gencode_version=46,
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations',
return_file=False):
"""
Get the coordinates, gene symbols and Ensembl IDs of exons of coding genes
on autosomes and sex chromosomes.
Args:
genome_build: the genome build (hg19 or hg38) to get coordinates for
gencode_version: the Gencode version to take coding exons from
gene_annotation_dir: the directory where a BED file of the coding exons
will be cached. Must be run on the login node to
generate this cache, if it doesn't exist. Because
generating the cache can take a long time, you
will probably want to leave this argument at its
default value.
return_file: If True, return a BED file path instead of a DataFrame.
Note: the start coordinates in the BED file are one less
than in the returned DataFrame, because BED files are
0-based while the returned DataFrame (and most sumstats)
are 1-based.
Returns:
A DataFrame with chrom, start, end, strand, gene, and Ensembl_IDs of
each exon as columns, or a BED file of the same if return_file=True.
The gene names are unaliased with unalias(), so you don't have to run
unalias() again. Exons in the pseudoautosomal regions have chrX as
their chromosome. Ensembl_IDs is a list column since gene symbols very
occasionally have multiple Ensembl IDs.
"""
check_valid_genome_build(genome_build)
coding_exons_file = f'{gene_annotation_dir}/coding_exons_{genome_build}_' \
f'gencode_v{gencode_version}.bed'
if os.path.exists(coding_exons_file):
return coding_exons_file if return_file else pl.read_csv(
coding_exons_file, separator='\t', has_header=False,
new_columns=['chrom', 'start', 'end', 'strand', 'gene',
'Ensembl_IDs']) \
.with_columns(pl.col.start + 1, # convert BED to one-based
pl.col.Ensembl_IDs.str.split(','))
os.makedirs(gene_annotation_dir, exist_ok=True)
coding_exons_intermediate_file = \
coding_exons_file.removesuffix('.bed') + '.intermediate.bed'
if not os.path.exists(coding_exons_intermediate_file):
raise_error_if_on_compute_node()
gencode_URL = (
f'https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/'
f'release_{gencode_version}'
f'{"/GRCh37_mapping" if genome_build == "hg19" else ""}/gencode.v'
f'{gencode_version}{"lift37" if genome_build == "hg19" else ""}.'
f'annotation.gtf.gz')
# Subtract 1 from the start coordinate because BED is 0-based, whereas
# GTF is 1-based
run(f'curl -fsSL {gencode_URL} | zcat | tr -d ";\\"" | awk \'$3 == '
f'"exon" && $0 ~ /protein_coding|IG_.*_gene|TR_.*_gene/ {{for '
f'(x = 1; x <= NF; x++) {{ if ($x == "gene_name") gene_name = '
f'$(x + 1); else if ($x == "gene_id") gene_id = $(x + 1); }} '
f'print $1, $4 - 1, $5, $7, gene_name, gene_id}}\' OFS="\t" | '
f'sort -k1,1V -k2,3n > {coding_exons_intermediate_file}')
coding_exons = pl.scan_csv(
coding_exons_intermediate_file, separator='\t', has_header=False,
new_columns=['chrom', 'start', 'end', 'strand', 'gene',
'Ensembl_IDs']) \
.filter(pl.col.chrom != 'chrM') \
.with_columns(Ensembl_IDs=pl.col.Ensembl_IDs.str.split_exact('.', 1)
.struct.field('field_0')) \
.unique(maintain_order=True) \
.collect()
# For hg19, 7 coding genes were subsequently merged into another gene
# (PRAMEF21 -> PRAMEF20, NBPF16 -> NBPF15, FOXD4L2 -> FOXD4L4,
# MRC1L1 -> MRC1, ANXA8L2 -> ANXA8L1, ASAH2C -> ASAH2B, CT45A4 -> CT45A3);
# to avoid the same gene being present in two different locations, do NOT
# unalias. Also, for both hg19 and hg38, ZNF475 is an alias of ZFP1, but
# also a gene in its own right; we can similarly ignore the alias.
assert len(coding_exons
.with_columns(unaliased_gene='gene')
.pipe(unalias, 'unaliased_gene',
gene_annotation_dir=gene_annotation_dir)
.filter(pl.col.unaliased_gene != pl.col.gene)
.unique('gene')) == (1 if genome_build == 'hg38' else 8)
# No exons should have more than one chromosome, except those in the
# pseudoautosomal regions
chrX_PAR1_end = 2_781_479 if genome_build == 'hg38' else 2_699_520
chrX_PAR2_start = 155_701_383 if genome_build == 'hg38' else 154_931_044
chrY_PAR1_end = 2_781_479 if genome_build == 'hg38' else 2_649_520
chrY_PAR2_start = 56_887_903 if genome_build == 'hg38' else 59_034_050
non_PAR_coding_exons = coding_exons.filter(
~(pl.col.chrom.eq('chrX') & (pl.col.end.lt(chrX_PAR1_end) |
pl.col.start.ge(chrX_PAR2_start))),
~(pl.col.chrom.eq('chrY') & (pl.col.end.lt(chrY_PAR1_end) |
pl.col.start.ge(chrY_PAR2_start))))
assert len(non_PAR_coding_exons.filter(
pl.col.chrom.n_unique().over('gene') != 1)) == 0
# No exons should have more than one strand
assert len(coding_exons.filter(
pl.col.strand.n_unique().over('gene') != 1)) == 0
# Occasionally, non-pseudoautosomal exons with the same gene symbol may
# appear multiple times under different Ensembl IDs. Merge them and make a
# list of their Ensembl IDs.
coding_exons = coding_exons \
.group_by('chrom', 'start', 'end', 'strand', 'gene',
maintain_order=True) \
.agg('Ensembl_IDs')
# Write to a file, and return
coding_exons \
.with_columns(pl.col.Ensembl_IDs.list.join(',')) \
.write_csv(coding_exons_file, separator='\t', include_header=False)
run(f"rm '{coding_exons_intermediate_file}'")
return coding_exons_file if return_file else \
coding_exons.with_columns(pl.col.start + 1)
def get_coding_TSSs(genome_build='hg38', gencode_version=46,
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations',
return_file=False):
"""
Get the coordinates, gene symbols and Ensembl IDs of transcription start
sites (TSSs) for coding genes on autosomes and sex chromosomes. Genes may
have multiple TSSs, one per transcript.
Args:
genome_build: the genome build (hg19 or hg38) to get coordinates for
gencode_version: the Gencode version to take coding genes from
gene_annotation_dir: the directory where a BED file of the coding TSSs
will be cached. Must be run on the login node to
generate this cache, if it doesn't exist. Because
generating the cache can take a long time, you
will probably want to leave this argument at its
default value.
return_file: If True, return a BED file path instead of a DataFrame.
The BED file will have start = bp - 1 and end = bp, since
BED files are 0-based while the returned DataFrame (and
most sumstats) are 1-based.
Returns:
A DataFrame with chrom, bp, strand, gene, and Ensembl_ID as columns,
or a BED file of the same if return_file=True. The gene names are
unaliased with unalias(), so you don't have to run unalias() again.
Genes in the pseudoautosomal regions have chrX as their chromosome.
"""
check_valid_genome_build(genome_build)
coding_TSSs_file = f'{gene_annotation_dir}/coding_TSSs_{genome_build}_' \
f'gencode_v{gencode_version}.bed'
if os.path.exists(coding_TSSs_file):
return coding_TSSs_file if return_file else pl.read_csv(
coding_TSSs_file, separator='\t', has_header=False,
new_columns=['chrom', 'start', 'end', 'strand', 'gene',
'Ensembl_ID']) \
.drop('start') \
.rename({'end': 'bp'})
os.makedirs(gene_annotation_dir, exist_ok=True)
coding_TSSs_intermediate_file = \
coding_TSSs_file.removesuffix('.bed') + '.intermediate.bed'
if not os.path.exists(coding_TSSs_intermediate_file):
raise_error_if_on_compute_node()
gencode_URL = (
f'https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/'
f'release_{gencode_version}'
f'{"/GRCh37_mapping" if genome_build == "hg19" else ""}/gencode.v'
f'{gencode_version}{"lift37" if genome_build == "hg19" else ""}.'
f'annotation.gtf.gz')
# Subtract 1 from the start coordinate because BED is 0-based, whereas
# GTF is 1-based; also uniquify (sort -u) since multiple transcripts of
# the same gene may share a TSS
run(f'curl -fsSL {gencode_URL} | zcat | tr -d ";\\"" | awk \'$3 == '
f'"transcript" && $0 ~ /protein_coding|IG_.*_gene|TR_.*_gene/ '
f'{{for (x = 1; x <= NF; x++) {{ if ($x == "gene_name") '
f'gene_name = $(x + 1); else if ($x == "gene_id") gene_id = '
f'$(x + 1); }} print $1, $7 == "+" ? $4 - 1 : $5 - 1, '
f'$7 == "+" ? $4 : $5, $7, gene_name, gene_id}}\' OFS="\t" | '
f'sort -k1,1V -k2,3n -u > {coding_TSSs_intermediate_file}')
coding_TSSs = pl.scan_csv(
coding_TSSs_intermediate_file, separator='\t', has_header=False,
new_columns=['chrom', 'start', 'end', 'strand', 'gene', 'Ensembl_ID']) \
.filter(pl.col.chrom != 'chrM') \
.with_columns(Ensembl_IDs=pl.col.Ensembl_IDs.str.split_exact('.', 1)
.struct.field('field_0')) \
.collect()
# For hg19, 7 coding genes were subsequently merged into another gene
# (PRAMEF21 -> PRAMEF20, NBPF16 -> NBPF15, FOXD4L2 -> FOXD4L4,
# MRC1L1 -> MRC1, ANXA8L2 -> ANXA8L1, ASAH2C -> ASAH2B, CT45A4 -> CT45A3);
# to avoid the same gene being present in two different locations, do NOT
# unalias. Also, for both hg19 and hg38, ZNF475 is an alias of ZFP1, but
# also a gene in its own right; we can similarly ignore the alias.
assert len(coding_TSSs
.with_columns(unaliased_gene='gene')
.pipe(unalias, 'unaliased_gene',
gene_annotation_dir=gene_annotation_dir)
.filter(pl.col.unaliased_gene != pl.col.gene)
.unique('gene')) == (1 if genome_build == 'hg38' else 8)
# No TSSs should have more than one chromosome, except those in the
# pseudoautosomal regions
chrX_PAR1_end = 2_781_479 if genome_build == 'hg38' else 2_699_520
chrX_PAR2_start = 155_701_383 if genome_build == 'hg38' else 154_931_044
chrY_PAR1_end = 2_781_479 if genome_build == 'hg38' else 2_649_520
chrY_PAR2_start = 56_887_903 if genome_build == 'hg38' else 59_034_050
non_PAR_coding_TSSs = coding_TSSs.filter(
~(pl.col.chrom.eq('chrX') & (pl.col.end.lt(chrX_PAR1_end) |
pl.col.start.ge(chrX_PAR2_start))),
~(pl.col.chrom.eq('chrY') & (pl.col.end.lt(chrY_PAR1_end) |
pl.col.start.ge(chrY_PAR2_start))))
assert len(non_PAR_coding_TSSs.filter(
pl.col.chrom.n_unique().over('gene') != 1)) == 0
# No TSSs should have more than one strand
assert len(coding_TSSs.filter(
pl.col.strand.n_unique().over('gene') != 1)) == 0
# Make sure no genes share a TSS
assert len(coding_TSSs.filter(
pl.struct('chrom', 'start', 'end').is_duplicated())) == 0
# Write to a file, and return
coding_TSSs \
.write_csv(coding_TSSs_file, separator='\t', include_header=False)
run(f"rm '{coding_TSSs_intermediate_file}'")
return coding_TSSs_file if return_file else \
coding_TSSs.drop('start').rename({'end': 'bp'})
# noinspection GrazieInspection
def get_nearest_gene(sumstats, genome_build, *, gencode_version=46,
SNP_col='SNP', chrom_col='CHROM', bp_col='BP',
gene_col='gene', gene_distance_col='gene_distance',
TSS_col='TSS', TSS_distance_col='TSS_distance', TSS=None,
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations',
signed=False, include_details=False):
"""
Gets the nearest coding gene(s) and/or TSS(s) for each variant in sumstats.
Given a DataFrame of sumstats with columns SNP_col, chrom_col and bp_col,
returns a DataFrame with those three columns plus (by default) four others,
described below: gene_col, gene_distance_col, TSS_col, and
TSS_distance_col. If include_details=True, include additional columns with
more details (see below).
Args:
sumstats: the summary statistics, as a polars DataFrame
genome_build: the genome build of sumstats
gencode_version: the Gencode version to take coding genes from
SNP_col: the name of the variant ID column in sumstats
chrom_col: the name of the chromosome column in sumstats
bp_col: the name of the base-pair column in sumstats
gene_col: the name of the column in the output DataFrame containing the
gene symbol(s) of the nearest gene(s) to each variant, as a
list
gene_distance_col: the name of the column in the output DataFrame
containing the distance from each variant to (each
of) the nearest gene(s) in gene_col, as a list. The
distance is 0 if the variant is inside the gene
(i.e. between the TSS and TES), and the distance to
whichever of the TSS and TES is closer, otherwise.
TSS_col: the name of the column in the output DataFrame containing the
gene symbol(s) of the gene(s) with the nearest TSS(s) to each
variant, as a list. For genes with multiple transcripts, only
the transcript(s) with the nearest TSS are counted.
TSS_distance_col: the name of the column in the output DataFrame
containing the distance from each variant to the
TSS(s) of (each of) the gene(s) in TSS_col, as a list
TSS: if None (the default), include both gene_col/gene_distance_col and
TSS_col/TSS_distance_col in the output. If False, include only
gene_col/gene_distance_col. If True, include only
TSS_col/TSS_distance_col.
gene_annotation_dir: the directory where the coding gene locations
returned by get_coding_genes() will be cached.
Must be run on the login node to generate this
cache, if it doesn't exist. Because generating the
cache can take a long time, you will probably want
to leave this argument at its default value.
signed: whether to report signed instead of unsigned distances
include_details: whether to include additional columns in the output:
For nearest genes (TSS=False or TSS=None):
- previous_gene_boundary: the base-pair position of the previous
gene boundary
- previous_genes: a list of the previous gene(s)
- next_gene_boundary: the base-pair position of the next gene
boundary
- next_genes: a list of the next gene(s)
- genes_inside: a list of gene(s) containing the variant
For nearest TSSs (TSS=True or TSS=None):
- previous_TSS: the base-pair position of the previous TSS
- previous_TSS_genes: a list of the previous TSS's gene(s)
- next_TSS: the base-pair position of the next TSS
- next_TSS_genes: a list of the next TSS's gene(s)
Returns:
A DataFrame with chrom_col, bp_col, gene_col/distance_col and/or
TSS_col/TSS_distance_col columns, and optionally the other columns
above if include_details=True.
"""
# Sanity-check inputs
check_valid_genome_build(genome_build)
if sumstats.is_empty():
raise ValueError(f'sumstats is empty!')
if SNP_col not in sumstats:
raise ValueError(f'{SNP_col!r} is not in sumstats; specify SNP_col')
if chrom_col not in sumstats:
raise ValueError(f'{chrom_col!r} is not in sumstats; specify '
f'chrom_col')
if bp_col not in sumstats:
raise ValueError(f'{bp_col!r} is not in sumstats; specify bp_col')
if TSS is not True and TSS is not False and TSS is not None:
raise ValueError(f'TSS must be True, False, or None')
include_nearest_gene = TSS is not True
include_nearest_TSS = TSS is not False
if include_nearest_gene and gene_col in sumstats:
raise ValueError(f'{gene_col!r} is already in sumstats; specify a '
f'different gene_col')
if include_nearest_gene and gene_distance_col in sumstats:
raise ValueError(f'{gene_distance_col!r} is already in sumstats; '
f'specify a different gene_distance_col')
if include_nearest_TSS and TSS_col in sumstats:
raise ValueError(f'{TSS_col!r} is already in sumstats; specify a '
f'different TSS_col')
if include_nearest_TSS and TSS_distance_col in sumstats:
raise ValueError(f'{TSS_distance_col!r} is already in sumstats; '
f'specify a different TSS_distance_col')
if include_details and TSS is True:
raise ValueError(f'include_details must be False when TSS=True')
# Standardize sumstats' chromosomes and sort
standardized_sorted_sumstats = sumstats \
.lazy() \
.select(SNP_col,
pl.col(chrom_col).pipe(standardize_chromosomes),
pl.col(chrom_col).alias('original_chrom'),
bp_col) \
.sort(chrom_col, bp_col)
if include_nearest_gene:
# Get coding genes
coding_genes = get_coding_genes(
genome_build=genome_build, gencode_version=gencode_version,
gene_annotation_dir=gene_annotation_dir)
# Get coding gene boundaries (starts and ends). Occasionally genes
# share the same boundary, so make a list of the genes that share each
# boundary. Almost all of these lists have just one gene, though.
gene_boundaries = pl.concat([coding_genes.select('chrom', key, 'gene')
.rename({'chrom': chrom_col, key: bp_col})
for key in ('start', 'end')]) \
.group_by(chrom_col, bp_col).agg('gene') \
.sort(chrom_col, bp_col)
# If genes are nested (e.g. GPR52 is entirely contained within
# RABGAP1L) or partially overlapping such that a gene boundary is
# contained within another gene, include the containing gene in the
# list as well. This is so that we can correctly infer which genes
# contain each variant below.
gene_boundaries = gene_boundaries \
.lazy() \
.drop('gene') \
.with_row_index() \
.join(gene_boundaries
.lazy()
.with_row_index()
.explode('gene')
.group_by('gene')
.agg(pl.col.index.min().alias('min_index'),
pl.col.index.max().alias('max_index'))
.with_columns(pl.int_ranges(
'min_index', pl.col.max_index + 1,
dtype=pl.UInt32).alias('index'))
.drop('min_index', 'max_index')
.explode('index')
.group_by('index')
.agg('gene'),
on='index') \
.sort('index') \
.drop('index') \
.collect()
# join_asof both forward, to the start of the next gene boundary, and
# backward, to the end of the previous gene boundary.
#
# If any of the next and previous gene(s) are the same, then those are
# the nearest genes, and since the variant is inside them, the distance
# is 0.
#
# Otherwise, get the distances to both the next and previous genes:
# when equal (very rare), take all the next and all the previous genes
# as nearest genes; otherwise, just take the ones on the closer side.
#
# Deduplicate and sort each variant's list of nearest genes at the end.
nearest_genes = standardized_sorted_sumstats \
.join_asof(gene_boundaries.lazy()
.rename({bp_col: 'previous_gene_boundary',
'gene': 'previous_genes'}),
left_on=bp_col, right_on='previous_gene_boundary',
by=chrom_col, strategy='backward',
check_sortedness=False) \
.join_asof(gene_boundaries.lazy()
.rename({bp_col: 'next_gene_boundary',
'gene': 'next_genes'}),
left_on=bp_col, right_on='next_gene_boundary',
by=chrom_col, strategy='forward',
check_sortedness=False) \
.with_columns(pl.col('previous_genes', 'next_genes')
.fill_null([])) \
.with_columns(distance_to_previous=pl.col(bp_col) -
pl.col.previous_gene_boundary,
distance_to_next=pl.col.next_gene_boundary -
pl.col(bp_col),
genes_inside=pl.col.previous_genes
.list.set_intersection(pl.col.next_genes)) \
.with_columns(pl.when(pl.col.genes_inside.list.len() > 0)
.then('genes_inside')
.when(pl.col.distance_to_previous ==
pl.col.distance_to_next)
.then(pl.concat_list('previous_genes', 'next_genes'))
.when((pl.col.distance_to_previous <
pl.col.distance_to_next) |
pl.col.distance_to_next.is_null())
.then('previous_genes')
.otherwise('next_genes')
.alias(gene_col),
pl.when(pl.col.genes_inside.list.len() > 0)
.then(0)
.when((pl.col.distance_to_previous <
pl.col.distance_to_next) |
pl.col.distance_to_next.is_null())
.then(-pl.col.distance_to_previous if signed else
'distance_to_previous')
.otherwise('distance_to_next')
.alias(gene_distance_col)) \
.with_columns(pl.col(gene_col).list.unique().list.sort()) \
.drop(chrom_col) \
.rename({'original_chrom': chrom_col}) \
.collect()
if include_details:
nearest_genes = nearest_genes \
.select(SNP_col, chrom_col, bp_col, gene_col,
gene_distance_col, 'previous_gene_boundary',
'previous_genes', 'next_gene_boundary', 'next_genes',
'genes_inside')
else:
nearest_genes = nearest_genes \
.select(SNP_col, chrom_col, bp_col, gene_col,
gene_distance_col)
if include_nearest_TSS:
# Get coding TSSs; note that since no genes share a TSS, we do not need
# to group_by(chrom_col, bp_col).agg('gene') as we did for
# gene_boundaries. Instead, just box each gene as a one-element list.
coding_TSSs = get_coding_TSSs(
genome_build=genome_build, gencode_version=gencode_version,
gene_annotation_dir=gene_annotation_dir) \
.with_columns(pl.col.gene.reshape((-1, 1))
.cast(pl.List(pl.String)))
# The logic for finding the nearest TSS is much simpler: just do two
# asof joins to find the nearest TSS in each direction, then take the
# closer of the two (or both, if equidistant)
nearest_TSSs = standardized_sorted_sumstats \
.join_asof(coding_TSSs.lazy()
.rename({'chrom': chrom_col, 'bp': 'previous_TSS',
'gene': 'previous_TSS_genes'}),
left_on=bp_col, right_on='previous_TSS', by=chrom_col,
strategy='backward', check_sortedness=False) \
.join_asof(coding_TSSs.lazy()
.rename({'chrom': chrom_col, 'bp': 'next_TSS',
'gene': 'next_TSS_genes'}),
left_on=bp_col, right_on='next_TSS', by=chrom_col,
strategy='forward', check_sortedness=False) \
.with_columns(pl.col('previous_TSS_genes', 'next_TSS_genes')
.fill_null([])) \
.with_columns(distance_to_previous=pl.col(bp_col) -
pl.col.previous_TSS,
distance_to_next=pl.col.next_TSS - pl.col(bp_col)) \
.with_columns(pl.when(pl.col.distance_to_previous ==
pl.col.distance_to_next)
.then(pl.concat_list('previous_TSS_genes',
'next_TSS_genes'))
.when((pl.col.distance_to_previous <
pl.col.distance_to_next) |
pl.col.distance_to_next.is_null())
.then('previous_TSS_genes')
.otherwise('next_TSS_genes')
.alias(TSS_col),
pl.when((pl.col.distance_to_previous <
pl.col.distance_to_next) |
pl.col.distance_to_next.is_null())
.then(-pl.col.distance_to_previous if signed else
'distance_to_previous')
.otherwise('distance_to_next')
.alias(TSS_distance_col)) \
.with_columns(pl.col(TSS_col).list.unique().list.sort()) \
.drop(chrom_col) \
.rename({'original_chrom': chrom_col}) \
.collect()
if include_details:
nearest_TSSs = nearest_TSSs \
.select(SNP_col, chrom_col, bp_col, TSS_col, TSS_distance_col,
'previous_TSS', 'previous_TSS_genes', 'next_TSS',
'next_TSS_genes')
else:
nearest_TSSs = nearest_TSSs \
.select(SNP_col, chrom_col, bp_col, TSS_col, TSS_distance_col)
if TSS is None:
# noinspection PyUnboundLocalVariable
return pl.concat([nearest_genes,
nearest_TSSs.drop(SNP_col, chrom_col, bp_col)],
how='horizontal')
elif TSS is False:
# noinspection PyUnboundLocalVariable
return nearest_genes
else:
# noinspection PyUnboundLocalVariable
return nearest_TSSs
def get_nearby_genes(sumstats, genome_build, *,
cache_directory='.', gencode_version=46,
max_distance=500_000, SNP_col='SNP', chrom_col='CHROM',
bp_col='BP', gene_col='gene',
gene_distance_col='gene_distance',
gene_annotation_dir=f'{get_base_data_directory()}/'
f'gene-annotations',
signed=False):
"""
Gets all coding genes within max_distance for each variant in sumstats, via
bedtools window.
Given a DataFrame of sumstats with columns SNP_col, chrom_col and bp_col,
returns a DataFrame with these three columns plus two others, described
below: gene_col and gene_distance_col.
Args:
sumstats: the summary statistics, as a polars DataFrame. Sumstats must
be sorted (this is not checked).
genome_build: the genome build of sumstats
cache_directory: the directory where the temporary files produced by
`bedtools window` will be stored
gencode_version: the Gencode version to take coding genes from
max_distance: the maximum distance from a variant to get nearby genes
SNP_col: the name of the variant ID column in sumstats
chrom_col: the name of the chromosome column in sumstats
bp_col: the name of the base-pair column in sumstats
gene_col: the name of the column in the output DataFrame containing the
gene symbol of each nearby gene, as a list
gene_distance_col: the name of the column in the output DataFrame
containing the distance from the variant to each of
the nearby genes in gene_col, as a list. The
distance is 0 if the variant is inside the gene
(i.e. between the TSS and TES), and the distance to
whichever of the TSS and TES is closer, otherwise.
gene_annotation_dir: the directory where the coding gene locations
returned by get_coding_genes() will be cached.
Must be run on the login node to generate this
cache, if it doesn't exist. Because generating the
cache can take a long time, you will probably want
to leave this argument at its default value.
signed: whether to report signed instead of unsigned distances
Returns:
A DataFrame with SNP_col, chrom_col, bp_col, gene_col,
gene_distance_col, and TSS_distance_col columns.
"""
# Sanity-check inputs
check_valid_genome_build(genome_build)
if sumstats.is_empty():
raise ValueError(f'sumstats is empty!')
if SNP_col not in sumstats:
raise ValueError(f'{SNP_col!r} not in sumstats; specify SNP_col')
if chrom_col not in sumstats:
raise ValueError(f'{chrom_col!r} not in sumstats; specify chrom_col')
if bp_col not in sumstats:
raise ValueError(f'{bp_col!r} not in sumstats; specify bp_col')
# Get coding genes BED file
coding_genes_bed_file = get_coding_genes(
genome_build=genome_build, gencode_version=gencode_version,
gene_annotation_dir=gene_annotation_dir, return_file=True)
# Make sumstats BED file
cache_prefix = f'{cache_directory}/get_nearby_genes'
sumstats_bed_file = f'{cache_prefix}.sumstats.bed'
sumstats \
.with_columns(pl.col(chrom_col).pipe(standardize_chromosomes),
original_chrom=pl.col(chrom_col)) \
.select(chrom_col, pl.col(bp_col).sub(1).alias('start'),
pl.col(bp_col).alias('end'), pl.col(SNP_col),
'original_chrom') \
.write_csv(sumstats_bed_file, separator='\t', include_header=False)
# Run bedtools window
nearby_genes_bed_file = f'{cache_prefix}.nearby_genes.bed'
run(f'bedtools window -a {sumstats_bed_file} -b {coding_genes_bed_file} '
f'-w {max_distance} > {nearby_genes_bed_file}')
# Load output; remember to add 1 to the gene start to convert to 1-based
nearby_genes = pl.scan_csv(
nearby_genes_bed_file, separator='\t', has_header=False,
new_columns=['__get_nearby_genes_variant_chrom',
'__get_nearby_genes_variant_start', bp_col, SNP_col,
chrom_col, '__get_nearby_genes_gene_chrom',
'__get_nearby_genes_gene_start',
'__get_nearby_genes_gene_end',
'__get_nearby_genes_gene_strand', gene_col,
'__get_nearby_genes_gene_Ensembl_ID'],
schema_overrides={chrom_col: sumstats[chrom_col].dtype}) \
.select(SNP_col, chrom_col, bp_col, gene_col,
pl.col.__get_nearby_genes_gene_start.add(1),
'__get_nearby_genes_gene_end') \
.with_columns(pl.when(pl.col(bp_col)
.is_between('__get_nearby_genes_gene_start',
'__get_nearby_genes_gene_end'))
.then(0)
.otherwise(
pl.when((pl.col(bp_col) -
pl.col.__get_nearby_genes_gene_start <
pl.col.__get_nearby_genes_gene_end -
pl.col(bp_col)) |
pl.col.__get_nearby_genes_gene_end.is_null())
.then(pl.col.__get_nearby_genes_gene_start -
pl.col(bp_col))
.otherwise(pl.col.__get_nearby_genes_gene_end -
pl.col(bp_col))
if signed else
pl.min_horizontal(
pl.col('__get_nearby_genes_gene_start',
'__get_nearby_genes_gene_end')
.sub(pl.col(bp_col)).abs()))
.alias(gene_distance_col)) \
.with_columns(pl.all()
.sort_by(pl.col(gene_distance_col).abs()
if signed else pl.col(gene_distance_col))
.over(SNP_col, chrom_col, bp_col)) \
.group_by(SNP_col, chrom_col, bp_col, maintain_order=True) \
.agg(pl.col(gene_col, gene_distance_col)) \
.collect()
# Align to sumstats, filling nulls with [] for variants with no genes
# within max_distance.
nearby_genes = sumstats \
.select(SNP_col, chrom_col, bp_col) \
.join(nearby_genes, on=(SNP_col, chrom_col, bp_col), how='left') \
.with_columns(pl.col(gene_col, gene_distance_col).fill_null([]))
# Clean up - but only at the end, so users can inspect what went wrong in
# case of an error
run(f"rm '{sumstats_bed_file}' '{nearby_genes_bed_file}'")
return nearby_genes
def get_nearest_variant(
genes, sumstats,
*,
cache_directory=os.environ['SCRATCH']
if os.environ.get('CLUSTER') == 'niagara' else '.',
gene_col='gene', gene_chrom_col='chrom', gene_start_col='start',
gene_end_col='end', SNP_col='SNP', sumstats_chrom_col='CHROM',
sumstats_bp_col='BP', distance_col='distance', signed=False):
"""
The reverse of `get_nearest_gene()`: for each gene in `genes`, gets the
distance to the nearest variant in `sumstats` via `bedtools closest -d`.
Based on Caleb Ji's `get_gene_distances` function.
Given a DataFrame of sumstats with columns `SNP_col`, `chrom_col` and
`bp_col`, and a DataFrame of genes with columns `gene_col`,
`gene_chrom_col`, `gene_start_col`, and `gene_end_col`, returns a DataFrame
with the four columns from the genes DataFrame plus two others: `SNP_col`,
the nearest variant, and `distance_col`, the distance from the gene to this
nearest variant.
Args:
genes: the list of genes, as a polars DataFrame. Must be sorted (this
is not checked). To use all protein-coding genes, set
`genes=get_coding_genes(genome_build)`.
sumstats: the summary statistics, as a polars DataFrame. Must be sorted
(this is not checked), and `genes` and `sumstats` must use
the same genome build (this is also not checked).
cache_directory: the directory where the temporary files produced by
`bedtools closest` will be stored
gene_col: the name of the gene symbol column in `genes`
gene_chrom_col: the name of the chromosome column in `genes`
gene_start_col: the name of the start column in `genes`
gene_end_col: the name of the end column in `genes`
SNP_col: the name of the variant ID column in `sumstats`.
sumstats_chrom_col: the name of the chromosome column in `sumstats`
sumstats_bp_col: the name of the base-pair column in `sumstats`
distance_col: the name of the column in the output DataFrame containing
the distance from each gene to its nearest variant. The
distance is 0 if the variant is inside the gene (i.e.
between the TSS and TES), and the distance to whichever
of the TSS and TES is closer, otherwise.
signed: whether to report signed instead of unsigned distances
Returns:
A DataFrame with `gene_col`, `gene_chrom_col`, `gene_start_col`,
`gene_end_col`, `SNP_col`, and `distance_col` columns. `SNP_col` and
`distance_col` will be `null` if the gene has no variants in `sumstats`
on the same chromosome.
"""
# Sanity-check inputs
if genes.is_empty():
error_message = 'genes is empty!'
raise ValueError(error_message)
if gene_col not in genes:
error_message = f'{gene_col!r} not in genes; specify gene_col'
raise ValueError(error_message)
if gene_chrom_col not in genes:
error_message = \
f'{gene_chrom_col!r} not in genes; specify gene_chrom_col'
raise ValueError(error_message)
if gene_start_col not in genes:
error_message = \
f'{gene_start_col!r} not in genes; specify gene_start_col'
raise ValueError(error_message)
if gene_end_col not in genes:
error_message = f'{gene_end_col!r} not in genes; specify gene_end_col'
raise ValueError(error_message)
if sumstats.is_empty():
error_message = 'sumstats is empty!'
raise ValueError(error_message)
if SNP_col not in sumstats:
error_message = f'{SNP_col!r} not in sumstats; specify SNP_col'
raise ValueError(error_message)
if sumstats_chrom_col not in sumstats:
error_message = (
f'{sumstats_chrom_col!r} not in sumstats; specify '
f'sumstats_chrom_col')
raise ValueError(error_message)
if sumstats_bp_col not in sumstats:
error_message = \
f'{sumstats_bp_col!r} not in sumstats; specify sumstats_bp_col'
raise ValueError(error_message)
# Make gene BED file
cache_prefix = f'{cache_directory}/get_nearest_variant'
gene_bed_file = f'{cache_prefix}.genes.bed'
genes \
.select(pl.col(gene_chrom_col).pipe(standardize_chromosomes),
gene_start_col, gene_end_col, gene_col) \
.write_csv(gene_bed_file, separator='\t', include_header=False)
# Make sumstats BED file
sumstats_bed_file = f'{cache_prefix}.sumstats.bed'
sumstats \
.with_columns(pl.col(sumstats_chrom_col).pipe(standardize_chromosomes),
original_chrom=pl.col(sumstats_chrom_col)) \
.select(sumstats_chrom_col,
pl.col(sumstats_bp_col).sub(1).alias('start'),
pl.col(sumstats_bp_col).alias('end'), pl.col(SNP_col),
'original_chrom') \
.write_csv(sumstats_bed_file, separator='\t', include_header=False)
# Run `bedtools closest`
distance_bed_file = f'{cache_prefix}.distance.bed'
run(f'bedtools closest {"-D ref" if signed else "-d"} -a {gene_bed_file} '
f'-b {sumstats_bed_file} > {distance_bed_file}')
nearest_variants = pl.scan_csv(
distance_bed_file, separator='\t', has_header=False,
new_columns=['__get_nearest_variant_gene_chrom', gene_start_col,
gene_end_col, '__get_nearest_variant_gene_strand',
gene_col, '__get_nearest_variant_gene_Ensembl_ID',
'__get_nearest_variant_variant_chrom',
'__get_nearest_variant_variant_start',
'__get_nearest_variant_variant_end', SNP_col,
gene_chrom_col, distance_col],
schema_overrides={gene_chrom_col: genes[gene_chrom_col].dtype}) \
.filter(~pl.col.distance.eq(-1)) \
.select([gene_col, gene_chrom_col, gene_start_col, gene_end_col,
SNP_col, distance_col]) \
.collect()
# Align with `genes`
nearest_variants = genes \
.select(gene_col, gene_chrom_col, gene_start_col, gene_end_col) \
.join(nearest_variants,
on=(gene_col, gene_chrom_col, gene_start_col, gene_end_col),
how='left')
# Clean up - but only at the end, so users can inspect what went wrong in
# case of an error
run(f"rm '{gene_bed_file}' '{sumstats_bed_file}' '{distance_bed_file}''")
return nearest_variants
def get_bim_or_pvar_file_type(bim_or_pvar_file):
"""
Gets whether a file is .bim or .pvar
Args:
bim_or_pvar_file: the bim or pvar file
Returns:
True if ends in .pvar, False if ends in .bim, raises an error otherwise
"""
is_pvar = bim_or_pvar_file.endswith('.pvar')
if not is_pvar and not bim_or_pvar_file.endswith('.bim'):
raise ValueError(f'bim_or_pvar_file "{bim_or_pvar_file}" must end '
f'with .bim or .pvar!')
return is_pvar
def read_bim_or_pvar(bim_or_pvar_file, categorical=False):
"""
Reads a plink variant file (.bim or .pvar) as a polars DataFrame.
Args:
bim_or_pvar_file: the .bim or .pvar file to read from; must end with
.bim or .pvar, and the extension determines which
file type it's parsed as
categorical: whether to read in the REF and ALT columns as Categorical
instead of String
Returns:
A polars DataFrame with the contents of the .bim or .pvar file,
with columns CHROM, SNP, CM (if present), BP, REF, ALT
"""
if not os.path.exists(bim_or_pvar_file):
raise FileNotFoundError(f'No such file or directory: '
f'{bim_or_pvar_file}')
is_pvar = get_bim_or_pvar_file_type(bim_or_pvar_file)
# If pvar, get the number of lines at the beginning starting with ##, and
# the number of lines after that starting with (should be 0 or 1)
if is_pvar:
num_double_comment_lines = int(run(
f"sed -n '/^##/!q;p' {bim_or_pvar_file} | wc -l",
stdout=subprocess.PIPE).stdout)
num_single_comment_lines = int(run(
f"sed -n '/^##/!{{/^#/p;q}}' {bim_or_pvar_file} | wc -l",
stdout=subprocess.PIPE).stdout)
if num_single_comment_lines > 1:
raise ValueError(f'pvar file "{bim_or_pvar_file}" has multiple '
f'header lines!')
skip_rows = num_double_comment_lines
has_header = bool(num_single_comment_lines)
else:
skip_rows = 0
has_header = False
dtypes = {'#CHROM' if has_header else 'column_1': pl.String}
if categorical:
if has_header:
ref_col = 'REF'
alt_col = 'ALT'
else:
first_line = pl.read_csv(bim_or_pvar_file, separator='\t',
has_header=has_header, n_rows=0)
if first_line.width not in (5, 6):
raise ValueError(f'bim_or_pvar_file "{bim_or_pvar_file}" has '
f'{first_line.width} columns, but should '
f'have 5 or 6!')
if first_line.width == 6:
ref_col = 'column_6'
alt_col = 'column_5'
else:
ref_col = 'column_5'
alt_col = 'column_4'
dtypes |= {ref_col: pl.Categorical('lexical'),
alt_col: pl.Categorical('lexical')}
variants = pl.read_csv(bim_or_pvar_file, separator='\t',
has_header=has_header, skip_rows=skip_rows,
schema_overrides=dtypes)
if has_header:
variants = variants.rename({'#CHROM': 'CHROM', 'POS': 'BP',
'ID': 'SNP'})
else:
if variants.width not in (5, 6):
raise ValueError(f'bim_or_pvar_file "{bim_or_pvar_file}" has '
f'{variants.width} columns, but should have 5 or '
f'6!')
if variants.width == 6:
variants = variants.rename({'column_1': 'CHROM', 'column_2': 'SNP',
'column_3': 'CM', 'column_4': 'BP',
'column_5': 'ALT', 'column_6': 'REF'})
else:
variants = variants.rename({'column_1': 'CHROM', 'column_2': 'SNP',
'column_3': 'BP', 'column_4': 'ALT',
'column_5': 'REF'})
return variants
def write_bim_or_pvar(variants, bim_or_pvar_file):
"""
Writes a polars DataFrame to bim or pvar format. If the filename ends in
.pvar, a header is added. If the centimorgan (CM) column isn't specified,
it's set to 0 (if bim) or omitted (if pvar).
Args:
variants: a polars DataFrame of variants with columns CHROM, SNP, CM,
BP, ALT, REF (CM is optional)
bim_or_pvar_file: the bim or pvar file to write to; must end with .bim
or .pvar, and the extension determines which file
type it's written as
"""
if variants.is_empty():
raise ValueError(f'variants is empty!')
is_pvar = get_bim_or_pvar_file_type(bim_or_pvar_file)
if is_pvar:
if 'CM' in variants:
columns = 'CHROM', 'BP', 'SNP', 'REF', 'ALT', 'CM'
else:
columns = 'CHROM', 'BP', 'SNP', 'REF', 'ALT'
else:
columns = 'CHROM', 'SNP', 'CM', 'BP', 'ALT', 'REF'
if 'CM' not in variants:
variants = variants.with_columns(CM=0)
variants = variants.select(columns) \
.rename({'CHROM': '#CHROM', 'BP': 'POS', 'SNP': 'ID'})
variants.write_csv(bim_or_pvar_file, separator='\t',
include_header=is_pvar)
def make_bim_or_pvar_IDs_unique(bim_or_pvar_file, new_bim_or_pvar_file,
separator='~'):
"""
Make the variant IDs of a bim or pvar file unique when there are
multiallelic variants, by replacing each SNP ID with separator.join(
(SNP, CHROM, REF, ALT)). Write the result to a new bim or pvar file.
On the off-chance a variant ID already contains the default separator, '~',
you will get an error and will need to specify a different string via the
separator argument.
Why include CHROM and not just ID/REF/ALT? Because the same variant may map
to multiple genomic locations (ncbi.nlm.nih.gov/snp/docs/rs_multi_mapping).
Why not include BP? Because it introduces a dependence on genome build. In
any case, it's vanishingly unlikely that a multi-mapping variant would have
the same base-pair position on two different chromosomes.
Args:
bim_or_pvar_file: the bim or pvar file to read from; must end with .bim
or .pvar, and the extension determines which file
type it's parsed and written as
new_bim_or_pvar_file: the bim or pvar file to write to after
uniquifying the variant IDs
separator: the separator to join the SNP, CHROM, REF, and ALT columns
with when uniquifying; no variant IDs may contain the
separator
"""
bim_or_pvar = read_bim_or_pvar(bim_or_pvar_file)
if bim_or_pvar['SNP'].str.contains(separator).any():
suffix = bim_or_pvar_file.split('.')[-1]
raise ValueError(f'Some variant IDs in the {suffix} file '
f'{bim_or_pvar_file} contain the separator '
f'{separator!r}; specify a different separator via '
f'the separator argument!')
bim_or_pvar \
.with_columns(pl.concat_str('SNP', 'CHROM', 'REF', 'ALT',
separator=separator)) \
.pipe(write_bim_or_pvar, new_bim_or_pvar_file)
def reverse_make_bim_or_pvar_IDs_unique(bim_or_pvar_file, new_bim_or_pvar_file,
separator='~'):
"""
Reverse the transformation of make_bim_or_pvar_IDs_unique() by removing the
CHROM, REF and ALT added to the SNP ID by that function. Write the result
to a new bim or pvar file.
Args:
bim_or_pvar_file: the bim or pvar file to read from; must end with .bim
or .pvar, and the extension determines which file
type it's parsed and written as
new_bim_or_pvar_file: the bim or pvar file to write to after removing
the CHROM, REF and ALT from the variant IDs
separator: the separator that the SNP, CHROM, REF, and ALT columns
were joined with when uniquifying
"""
bim_or_pvar = read_bim_or_pvar(bim_or_pvar_file)
if not bim_or_pvar['SNP'].str.contains(separator).all():
suffix = bim_or_pvar_file.split('.')[-1]
raise ValueError(f'not all variant IDs in the {suffix} file '
f'{bim_or_pvar_file} contain the separator '
f'{separator!r}; specify a different separator via '
f'the separator argument!')
bim_or_pvar \
.with_columns(SNP=pl.col.SNP.str.split_exact(separator, 1)
.struct.field('field_0')) \
.pipe(write_bim_or_pvar, new_bim_or_pvar_file)
def merge_bfiles(bfiles, merged_bfile):
"""
Merge multiple bed/bim/fam filesets (bfiles) into one big fileset
(merged_bfile).
As an optimization, just concatenates the raw bytes of the bed files rather
than calling plink. As a result, all filesets in bfiles must have the same
sample info (.fam file) and non-overlapping variants, though this is not
checked for speed. This is almost always true when merging, e.g. in the
very common case where you want to merge genetic data with one fileset per
chromosome into a single fileset.
Supports multiallelic variants, unlike plink 1.9's --merge-list and plink
2.0's --pmerge-list. (More precisely, plink 2.0 doesn't support
--pmerge-list on files with "split" multiallelic variants, but also doesn't
yet, as of September 2023, support merging multiallelic variants with
--make-pgen multiallelics=+. Even once this is supported, merge_bfiles()
should still be much faster than --pmerge-list.)
Args:
bfiles: a list/tuple/etc. of prefixes of bed/bim/fam filesets to merge
merged_bfile: the prefix of the merged bed/bim/fam fileset to be
created; {merged_bfile}.{bed,bim,fam} must not exist yet
"""
if isinstance(bfiles, str):
raise ValueError(f'bfiles must be a list/tuple/etc. of filesets, not '
f'a single string!')
if len(bfiles) < 2:
raise ValueError(f'bfiles must contain two or more filesets, not '
f'{len(bfiles):,}!')
for suffix in 'bed', 'bim', 'fam':
output_plink_file = f'{merged_bfile}.{suffix}'
if os.path.exists(output_plink_file):
raise FileExistsError(f'{output_plink_file} already exists!')
for bfile in bfiles:
input_plink_file = f'{bfile}.{suffix}'
if not os.path.exists(input_plink_file):
raise FileNotFoundError(f'No such file or directory: '
f'{input_plink_file}')
# Concatenate bed files: the first three bytes are the header, which is
# always 01101100 00011011 00000001, and the rest is the data. Take the
# header from the first file, and just the part after the header (tail
# -qc+4) for the remaining files.
run(f"(cat '{bfiles[0]}.bed'; tail -qc+4 " + ' '.join(
f"'{bfile}.bed'" for bfile in bfiles[1:]) +
f") > '{merged_bfile}.bed'")
# Concatenate bim files
run(f"cat " + ' '.join(f"'{bfile}.bim'" for bfile in bfiles) +
f" > '{merged_bfile}.bim'")
# Take the first fam file (since they're assumed to be all the same)
run(f"cp '{bfiles[0]}.fam' '{merged_bfile}.fam'")
def merge_pfiles(pfiles, merged_pfile, temp_dir=os.environ.get('SCRATCH',
'../../../../Desktop'),
num_threads=1, memory=2000):
"""
Merge multiple pgen/pvar/psam filesets (pfiles) into one big fileset
(merged_pfile) using plink 2.0's --pmerge-list.
--pmerge-list doesn't support multiallelic variants. (More precisely,
plink2 doesn't support --pmerge-list on files with "split" multiallelic
variants, but also doesn't yet, as of September 2023, support merging
multiallelic variants with --make-pgen multiallelics=+.) We circumvent this
by making the pvar IDs unique with make_bim_or_pvar_IDs_unique() before
merging, and reset them after with reverse_make_bim_or_pvar_IDs_unique().
Args:
pfiles: a list/tuple/etc. of prefixes of pgen/pvar/psam filesets to
merge
merged_pfile: the prefix of the merged pgen/pvar/psam fileset to be
created; {merged_bfile}.{pgen,pvar,psam} must not exist
yet
temp_dir: a directory to store the temporary uniquified pvars created
by make_bim_or_pvar_IDs_unique(), which will be deleted at
the end of the run if successful
num_threads: the number of threads to use for merging; if None, use all
available cores
memory: the number of megabytes of memory for plink to reserve during
merging via the --memory flag
"""
temp_pfiles = [os.path.join(temp_dir, f'{os.path.basename(pfile)}.tmp')
for pfile in pfiles]
for pfile, temp_pfile in zip(pfiles, temp_pfiles):
make_bim_or_pvar_IDs_unique(f'{pfile}.pvar', f'{temp_pfile}.pvar')
pmerge_list = ';'.join(f'echo {pfile}.pgen {temp_pfile}.pvar {pfile}.psam'
for pfile, temp_pfile in zip(pfiles, temp_pfiles))
run(f'plink2 '
f'{f"--threads {num_threads} " if num_threads is not None else ""}'
f'--memory {memory} '
f'--pmerge-list <({pmerge_list}) '
f'--make-pgen '
f'--out {merged_pfile}')
reverse_make_bim_or_pvar_IDs_unique(f'{merged_pfile}.pvar',
f'{merged_pfile}.pvar')
for temp_pfile in temp_pfiles:
run(f'rm "{temp_pfile}.pvar"')
run(f'rm "{merged_pfile}-merge.pgen" "{merged_pfile}-merge.pvar" '
f'"{merged_pfile}-merge.psam"')
def flip_alleles(sumstats, *, flip_col='FLIP', ref_col='REF', alt_col='ALT',
beta_col='BETA', OR_col='OR', Z_col='Z', AAF_col='AAF',
AAF_cases_col='AAF_CASES', AAF_controls_col='AAF_controls'):
"""
Flips alleles in sumstats where flip_col is True, and also flips effect
sizes, odds ratios, Z-score, and allele frequencies for these variants.
flip_col, ref_col, alt_col, and at least one of beta_col, OR_col and Z_col
are required. AAF_col, AAF_cases_col, AAF_controls_col, and the remainder
of beta_col, OR_col and Z_col are optional: if a column name doesn't appear
in sumstats, it won't be flipped, but there won't be an error.
Args:
sumstats: a polars DataFrame of summary statistics
flip_col: the name of a boolean column in sumstats saying which alleles
to flip
ref_col: the name of the reference allele column in sumstats
alt_col: the name of the alternate allele column in sumstats
beta_col: the name of the effect size column in sumstats
OR_col: the name of the odds ratio column in sumstats
Z_col: the name of the Z-score column in sumstats
AAF_col: the name of the alternate allele frequency column in sumstats
AAF_cases_col: the name of the case alternate allele frequency column
in sumstats
AAF_controls_col: the name of the control alternate allele frequency
column in sumstats
Returns:
sumstats with alleles flipped where flip_col is True.
"""
if ref_col not in sumstats:
raise ValueError(f'ref_col "{ref_col}" not in sumstats! Did you '
f'forget to specify a custom ref_col?')
if alt_col not in sumstats:
raise ValueError(f'alt_col "{alt_col}" not in sumstats! Did you '
f'forget to specify a custom alt_col?')
if beta_col not in sumstats and OR_col not in sumstats and \
Z_col not in sumstats:
raise ValueError(f'none of beta_col "{beta_col}", OR_col "{OR_col}", '
f'and Z_col "{Z_col}" are in sumstats! Did you '
f'forget to specify a custom beta_col, OR_col or '
f'Z_col?')
flip_transformations = {
ref_col: pl.col(alt_col), alt_col: pl.col(ref_col),
beta_col: -pl.col(beta_col), OR_col: 1 / pl.col(OR_col),
Z_col: -pl.col(Z_col), AAF_col: 1 - pl.col(AAF_col),
AAF_cases_col: 1 - pl.col(AAF_cases_col),
AAF_controls_col: 1 - pl.col(AAF_controls_col)}
return sumstats.with_columns(**{
column_name: pl.when(pl.col(flip_col)).then(transformation)
.otherwise(pl.col(column_name))
for column_name, transformation in flip_transformations.items()
if column_name in sumstats})
def harmonize_sumstats_to_bim_or_pvar(
sumstats, bim_or_pvar, *, SNP_col='SNP', ref_col='REF', alt_col='ALT',
beta_col='BETA', OR_col='OR', Z_col='Z', AAF_col='AAF',
AAF_cases_col='AAF_CASES', AAF_controls_col='AAF_controls',
chrom_col='CHROM', convert_chromosomes=False):
"""
Harmonizes sumstats to a bim or pvar file by flipping single nucleotide
variants' alleles with flip_alleles() as necessary to match the bim/pvar.
Indels that don't match, or single-nucleotide variants that don't match
even with flipping, are removed. If convert_chromosomes=True, converts
chromosomes to plink format ('1', '2', ..., '22', 'X', 'Y').
flip_col, ref_col, alt_col, and at least one of beta_col, OR_col and Z_col
are required. AAF_col, AAF_cases_col, AAF_controls_col, and the remainder
of beta_col, OR_col and Z_col are optional: if a column name doesn't appear
in sumstats, it won't be flipped, but there won't be an error.
Args:
sumstats: a polars DataFrame of sumstats
bim_or_pvar: a polars DataFrame of the contents of a bim or pvar file
created with get_rs_numbers_bim_or_pvar() and loaded into
memory with read_bim_or_pvar()
SNP_col: the name of the variant ID column in sumstats
ref_col: the name of the reference allele column in sumstats
alt_col: the name of the alternate allele column in sumstats
beta_col: the name of the effect size column in sumstats
OR_col: the name of the odds ratio column in sumstats
Z_col: the name of the Z-score column in sumstats
AAF_col: the name of the alternate allele frequency column in sumstats
AAF_cases_col: the name of the case alternate allele frequency column
in sumstats
AAF_controls_col: the name of the control alternate allele frequency
column in sumstats
chrom_col: the name of the chromosome column in sumstats; only used if
convert_chromosomes=True
convert_chromosomes: whether to convert the chromosomes in chrom_col to
plink format
Returns:
Sumstats with alleles flipped and chromosomes optionally converted to
plink format.
"""
# Use 0 as a dummy value just to check whether it's not null after joining.
return sumstats \
.lazy() \
.join(bim_or_pvar.lazy().select('SNP', 'REF', 'ALT',
matches_without_flips=0),
left_on=(SNP_col, ref_col, alt_col),
right_on=('SNP', 'REF', 'ALT'), how='left') \
.with_columns(pl.col.matches_without_flips.is_not_null()) \
.join(bim_or_pvar.lazy().select('SNP', 'REF', 'ALT',
matches_with_flips=0),
left_on=(SNP_col, ref_col, alt_col),
right_on=('SNP', 'ALT', 'REF'), how='left') \
.with_columns(pl.col.matches_with_flips.is_not_null() &
pl.col(ref_col).str.len_bytes().eq(1) &
pl.col(alt_col).str.len_bytes().eq(1)) \
.pipe(flip_alleles, flip_col='matches_with_flips', ref_col=ref_col,
alt_col=alt_col, beta_col=beta_col, OR_col=OR_col, Z_col=Z_col,
AAF_col=AAF_col, AAF_cases_col=AAF_cases_col,
AAF_controls_col=AAF_controls_col) \
.filter(pl.col.matches_without_flips | pl.col.matches_with_flips) \
.drop('matches_without_flips', 'matches_with_flips') \
.pipe(lambda df: df.with_columns(pl.col(chrom_col).pipe(
standardize_chromosomes, omit_chr_prefix=True))
if convert_chromosomes else df) \
.collect()
def ld_clump(sumstats, *,
cache_directory=os.environ['SCRATCH']
if os.environ.get('CLUSTER') == 'niagara'
else '.',
pfile=None, bfile=None,
clump_p1=5e-8, clump_p2=0.001, clump_r2=0.001, clump_kb=5000,
SNP_col='SNP', chrom_col='CHROM', bp_col='BP', ref_col='REF',
alt_col='ALT', p_col='P', separator='~', num_threads=1,
memory=2000, verbose=True):
"""
Performs linkage disequilibrium (LD) clumping on summary statistics using
plink's --clump (cog-genomics.org/plink/2.0/postproc#clump) and the LD info
from a plink pgen/pvar/psam ("pfile") or bed/bim/bam ("bfile") fileset.
Specify exactly one of pfile or bfile. pfile/bfile can be a list/tuple,
e.g. if there's 1 fileset per chromosome.
For example, try (on Niagara):
import os
import polars as pl
from utils import ld_clump
sumstats_file = '/scratch/w/wainberg/wainberg/sumstats/daner_MDDwoBP_' \
'20201001_2015iR15iex_HRC_MDDwoBP_iPSYCH2015i_' \
'UKBtransformed_Wray_FinnGen_MVPaf_2_HRC_MAF01.gz'
sumstats = pl.read_csv(sumstats_file, separator=' ', null_values='-')
clumped_variants = ld_clump(
sumstats, pfile='scratch/wainberg/1000G/European_autosomal_chrX_hg38',
chrom_col='CHR', ref_col='A2', alt_col='A1')
For information on LD clumping, see:
- cog-genomics.org/plink/1.9/postproc#clump
- zzz.bwh.harvard.edu/plink/clump.shtml
- explodecomputer.github.io/EEPE_2016/worksheets_win/bioinformatics.html
sumstats will be harmonized to pfile/bfile via
harmonize_sumstats_to_bim_or_pvar() and then written to the temp file
f'{cache_prefix}.sumstats.tmp' prior to clumping. This means that it's okay
if some ref and alt alleles are flipped between the sumstats file and the
pfile/bfile, and it does not matter which column is ref_col and which is
alt_col (though you may as well specify ref_col='REF', alt_col='ALT' or
ref_col='A2', alt_col='A1' for clarity).
ld_clump() ensures all variant IDs are unique by setting IDs for both
sumstats and pfile/bfile to ID + '~' + CHROM + '~' + REF + '~' + ALT. This
requires making a temporary pvar file (or files, if pfile/bfile is a list).
(On the off-chance a variant ID already contains '~', you will get an error
and will need to specify a different string via the separator argument.)
Why include CHROM and not just ID/REF/ALT? Because the same variant may map
to multiple genomic locations (ncbi.nlm.nih.gov/snp/docs/rs_multi_mapping).
Why not include BP? Because it would require the genome builds to match,
and it's vanishingly unlikely that a multi-mapping variant would have the
same base-pair position on two different chromosomes.
Args:
sumstats: the summary statistics, as a polars DataFrame
cache_directory: the directory where temporary sumstats, pvar files and
--clump results will be stored
pfile: the plink 2.x pgen/pvar/psam fileset LD info is taken from; can
be a string (assumed to be a single bfile for all chromosomes)
or a list/tuple of bfiles (assumed to be one per chromosome).
Exactly one of pfile and bfile must be specified.
bfile: the plink 1.x bed/bim/bam fileset LD info is taken from; can be
a string (assumed to be a single bfile for all chromosomes) or a
list/tuple of bfiles (assumed to be one per chromosome).
clump_p1: the value of --clump-p1 passed to --clump; for definitions of
--clump args, see cog-genomics.org/plink/2.0/postproc#clump
clump_p2: the value of --clump-p2 passed to --clump
clump_r2: the value of --clump-r2 passed to --clump
clump_kb: the value of --clump-kb passed to --clump
SNP_col: the name of the variant ID column in sumstats
chrom_col: the name of the chromosome column in sumstats
bp_col: the name of the base-pair column in sumstats; does NOT have to
be in the same genome build as pfile/bfile
ref_col: the name of the reference allele column in sumstats
alt_col: the name of the alternate column in sumstats
p_col: the name of the p-value column in sumstats
separator: the string to place between the variant ID, chromosome,
reference allele, and alternate allele when creating the
temporary sumstats and pvar files; no variant IDs may
contain the separator
num_threads: the number of threads to use for LD clumping; if None, use
all available cores
memory: the number of megabytes of memory for plink to reserve during
LD-clumping via the --memory flag
verbose: whether to print details of the LD clumping process
Returns:
A polars DataFrame with one row per variant in an LD clump, with the
following columns:
- SNP_col: the variant's ID
- ref_col: the variant's reference allele
- alt_col: the variant's alternate allele
- chrom_col: the variant's chromosome, in the same format as sumstats
- bp_col: the variant's base-pair position, in the same genome build as
sumstats
- f'{SNP_col}_lead': the ID of the lead variant, i.e. the lowest
p-value variant in the clump
- f'{ref_col}_lead': the lead variant's reference allele
- f'{alt_col}_lead': the lead variant's alternate allele
- f'{bp_col}_lead': the lead variant's base-pair position, in the same
genome build as sumstats
- 'is_lead': whether the variant is the lead variant at its locus
(equivalent to testing whether SNP_col ==
f'{SNP_col}_lead' and similarly for ref_col, alt_col and
bp_col)
- 'clump_start': the base-pair position of the start of the clump, in
the same genome build as sumstats
- 'clump_end': the base-pair position of the end of the clump, in the
same genome build as sumstats
"""
if sumstats.is_empty():
raise ValueError(f'sumstats is empty!')
if SNP_col not in sumstats:
raise ValueError(f'SNP_col "{SNP_col}" not in sumstats! Did you '
f'forget to specify a custom SNP_col?')
if chrom_col not in sumstats:
raise ValueError(f'chrom_col "{chrom_col}" not in sumstats! Did you '
f'forget to specify a custom chrom_col?')
if bp_col not in sumstats:
raise ValueError(f'bp_col "{bp_col}" not in sumstats! Did you '
f'forget to specify a custom bp_col?')
if ref_col not in sumstats:
raise ValueError(f'ref_col "{ref_col}" not in sumstats! Did you '
f'forget to specify a custom ref_col?')
if alt_col not in sumstats:
raise ValueError(f'alt_col "{alt_col}" not in sumstats! Did you '
f'forget to specify a custom alt_col?')
if p_col not in sumstats:
raise ValueError(f'p_col "{p_col}" not in sumstats!')
if sumstats[p_col].min() > clump_p1:
raise ValueError('All p-values in sumstats are > clump_p1')
if sumstats[SNP_col].str.contains(separator).any():
raise ValueError(f'Some variant IDs in sumstats contain the '
f'separator "{separator}"; specify a different '
f'separator via the separator argument!')
if pfile is None and bfile is None:
raise ValueError('Must specify either pfile or bfile (but not both)')
if pfile is not None and bfile is not None:
raise ValueError('Do not specify both pfile and bfile')
# Since only one of pfile and bfile is specified, refer to it with a single
# variable, filesets. If there's only one fileset, box it in a tuple
filesets = pfile if pfile is not None else bfile
if isinstance(filesets, str):
filesets = (filesets,)
suffix = 'pvar' if pfile is not None else 'bim'
# Load bim/pvar files from each fileset in pfile/bfile
bim_or_pvars = {fileset: read_bim_or_pvar(f'{fileset}.{suffix}')
for fileset in filesets}
for fileset, bim_or_pvar in bim_or_pvars.items():
if bim_or_pvar['SNP'].str.contains(separator).any():
raise ValueError(f'Some variant IDs in the {suffix} file '
f'{fileset}.{suffix} contain the separator '
f'{separator!r}; specify a different separator '
f'via the separator argument!')
# Harmonize sumstats to bim/pvar
old_length = len(sumstats)
sumstats = sumstats \
.with_columns(original_ref=ref_col, original_alt=alt_col) \
.pipe(harmonize_sumstats_to_bim_or_pvar,
pl.concat(bim_or_pvars.values()),
SNP_col=SNP_col, ref_col=ref_col, alt_col=alt_col,
chrom_col=chrom_col, convert_chromosomes=True)
num_filtered = old_length - len(sumstats)
if verbose:
print(f'Removing {num_filtered:,} {plural("variant", num_filtered)} '
f'in the sumstats '
f'({100 * num_filtered / old_length:.2f}%) '
f'that {"was" if num_filtered == 1 else "were"} not found in '
f'the {suffix} file')
# Create temporary sumstats file; map chromosomes to numbers to match plink
cache_prefix = f'{cache_directory}/ld_clump'
temp_sumstats_file = f'{cache_prefix}.sumstats.tmp'
sumstats \
.with_columns(pl.concat_str(SNP_col, chrom_col, ref_col, alt_col,
separator=separator)) \
.write_csv(temp_sumstats_file, separator='\t')
# Run --clump on each fileset. The default settings are:
# --clump-p1 5e-8: start with genome-wide-significant SNPs as lead SNPs
# --clump-p2 0.001: clumps will include all SNPs with p < 0.001...
# --clump-r2 0.01: ...r2 > 0.01 with the lead SNP...
# --clump-kb 5000: ...and within 5 MB of the lead SNP
# Create a temporary pvar file for each fileset first, with ID set to
# SNP + '_' + CHROM + '_' + REF + '_' + ALT.
for file_index, (fileset, bim_or_pvar) in enumerate(
bim_or_pvars.items(), start=1):
temp_pvar_file = f'{cache_prefix}.pvar' if len(filesets) == 1 else \
f'{cache_prefix}_{file_index}.pvar'
bim_or_pvar \
.with_columns(pl.concat_str('SNP', 'CHROM', 'REF', 'ALT',
separator=separator)) \
.pipe(write_bim_or_pvar, temp_pvar_file)
run(f'plink2 '
f'{f"--threads {num_threads} " if num_threads is not None else ""}'
f'--memory {memory} '
f'--seed 0 ' +
(f'--pfile {fileset.removesuffix(".pgen")} '
if pfile is not None else
f'--bfile {fileset.removesuffix(".bed")} ') +
f'--pvar {temp_pvar_file} '
f'--no-psam-pheno ' # optimization: avoids loading phenotypes
f'--clump {temp_sumstats_file} '
f'--clump-p1 {clump_p1} '
f'--clump-p2 {clump_p2} '
f'--clump-r2 {clump_r2} '
f'--clump-kb {clump_kb} '
f'--clump-snp-field {SNP_col} '
f'--clump-field {p_col} '
f'--out {cache_prefix}' +
('' if len(filesets) == 1 else f'_{file_index}'))
# Aggregate clumping results, if more than one fileset
clumping_results_file = f'{cache_prefix}.clumps'
if len(filesets) > 1:
run(f"(cat '{cache_prefix}_1.clumps'; tail -qn+2 " + ' '.join(
f"'{cache_prefix}_{file_index}.clumps'"
for file_index in range(2, len(filesets) + 1)) +
f") > '{clumping_results_file}'")
# Load clumping results; map each clumped variant to its lead variant
# (also remember to add a mapping from each lead variant to itself, and
# unescape the double-separator back to comma at the right moment)
if not os.path.exists(clumping_results_file):
raise RuntimeError(f'Clumping results file {clumping_results_file} is '
f'empty - this is a bug in ld_clump()!')
if os.path.exists(f'{clumping_results_file}.missing_id'):
raise RuntimeError(f'Some variants in sumstats were missing from '
f'{"pfile" if pfile is not None else "bfile"} - '
f'this is a bug in ld_clump()')
clumped_variants = pl.read_csv(clumping_results_file, separator='\t',
columns=['ID', 'SP2']) \
.with_columns(pl.when(pl.col.SP2 != '.').then(pl.col.SP2)
.str.split(',')) \
.explode('SP2') \
.pipe(lambda df: pl.concat([
df.filter(pl.col.SP2 != 'NONE')
.with_columns(SP2=pl.col.SP2.str.split_exact('(', 1)
.struct.field('field_0')),
# add a mapping from each lead variant to itself
pl.DataFrame({'ID': df['ID'].unique()}).with_columns(SP2='ID')])) \
.rename({'ID': 'lead_variant', 'SP2': SNP_col}) \
.select(SNP_col, 'lead_variant') \
.with_columns(pl.col(SNP_col).str.split_exact('~', 3).struct
.rename_fields([SNP_col, chrom_col, ref_col, alt_col]),
pl.col('lead_variant').str.split_exact('~', 3).struct
.rename_fields([f'{SNP_col}_lead', f'{chrom_col}_lead',
f'{ref_col}_lead', f'{alt_col}_lead'])) \
.unnest(SNP_col, 'lead_variant') \
.pipe(lambda df: df if sumstats[chrom_col].dtype == pl.String else df
.with_columns(pl.col(chrom_col, f'{chrom_col}_lead').cast(int))) \
.drop(f'{chrom_col}_lead')
assert not clumped_variants.select(SNP_col, chrom_col, ref_col, alt_col) \
.is_duplicated().any()
# Add the base-pair extent of each clump
clumped_variants = clumped_variants \
.join(sumstats.select(SNP_col, chrom_col, ref_col, alt_col, bp_col,
'original_ref', 'original_alt'),
on=(SNP_col, chrom_col, ref_col, alt_col), how='left') \
.join(sumstats.select(SNP_col, chrom_col, ref_col, alt_col, bp_col,
'original_ref', 'original_alt')
.rename({bp_col: f'{bp_col}_lead',
'original_ref': 'original_ref_lead',
'original_alt': 'original_alt_lead'}),
left_on=(f'{SNP_col}_lead', chrom_col, f'{ref_col}_lead',
f'{alt_col}_lead'),
right_on=(SNP_col, chrom_col, ref_col, alt_col), how='left') \
.with_columns(clump_start=pl.min(bp_col).over(f'{SNP_col}_lead'),
clump_end=pl.max(bp_col).over(f'{SNP_col}_lead'),
is_lead=pl.col(SNP_col).eq(pl.col(f'{SNP_col}_lead')) &
pl.col(ref_col).eq(pl.col(f'{ref_col}_lead')) &
pl.col(alt_col).eq(pl.col(f'{alt_col}_lead')) &
pl.col(bp_col).eq(pl.col(f'{bp_col}_lead'))) \
.with_columns(pl.col('original_ref').alias(ref_col),
pl.col('original_alt').alias(alt_col),
pl.col('original_ref_lead').alias(f'{ref_col}_lead'),
pl.col('original_alt_lead').alias(f'{alt_col}_lead')) \
.select(SNP_col, ref_col, alt_col, chrom_col, bp_col,
f'{SNP_col}_lead', f'{ref_col}_lead', f'{alt_col}_lead',
f'{bp_col}_lead', 'is_lead', 'clump_start', 'clump_end')
assert clumped_variants.null_count().sum_horizontal().item() == 0
# Sort by chromosome and base-pair position
clumped_variants = clumped_variants.pipe(
sort_sumstats, chrom_col=chrom_col, bp_col=bp_col)
# Clean up - but only at the end, so users can inspect what went wrong in
# case of an error
run(f"rm '{temp_sumstats_file}'")
if len(filesets) == 1:
run(f"rm '{cache_prefix}.pvar' '{cache_prefix}.clumps' "
f"'{cache_prefix}.log'")
else:
run(f"rm '{cache_prefix}.clumps' " +
' '.join(f"'{cache_prefix}_{file_index}.{suffix}'"
for file_index in range(1, len(filesets) + 1)
for suffix in ('pvar', 'clumps', 'log')))
return clumped_variants
def sort_sumstats(sumstats, *, chrom_col='CHROM', bp_col='BP'):
"""
Sorts summary statistics by chromosome and base-pair position.
Args:
sumstats: a polars DataFrame of summary statistics
chrom_col: the name of the chromosome column in sumstats
bp_col: the name of the base-pair position column in sumstats
Returns:
sumstats, sorted by chromosome and base-pair position.
"""
assert 'numeric_chrom' not in sumstats
return sumstats \
.with_columns(numeric_chrom=pl.col(chrom_col).pipe(
standardize_chromosomes, return_numeric=True)) \
.sort('numeric_chrom', bp_col) \
.drop('numeric_chrom')
# noinspection PyShadowingBuiltins
def munge_sumstats(raw_sumstats_files, munged_sumstats_file, REF, ALT, P, *,
SNP=None, CHROM=None, BP=None, BETA=None, OR=None, SE=None,
N=None, N_CASES=None, N_CONTROLS=None, NEFF=None, AAF=None,
AAF_CASES=None, AAF_CONTROLS=None, INFO=None,
quantitative=False, dbSNP=None, preamble=None,
separator='\t', verbose=True, **read_csv_kwargs):
"""
"Munges" the sumstats file(s) raw_sumstats_files into a consistent format,
saving to munged_sumstats_file (if not None) and returning the sumstats.
If dbSNP is not None, infers rs numbers. You should always do this if CHROM
and BP are available, even if the sumstats already came with rs numbers!
Make sure to set schema_overrides={'CHROM': str} explicitly when the
sumstats include sex chromosomes and there is no chr prefix (e.g. when it's
'1' and 'X' rather than 'chr1' and 'chrX').
Performs the following steps, in order, on each file in raw_sumstats_files:
1. Reads in the file with pl.read_csv(). Use separator to specify the
separator (tab by default, not comma!) and **kwargs to specify any other
arguments you want. Use separator='whitespace' to use any number of
consecutive whitespace characters as the separator.
2. Creates columns for each of the arguments from SNP to INFO. Specify a
single column name or polars expression of other columns, e.g.:
- AAF='allele_frequency'
- AAF=(pl.col.FCAS * pl.col.NCAS + pl.col.FCON * pl.col.NCON) /
(pl.col.NCAS + pl.col.NCON)
- AAF=pl.when(pl.col.A1 == pl.col.MINOR_ALLELE).then(pl.col.MAF)
.otherwise(1 - pl.col.MAF)
Certain columns are required: REF, ALT, P, CHROM + BP (if dbSNP is not
None), SNP (if dbSNP is None), and either BETA or OR (but not both). N
is required for quantitative traits (quantitative=True) and at least one
of NEFF, AAF, and N_CASES + N_CONTROLS is required for case-control
traits. Optional columns that are None won't be included in the output
unless they can be inferred from columns that are specified.
3. QC: removes variants with missing data in any column, non-ACGT alleles,
P outside (0, 1], OR <= 0, SE <= 0, N/N_CASES/N_CONTROLS/NEFF <= 0, AAF/
AAF_CASES/AAF_CONTROLS outside (0, 1), INFO outside (0, 1].
If verbose=True, prints stats on how many were removed.
4. Standardizes chromosomes to chr1, chr2, ... chr22, chrX, chrY
5. Converts variants to their minimal representations with
get_minimal_representations().
6. If dbSNP is not None, infers rs numbers based on CHROM and BP. If
verbose=True, prints the number and % of variants that had rs numbers in
dbSNP and the number and % that had to be flipped in order to match; if
SNP is not None as well, prints the number and % of rs numbers (variant
IDs starting with rs) in the SNP column that matched the ones from
dbSNP. These messages are useful for catching errors: for instance, you
may have specified a different genome build of dbSNP than the chrom/bp
positions in raw_sumstats_files, in which case most variants won't have
rs numbers, or you may have erroneously specified the alt columns as ref
and vice versa.
7. If dbSNP is not None, flips alleles for single-nucleotide variants if
necessary to match dbSNP, and sets BETA = -BETA, OR = 1 / OR, AAF = 1 -
AAF, AAF_CASES = 1 - AAF_CASES, and AAF_CONTROLS = 1 - AAF_CONTROLS for
these variants. Only REF/ALT flips are considered, NOT strand flips,
because modern GWAS datasets don't have strand flips - e.g. there isn't
a single non-ambiguous variant (i.e. not A/T, T/A, C/G, or G/C) in the
Als et al. 2023 depression GWAS or the Watson et al. 2019 anorexia GWAS
that matches dbSNP with a strand flip.
8. Removes variants without rs numbers in dbSNP and variants that aren't
unique (based on SNP + REF + ALT, + CHROM/BP if they're not None).
9. Sorts by CHROM and then BP, if those are not None.
10. Saves to munged_sumstats_file, unless munged_sumstats_file is None. If
munged_sumstats_file ends in .gz, sumstats will be block-gzipped with
bgzip for convenience.
11. Returns the sumstats.
For example, let's munge the depression sumstats on Niagara at
/scratch/w/wainberg/wainberg/sumstats/daner_MDDwoBP_20201001_2015iR15iex_
HRC_MDDwoBP_iPSYCH2015i_UKBtransformed_Wray_FinnGen_MVPaf_2_HRC_MAF01.gz.
The header line is:
CHR SNP BP A1 A2 FRQ_A_294322 FRQ_U_741438 INFO OR SE P ngt Direction \
HetISqt HetDf HetPVa Nca Nco Neff_half
So, filling in all the matching columns from left to right, we have:
munge_sumstats(
raw_sumstats_files='sumstats/daner_MDDwoBP_20201001_2015iR15iex_HRC_'
'MDDwoBP_iPSYCH2015i_UKBtransformed_Wray_FinnGen_'
'MVPaf_2_HRC_MAF01.gz',
munged_sumstats_file='MD.gz',
CHROM='CHR', SNP='SNP', BP='BP', ALT='A1', REF='A2', INFO='INFO',
OR='OR', SE='SE', P='P', N_CASES='Nca', N_CONTROLS='Nco', ...
To complete the picture, we must also specify AAF, which is a function of
the case and control Ns and frequencies. We can also specify NEFF, which
is two times the Neff_half column. We must also specify that the input
sumstats are space-delimited:
..., AAF=(pl.col.FRQ_A_294322 * pl.col.Nca + pl.col.FRQ_U_741438 *
pl.col.Nco) / (pl.col.Nca + pl.col.Nco)',
NEFF=2 * pl.col.Neff_half, separator=' ')
Args:
raw_sumstats_files: the filename of the input sumstats, or a
list/tuple/etc. of filenames to be concatenated
munged_sumstats_file: if not None, the output sumstats filename to save
to
REF: the reference allele column/expression; mandatory
ALT: the alternate allele column/expression; mandatory
P: the p-value column/expression; mandatory
SNP: The variant ID column in raw_sumstats_files, or a polars
expression to generate a variant ID column; mandatory when dbSNP
is None.
CHROM: the chromosome column/expression; mandatory when dbSNP is not
None or BP is not None
BP: the base-pair position column/expression; mandatory when dbSNP is
not None or CHROM is not None
BETA: The effect size column/expression. Specify exactly one of BETA
and OR; if OR is specified, BETA will be calculated as log(OR).
OR: The odds ratio column/expression. Specify exactly one of BETA/OR.
Not explicitly disallowed for quantitative traits!
SE: The standard error column/expression. If None, will be back-
calculated from BETA and P: SE = |BETA| / p_to_abs_z(P).
p_to_abs_z() is an awk implementation of utils.py's p_to_abs_z()
function. log(OR) is substituted for BETA if BETA is None.
N: The sample size column/expression. Mandatory for quantitative traits
(quantitative=True). For case-control traits (quantitative=False),
will be calculated as N_CASES + N_CONTROLS if both N_CASES and
N_CONTROLS are not None. Can be a number, in which case N is
assumed to be that number for all variants.
N_CASES: The number of cases column/expression. Must be None for
quantitative traits (quantitative=True). Can be a number, in
which case N_CASES is assumed to be the same for all variants.
N_CONTROLS: The number of controls column/expression. Must be None for
quantitative traits (quantitative=True). Can be a number,
in which case N_CONTROLS is assumed to be the same for all
variants.
NEFF: The effective sample size column/expression. For quantitative
traits (quantitative=True), will be set to N if missing. For
case-control traits (quantitative=False), will be set to
((4 / (2 * AAF * (1 - AAF) * INFO)) - BETA^2) / SE^2 if AAF is
not None (substituting log(OR) for BETA if BETA is None, and
skipping the "* INFO" if INFO is None), and
4 / (1 / N_CASES + 1 / N_CONTROLS) if both N_CASES and N_CONTROLS
are not None; see comment below for justification. As a result,
for case-control traits, either NEFF or AAF or both N_CASES and
N_CONTROLS are mandatory. Can be a number, in which case NEFF is
assumed to be the same for all variants.
AAF: the alternate allele frequency column/expression
AAF_CASES: the case alternate allele frequency column/expression
AAF_CONTROLS: the control alternate allele frequency column/expression
INFO: the imputation INFO score column/expression
quantitative: is the trait quantitative (True) or case-control (False)?
dbSNP: a DataFrame returned by load_dbSNP(); must match the genome
build of raw_sumstats_files' BP column! If not specified, do not
infer rs numbers.
preamble: an optional function to apply immediately after loading each
file in raw_sumstats_files (so use the original column names
in the expression, not the new column names). For instance,
to filter to variants with minor allele frequencies between
0.01 and 0.99, use preamble=lambda df: df.filter(
pl.col.all_maf.between(0.01, 0.99)).
separator: the separator to use when parsing each raw sumstats file
with pl.read_csv(), or 'whitespace' to use any number of
consecutive whitespace characters as the separator, like
delim_whitespace=True in pandas.read_table().
verbose: whether to print details of the munging process
**read_csv_kwargs: keyword arguments to pl.read_csv(). By default,
null_values='NA' is set.
Returns:
A DataFrame of the munged sumstats.
"""
# Check inputs
if REF is None:
raise ValueError('REF is always mandatory')
if ALT is None:
raise ValueError('ALT is always mandatory')
if P is None:
raise ValueError('P is always mandatory')
if dbSNP is None and SNP is None:
raise ValueError('SNP is mandatory when dbSNP is None')
if CHROM is None and BP is not None:
raise ValueError('CHROM is mandatory when BP is not None')
if CHROM is not None and BP is None:
raise ValueError('BP is mandatory when CHROM is not None')
if dbSNP is not None and BP is None:
raise ValueError('CHROM and BP are mandatory when dbSNP is not None')
if BETA is None and OR is None:
raise ValueError('Neither BETA nor OR specified; specify exactly one')
if BETA is not None and OR is not None:
raise ValueError('Both BETA and OR specified; specify exactly one')
if quantitative and N is None:
raise ValueError('N is mandatory when quantitative=True')
if quantitative and N_CASES is not None:
raise ValueError('N_CASES cannot be specified when quantitative=True')
if quantitative and N_CONTROLS is not None:
raise ValueError('N_CONTROLS cannot be specified when '
'quantitative=True')
if not quantitative and NEFF is None and AAF is None and \
(N_CASES is None or N_CONTROLS is None):
raise ValueError('quantitative=True and NEFF is None, so either AAF '
'or both of N_CASES and N_CONTROLS need to be '
'specified to infer it')
# Wrap string arguments in pl.col()
if isinstance(SNP, str):
SNP = pl.col(SNP)
if isinstance(REF, str):
REF = pl.col(REF)
if isinstance(ALT, str):
ALT = pl.col(ALT)
if isinstance(P, str):
P = pl.col(P)
if isinstance(CHROM, str):
CHROM = pl.col(CHROM)
if isinstance(BP, str):
BP = pl.col(BP)
if isinstance(BETA, str):
BETA = pl.col(BETA)
if isinstance(OR, str):
OR = pl.col(OR)
if isinstance(SE, str):
SE = pl.col(SE)
if isinstance(N, str):
N = pl.col(N)
if isinstance(N_CASES, str):
N_CASES = pl.col(N_CASES)
if isinstance(N_CONTROLS, str):
N_CONTROLS = pl.col(N_CONTROLS)
if isinstance(NEFF, str):
NEFF = pl.col(NEFF)
if isinstance(AAF, str):
AAF = pl.col(AAF)
if isinstance(AAF_CASES, str):
AAF_CASES = pl.col(AAF_CASES)
if isinstance(AAF_CONTROLS, str):
AAF_CONTROLS = pl.col(AAF_CONTROLS)
if isinstance(INFO, str):
INFO = pl.col(INFO)
# If SE is missing, back-calculate it from BETA and P: |Z| = |BETA| / SE,
# so SE = |BETA| / |Z|. We can get |Z| from P via p_to_abs_z().
if SE is None:
# noinspection PyUnresolvedReferences
SE = (BETA if BETA is not None else OR.log()).abs() / p_to_abs_z(P)
# Try to infer N and NEFF if not specified
#
# For quantitative traits, NEFF = N. For binary traits, NEFF can be given
# in two mathematically equivalent ways: 4/(1/N_CASES + 1/N_CONTROLS)
# (cell.com/ajhg/pdfExtended/S0002-9297(21)00145-2) and 4v(1-v)N where v =
# N_CASES/N (medrxiv.org/content/10.1101/2021.09.22.21263909v1.full-text).
#
# But this doesn't account for varying case-control ratios across the
# cohorts in a meta-analysis, which biases estimates that depend on N
# (medrxiv.org/content/10.1101/2021.09.22.21263909v1.full).
#
# Instead, biorxiv.org/content/10.1101/2021.03.29.437510v4.full suggests
# a per-variant Neff = (4 / (2 * MAF * (1 - MAF) * INFO) - BETA^2) / SE^2,
# bounded to between 0.5 and 1.1 times "the total (effective) sample size",
# which seems to be max(4/(1/N_CASES + 1/N_CONTROLS)) across SNPs:
# github.com/privefl/paper-misspec/blob/main/code/investigate-misspec-N.R
# To get a global Neff, they take the 80th %ile of the per-variant Neffs.
# This formula is based on Equation 1 of the "New formula used in LDpred2"
# section of sciencedirect.com/science/article/abs/pii/S0002929721004201.
#
# github.com/GenomicSEM/GenomicSEM/wiki/2.1-Calculating-Sum-of-Effective-
# Sample-Size-and-Preparing-GWAS-Summary-Statistics drops the BETA^2 and
# INFO, i.e. Neff = 4 / (2 * MAF * (1 - MAF)) / SE^2, bounded to between
# 0.5 and 1.1 times 4/(1/N_CASES + 1/N_CONTROLS), where N_CASES/N_CONTROLS
# are again global rather than per-SNP. This formula is based on an old
# version of the above formula that was used in the original LDpred2 paper.
#
# Here, we use Neff = (4 / (2 * AAF * (1 - AAF) * INFO) - BETA^2) / SE^2
# (dropping the INFO part if not available), without bounding (because it
# seems like a hack). If MAF not available, but N_CASES and N_CONTROLS are
# available, fall back to using 4/(1/N_CASES + 1/N_CONTROLS). Note that
# even if AAF = 1 - MAF instead of MAF, AAF * (1 - AAF) == MAF * (1 - MAF).
# "((4 / (2 * ..." has been simplified to "((2 / (..." below.
if N is None and N_CASES is not None and N_CONTROLS is not None:
N = N_CASES + N_CONTROLS
if NEFF is None:
if quantitative and N is not None:
NEFF = N
elif AAF is not None:
if INFO is not None:
NEFF = (2 / (AAF * (1 - AAF) * INFO) - (
BETA if BETA is not None else OR.log()) ** 2) / SE ** 2
else:
NEFF = (2 / (AAF * (1 - AAF)) - (
BETA if BETA is not None else OR.log()) ** 2) / SE ** 2
elif N_CASES is not None and N_CONTROLS is not None:
NEFF = 4 / (1 / N_CASES + 1 / N_CONTROLS)
else:
raise ValueError('NEFF is unexpectedly None; internal error')
# For case-control traits, output ORs even if input has betas, & vice versa
if not quantitative and OR is None:
OR = BETA.exp()
if quantitative and BETA is None:
BETA = OR.log()
# Ensure REF and ALT are capitalized
REF = REF.str.to_uppercase()
ALT = ALT.str.to_uppercase()
# Enumerate columns to include, and the formula for each column
column_formulas = {
column_name: formula for column_name, formula in {
'SNP': SNP, 'CHROM': CHROM, 'BP': BP, 'REF': REF, 'ALT': ALT,
'AAF': AAF, 'AAF_CASES': AAF_CASES, 'AAF_CONTROLS': AAF_CONTROLS,
'INFO': INFO, 'BETA' if quantitative else 'OR':
BETA if quantitative else OR, 'SE': SE, 'P': P, 'N': N,
'N_CASES': N_CASES, 'N_CONTROLS': N_CONTROLS, 'NEFF': NEFF}.items()
if formula is not None}
# Load files in raw_sumstats_files; apply the preamble function (if not
# None) and select columns
raw_sumstats_files = (raw_sumstats_files,) \
if isinstance(raw_sumstats_files, str) else tuple(raw_sumstats_files)
columns = tuple(set.union(*(set(formula.meta.root_names())
if isinstance(formula, pl.Expr) else formula
for formula in column_formulas.values()
if not isinstance(formula, (int, float)))))
default_read_csv_kwargs = dict(null_values='NA')
read_csv_kwargs = default_read_csv_kwargs | read_csv_kwargs \
if read_csv_kwargs is not None else default_read_csv_kwargs
sumstats = pl.concat((
read_csv_delim_whitespace(raw_sumstats_file,
columns=columns,
**read_csv_kwargs)
if separator == 'whitespace' else
# Use scan_csv() instead of read_csv() when possible, i.e. when the raw
# sumstats file isn't gzipped and delim_whitespace=False.
pl.read_csv(raw_sumstats_file, columns=columns,
separator=separator,
**read_csv_kwargs)
if raw_sumstats_file.endswith('.gz') else
pl.scan_csv(raw_sumstats_file,
separator=separator,
**read_csv_kwargs)
.select(columns))
.lazy()
.pipe(
preamble if preamble is not None else lambda df: df)
.select(**column_formulas)
for raw_sumstats_file in raw_sumstats_files) \
.collect()
# QC: remove variants with missing data in any column, non-ACGT alleles, P
# outside (0, 1], OR <= 0, SE <= 0, N/N_CASES/N_CONTROLS/NEFF <= 0, AAF/
# AAF_CASES/AAF_CONTROLS outside (0, 1), INFO outside (0, 1].
# If verbose=True, print stats on how many were removed.
if verbose:
num_initial_variants = len(sumstats)
# noinspection PyUnresolvedReferences
filters = {filter_name: filter for filter_name, filter in ({
'non-ACGT alleles':
sumstats[
'REF'].str.contains(
'[^ACGT]') |
sumstats[
'ALT'].str.contains(
'[^ACGT]'),
'P outside (0, 1]': ~
sumstats[
'P'].is_between(
0, 1,
closed='right'),
'OR <= 0':
sumstats[
'OR'] <= 0 if 'OR' in sumstats else None,
'SE <= 0':
sumstats[
'SE'] <= 0 if 'SE' in sumstats else None,
'N <= 0':
sumstats[
'N'] <= 0 if 'N' in sumstats else None,
'N_CASES <= 0':
sumstats[
'N_CASES'] <= 0
if 'N_CASES' in sumstats else None,
'N_CONTROLS <= 0':
sumstats[
'N_CONTROLS'] <= 0
if 'N_CONTROLS' in sumstats else None,
'NEFF <= 0':
sumstats[
'NEFF'] <= 0 if 'NEFF' in sumstats else None,
'AAF outside (0, 1)': ~
sumstats[
'AAF'].is_between(
0, 1,
closed='none')
if 'AAF' in sumstats else None,
'AAF_CASES outside (0, 1)': ~
sumstats[
'AAF_CASES'].is_between(
0, 1,
closed='none') if 'AAF_CASES' in sumstats else None,
'AAF_CONTROLS outside (0, 1)': ~
sumstats[
'AAF_CONTROLS'].is_between(
0, 1,
closed='none') if 'AAF_CONTROLS' in sumstats else None,
'INFO <= 0':
sumstats[
'INFO'] <= 0 if 'INFO' in sumstats else None} | {
f'null {column_name}':
sumstats[
column_name].is_null()
for
column_name
in
column_formulas}).items()
if filter is not None and filter.any()}
if len(filters) > 0:
union_of_filters = reduce(lambda a, b: a | b, filters.values())
total_num_filtered = union_of_filters.sum()
if total_num_filtered == len(sumstats):
raise ValueError('All variants would be filtered out!')
if total_num_filtered > 0:
if verbose:
last_filter_name = tuple(filters)[-1] \
if len(filters) > 1 else None
print(f'Removing ' + ', '.join(
f'{"and " if filter_name == last_filter_name else ""}'
f'{num_filtered:,} {plural("variant", num_filtered)} with '
f'{filter_name}'
for filter_name, filter, num_filtered in
((filter_name, filter, filter.sum())
for filter_name, filter in filters.items())) +
f', for a total of {total_num_filtered:,} '
f'{plural("variant", total_num_filtered)} '
f'({100 * total_num_filtered / len(sumstats):.2f}%)')
sumstats = sumstats.filter(~union_of_filters)
elif verbose:
print('All variants pass initial QC filters!')
# Standardize chromosomes
if CHROM is not None:
sumstats = sumstats \
.with_columns(pl.col.CHROM.pipe(standardize_chromosomes))
# Convert variants to their minimal representations
sumstats = sumstats.pipe(get_minimal_representations,
bp_col='BP' if BP is not None else None)
# If dbSNP is not None, infer rs numbers and flip variants to match dbSNP;
# if verbose=True, print stats
if dbSNP is not None:
if verbose and SNP is not None:
sumstats = sumstats.rename({'SNP': 'original_SNP'})
sumstats = sumstats \
.pipe(get_rs_numbers, dbSNP=dbSNP, verbose=verbose) \
.select('SNP', pl.exclude('SNP')) \
.pipe(flip_alleles)
if verbose:
num_nonmissing_rs = sumstats['SNP'].is_not_null().sum()
num_flipped = sumstats['FLIP'].sum()
percent_flipped = 100 * num_flipped / num_nonmissing_rs
print(f'{num_flipped:,} of these {num_nonmissing_rs:,} variants '
f'({percent_flipped:.2f}%) had to be flipped in order to '
f'match dbSNP')
if SNP is not None:
matching_rs = sumstats \
.select('original_SNP', 'SNP') \
.filter(pl.col.original_SNP.str.starts_with('rs'),
pl.col.SNP.is_not_null()) \
.with_columns(pl.col.SNP.str.split(',')) \
.with_columns(match=pl.col.SNP.list.contains(
pl.col.original_SNP))
num_matching_rs = matching_rs['match'].sum()
percent_matching_rs = 100 * num_matching_rs / len(matching_rs)
likely_merges = matching_rs \
.filter(~pl.col.match, ~pl.col.SNP.list.len().eq(1)) \
.select(pl.col.SNP.list.get(0).str.slice(2).cast(int) <
pl.col.original_SNP.str.slice(2).cast(int)) \
.to_series()
num_likely_merges = likely_merges.sum()
percent_likely_merges = 100 * num_likely_merges / \
len(likely_merges)
print(f'{len(matching_rs):,} variant IDs in the SNP column/'
f'expression you specified start with "rs" and have at '
f'least one rs number in dbSNP; {num_matching_rs:,} '
f'({percent_matching_rs:.2f}%) of these match dbSNP. Of '
f'those that didn\'t, {len(likely_merges):,} have '
f'exactly one matching rs number, and for '
f'{num_likely_merges:,} ({percent_likely_merges:.2f}%) '
f'of these, the matching rs number is numerically '
f'smaller than the original, which suggests the two rs '
f'numbers might have been merged.')
sumstats = sumstats.drop('original_SNP')
sumstats = sumstats.drop('FLIP')
# Remove variants without rs numbers in dbSNP
variants_with_rs_numbers = sumstats['SNP'].is_not_null()
num_without_rs_numbers = len(sumstats) - variants_with_rs_numbers.sum()
if num_without_rs_numbers > 0:
if verbose:
print(f'Removing {num_without_rs_numbers:,} '
f'{plural("variant", num_without_rs_numbers)} '
f'({100 * num_without_rs_numbers / len(sumstats):.2f}%) '
f'without rs numbers in dbSNP')
sumstats = sumstats.filter(variants_with_rs_numbers)
# Remove non-unique variants
variant_columns = ['SNP', 'REF', 'ALT'] + \
(['CHROM', 'BP'] if BP is not None else [])
unique_mask = sumstats.select(variant_columns).is_unique()
num_non_unique = len(sumstats) - unique_mask.sum()
if num_non_unique > 0:
if verbose:
print(f'Removing {num_non_unique:,} non-unique '
f'{plural("variant", num_non_unique)} '
f'({100 * num_non_unique / len(sumstats):.2f}%)')
sumstats = sumstats.filter(unique_mask)
# Print how many variants were retained
if verbose:
# noinspection PyUnboundLocalVariable
print(f'{len(sumstats):,} of {num_initial_variants:,} variants '
f'({100 * len(sumstats) / num_initial_variants:.2f}%) were '
f'retained')
# Sort
if BP is not None:
sumstats = sumstats.pipe(sort_sumstats)
# Save, if munged_sumstats_file is not None; block-gzip if
# munged_sumstats_file ends in .gz
if munged_sumstats_file is not None:
sumstats \
.with_columns(pl.selectors.float()
.map_elements('{:.12g}'.format,
return_dtype=pl.String)) \
.write_csv(munged_sumstats_file.removesuffix('.gz'),
separator='\t')
if munged_sumstats_file.endswith('.gz'):
run(f'bgzip -f {munged_sumstats_file.removesuffix(".gz")}')
# Return the sumstats, even if saving
return sumstats
def munge_regenie_sumstats(raw_sumstats_files, munged_sumstats_file, *,
quantitative=False, dbSNP=None):
"""
A wrapper for munge_sumstats() when all sumstats are in regenie format.
Args:
raw_sumstats_files: a sumstats file (or list thereof) in regenie format
munged_sumstats_file: a file path where the munged sumstats file will
be output
quantitative: are the sumstats files for quantitative traits?
dbSNP: a DataFrame returned by load_dbSNP(); must match the genome
build of raw_sumstats_files' BP column! If not specified, do not
infer rs numbers.
"""
munge_sumstats(raw_sumstats_files, munged_sumstats_file,
CHROM=pl.col.CHROM.cast(pl.String).replace({'23': 'X'}),
BP='GENPOS', SNP='ID', REF='ALLELE0', ALT='ALLELE1',
AAF='A1FREQ',
AAF_CASES=None if quantitative else 'A1FREQ_CASES',
AAF_CONTROLS=None if quantitative else 'A1FREQ_CONTROLS',
INFO='INFO', N='N',
N_CASES=None if quantitative else 'N_CASES',
N_CONTROLS=None if quantitative else 'N_CONTROLS',
BETA='BETA', SE='SE', P=10 ** -pl.col.LOG10P,
quantitative=quantitative,
separator=' ', dbSNP=dbSNP)
def get_minimal_representations_numpy(ref_vector, alt_vector, bp_vector):
"""
Converts variants - represented as 1D NumPy arrays ref_vector, alt_vector,
and optionally bp_vector - to their minimal representations.
Modifies bp_vector, ref_vector and alt_vector in-place.
cureffi.org/2014/04/24/converting-genetic-variants-to-their-minimal-
representation explains what a minimal representation is.
github.com/ericminikel/minimal_representation/blob/master/
minimal_representation.py contains the code this function is based on.
This operation only affects indels: by definition, single-nucleotide
variants are already in their minimal representation.
Args:
ref_vector: the reference alleles of the variants
alt_vector: the alternate alleles of the variants
bp_vector: the base pairs of the variants; may be None
"""
import numpy as np
assert isinstance(ref_vector, np.ndarray) and ref_vector.ndim == 1
assert isinstance(alt_vector, np.ndarray) and alt_vector.ndim == 1
assert len(alt_vector) == len(ref_vector)
if bp_vector is not None:
assert isinstance(bp_vector, np.ndarray) and bp_vector.ndim == 1
assert len(bp_vector) == len(ref_vector)
assert bp_vector.dtype in ('int32', int), bp_vector.dtype
int_type = 'int' if bp_vector.dtype == 'int32' else 'long'
# noinspection PyUnboundLocalVariable
cython_inline(rf'''
def get_minimal_representations_cython(str[:] ref_vector, str[:] alt_vector
{f', {int_type}[:] bp_vector' if bp_vector is not None else ''}):
cdef long i, min_len, ref_len, alt_len, ref_start, ref_end, \
alt_start, alt_end
{f'cdef {int_type} bp' if bp_vector is not None else ''}
cdef str ref, alt
for i in range(ref_vector.shape[0]):
ref = ref_vector[i]
alt = alt_vector[i]
ref_len = len(ref)
alt_len = len(alt)
min_len = ref_len if ref_len < alt_len else alt_len
if min_len == 1:
continue
{f'bp = bp_vector[i]' if bp_vector is not None else ''}
ref_start = 0
alt_start = 0
ref_end = ref_len - 1
alt_end = alt_len - 1
while ref[ref_end] == alt[alt_end]:
ref_end -= 1
alt_end -= 1
min_len -= 1
if min_len == 1: break
else:
while ref[ref_start] == alt[alt_start]:
ref_start += 1
alt_start += 1
{'bp += 1' if bp_vector is not None else ''}
min_len -= 1
if min_len == 1: break
ref_vector[i] = ref[ref_start:ref_end + 1]
alt_vector[i] = alt[alt_start:alt_end + 1]
{'bp_vector[i] = bp' if bp_vector is not None else ''}
''')['get_minimal_representations_cython'](
**{'ref_vector': ref_vector, 'alt_vector': alt_vector} |
(({'bp_vector': bp_vector}) if bp_vector is not None else {}))
def get_minimal_representations(df, *, ref_col='REF', alt_col='ALT',
bp_col='BP'):
"""
A wrapper for get_minimal_representations_numpy() for polars DataFrames.
Args:
df: a polars DataFrame
ref_col: the name of the reference allele column in df
alt_col: the name of the alternate allele column in df
bp_col: the name of the base-pair column in df; may be None
Returns:
A same-sized df with each variant converted to its minimal
representation.
"""
import pyarrow as pa
if not isinstance(df, pl.DataFrame):
raise ValueError('df must be a polars DataFrame!')
if df.is_empty():
raise ValueError(f'df is empty!')
if ref_col not in df:
raise ValueError(f'{ref_col!r} not in df; specify ref_col')
if alt_col not in df:
raise ValueError(f'{alt_col!r} not in df; specify alt_col')
if bp_col is not None and bp_col not in df:
raise ValueError(f'{bp_col!r} not in df; specify bp_col')
ref_vector = df[ref_col].to_numpy()
alt_vector = df[alt_col].to_numpy()
bp_vector = df[bp_col].to_numpy(writable=True) \
if bp_col is not None else None
get_minimal_representations_numpy(ref_vector, alt_vector, bp_vector)
df = df.with_columns(pl.from_arrow(pa.array(
ref_vector, type=pa.large_utf8())).alias(ref_col),
pl.from_arrow(pa.array(
alt_vector, type=pa.large_utf8())).alias(alt_col),
**{bp_col: pl.from_numpy(bp_vector)[:, 0]}
if bp_col is not None else {})
return df
def get_minimal_representations_awk(include_bp=True):
"""
An awk implementation of get_minimal_representations(). Assumes the
variant's reference allele, alternate allele and base-pair position are
stored in the variables ref, alt and bp.
Args:
include_bp: if False, do not include the correction for base-pair
position (include_bp=False is useful when there's no
base-pair column)
Returns:
A code string to be integrated into a larger awk command.
"""
import re
return re.sub(r'\s+', ' ', f'''
if (length(ref) > 1 || length(alt) > 1) {{
while (length(alt) > 1 && length(ref) > 1 &&
substr(alt, length(alt)) == substr(ref, length(ref))) {{
alt = substr(alt, 1, length(alt) - 1);
ref = substr(ref, 1, length(ref) - 1);
}}
while (length(alt) > 1 && length(ref) > 1 &&
substr(alt, 1, 1) == substr(ref, 1, 1)) {{
alt = substr(alt, 2);
ref = substr(ref, 2);
{"bp++;" if include_bp else ""}
}}
}}
''').strip().rstrip()
def load_dbSNP(genome_build, *,
dbSNP_dir=f'{get_base_data_directory()}/dbSNP'):
"""
Load all autosomal + chrX/Y/M variants in dbSNP, caching intermediate and
final results in the cache directory dbSNP_dir. Multiallelic variants are
split, then variants are converted to their minimal representations.
Variants that map to multiple genomic locations
(ncbi.nlm.nih.gov/snp/docs/rs_multi_mapping) are removed.
Must be run on a node with a large amount of memory! 120 GiB should be
sufficient. (You will need more - closer to 188 GiB, the size of a
Niagara compute node - if the dbSNP cache has not been generated. Also,
you will need to first run it on a login node, to download the files, then
when it runs out of memory, switch to a compute node. It is hacky.)
Args:
genome_build: the genome build (hg38 or hg19)
dbSNP_dir: the cache directory to store intermediate and final results.
Because generating the cache can take a long time, you will
probably want to leave this argument at its default value.
Returns:
A polars DataFrame with columns CHROM, BP, REF, ALT and RSID.
"""
check_valid_genome_build(genome_build)
dbSNP_file = os.path.join(dbSNP_dir, f'{genome_build}.tsv')
read_csv_kwargs = dict(separator='\t', schema_overrides={
'CHROM': pl.Categorical, 'BP': pl.Int32,
'REF': pl.Categorical('lexical'), 'ALT': pl.Categorical('lexical')})
if os.path.exists(dbSNP_file):
return pl.read_csv(dbSNP_file, **read_csv_kwargs)
else:
full_dbSNP_file = os.path.join(
dbSNP_dir, f'GCF_000001405.'
f'{40 if genome_build == "hg38" else 25}.gz')
if not os.path.exists(full_dbSNP_file):
raise_error_if_on_compute_node()
run(f'mkdir -p {dbSNP_dir} && '
f'wget https://ftp.ncbi.nih.gov/snp/latest_release/VCF/'
f'{os.path.basename(full_dbSNP_file)} -O {full_dbSNP_file}')
dbSNP_file_with_multimapping = os.path.join(
dbSNP_dir, f'{genome_build}_with_multimapping.tsv')
if not os.path.exists(dbSNP_file_with_multimapping):
dbSNP_chromosome_IDs = {
'NC_000001.11': 'chr1', 'NC_000002.12': 'chr2',
'NC_000003.12': 'chr3', 'NC_000004.12': 'chr4',
'NC_000005.10': 'chr5', 'NC_000006.12': 'chr6',
'NC_000007.14': 'chr7', 'NC_000008.11': 'chr8',
'NC_000009.12': 'chr9', 'NC_000010.11': 'chr10',
'NC_000011.10': 'chr11', 'NC_000012.12': 'chr12',
'NC_000013.11': 'chr13', 'NC_000014.9': 'chr14',
'NC_000015.10': 'chr15', 'NC_000016.10': 'chr16',
'NC_000017.11': 'chr17', 'NC_000018.10': 'chr18',
'NC_000019.10': 'chr19', 'NC_000020.11': 'chr20',
'NC_000021.9': 'chr21', 'NC_000022.11': 'chr22',
'NC_000023.11': 'chrX', 'NC_000024.10': 'chrY',
'NC_012920.1': 'chrM'
} if genome_build == 'hg38' else {
'NC_000001.10': 'chr1', 'NC_000002.11': 'chr2',
'NC_000003.11': 'chr3', 'NC_000004.11': 'chr4',
'NC_000005.9': 'chr5', 'NC_000006.11': 'chr6',
'NC_000007.13': 'chr7', 'NC_000008.10': 'chr8',
'NC_000009.11': 'chr9', 'NC_000010.10': 'chr10',
'NC_000011.9': 'chr11', 'NC_000012.11': 'chr12',
'NC_000013.10': 'chr13', 'NC_000014.8': 'chr14',
'NC_000015.9': 'chr15', 'NC_000016.9': 'chr16',
'NC_000017.10': 'chr17', 'NC_000018.9': 'chr18',
'NC_000019.9': 'chr19', 'NC_000020.10': 'chr20',
'NC_000021.8': 'chr21', 'NC_000022.10': 'chr22',
'NC_000023.10': 'chrX', 'NC_000024.9': 'chrY',
'NC_012920.1': 'chrM'}
dbSNP_awk_array = ''.join(
f'dbSNP_map["{ID}"] = "{chrom}"; '
for ID, chrom in dbSNP_chromosome_IDs.items())
# Memory optimization: seen is deleted after every chromosome
# (assumes chromosomes are contiguous in dbSNP, which they are)
run(f'awk \'BEGIN {{while ((getline line) && (line ~ /^##/)); '
f'print "CHROM", "BP", "REF", "ALT", "RSID"; {dbSNP_awk_array}'
f'}} $1 in dbSNP_map {{split($5, alts, ","); chrom = '
f'dbSNP_map[$1]; if (chrom != prev_chrom) delete seen; '
f'prev_chrom = chrom; for (i in alts) {{ bp = $2; ref = $4; '
f'alt = alts[i]; {get_minimal_representations_awk()}; '
f'if (!seen[bp":"ref":"alt":"$3]++) print chrom, bp, ref, '
f'alt, $3}}}}\' OFS="\t" <(zcat {full_dbSNP_file}) > '
f'{dbSNP_file_with_multimapping}')
# Remove multi-mapping variants; this is equivalent (aside from
# sorting) to dbSNP = dbSNP.unique(['RSID', 'REF', 'ALT'], keep='none')
dbSNP = pl.read_csv(dbSNP_file_with_multimapping, **read_csv_kwargs)
dbSNP = dbSNP \
.cast({'CHROM': pl.Enum([f'chr{i}' for i in range(1, 23)] +
['chrX', 'chrY', 'chrM'])}) \
.sort('RSID', 'REF', 'ALT') \
.with_columns(rle_id=pl.struct('RSID', 'REF', 'ALT').rle_id()) \
.filter(pl.col.rle_id.ne(pl.col.rle_id.shift(fill_value=-1)) &
pl.col.rle_id.ne(pl.col.rle_id.shift(-1, fill_value=-1))) \
.drop('rle_id') \
.sort('CHROM', 'BP', 'REF', 'ALT', 'RSID')
dbSNP.write_csv(dbSNP_file, separator='\t')
return dbSNP
def check_valid_dbSNP(dbSNP):
"""
Checks if a dbSNP instance is valid.
Args:
dbSNP: a dbSNP instance
"""
if not isinstance(dbSNP, pl.DataFrame):
raise TypeError(f'dbSNP must be a DataFrame returned by '
f'load_dbSNP(), but has type {type(dbSNP).__name__}')
if dbSNP.columns != ['CHROM', 'BP', 'REF', 'ALT', 'RSID']:
raise ValueError(f"dbSNP must be a DataFrame returned by load_dbSNP() "
f"with columns ['CHROM', 'BP', 'REF', 'ALT', 'RSID']")
if dbSNP.height < 1_000_000_000:
raise ValueError(f'dbSNP must be a DataFrame returned by '
f'load_dbSNP() and should have at least a billion '
f'rows, but yours has only {dbSNP.height} rows')
def get_rs_numbers(df, dbSNP, *, chrom_col='CHROM', bp_col='BP', ref_col='REF',
alt_col='ALT', rs_col='SNP', flip_col='FLIP',
fall_back_to_old_IDs=False, include_merged_IDs=False,
verbose=True):
"""
Given a DataFrame of variants with chrom_col, bp_col, ref_col, and alt_col
columns, adds a column rs_col to the DataFrame with the rs numbers (joined
with commas, in the rare case a variant has multiple). Also adds a column
flip_col, saying which variants needed to have their ref and alt alleles
flipped to match dbSNP; flipping is only attempted for single-nucleotide
variants.
If rs_col is already a column of df, variants present in dbSNP will have
their IDs in rs_col overwritten with their rs numbers. rs numbers for
variants not present in dbSNP (including multi-mapping variants; see
load_dbSNP()) will be set to null, unless fall_back_to_old_IDs=True, in
which case their original IDs will be retained.
If rs_col is not already a column of df, variants missing from dbSNP will
always have their rs numbers set to null, as there are no old IDs to fall
back to.
Args:
df: a DataFrame with chrom_col, bp_col, ref_col, and alt_col columns
dbSNP: a DataFrame returned by load_dbSNP(); must match the genome
build of df's bp_col!
chrom_col: the name of the chromosome column in df
bp_col: the name of the base-pair column in df
ref_col: the name of the reference allele column in df
alt_col: the name of the alternate allele column in df
rs_col: the name of the rs number column to be added to df
flip_col: the name of the flip column to be added to df. True where
alleles had to be flipped to match dbSNP, False where they
matched without flipping, null if the variant didn't match
dbSNP either with or without flipping
fall_back_to_old_IDs: if True, variants missing from dbSNP will retain
their original variant IDs, instead of having
them set to null. Requires chrom_col and bp_col
to already be present in df.
include_merged_IDs: determines how to handle the case where multiple rs
numbers have been merged into a single rs number.
If True, include all of them as a List[String]
column; if False, only include the merged ID (which
we assume is the lowest-numbered one).
verbose: whether to print what's happening at each step
Returns: df with two additional columns: rs_col, containing the rs numbers,
and flip_col, containing which variants were flipped.
"""
if df.is_empty():
raise ValueError(f'df is empty!')
if chrom_col not in df:
raise ValueError(f'"{chrom_col}" not in df; specify chrom_col')
if bp_col not in df:
raise ValueError(f'"{bp_col}" not in df; specify bp_col')
if ref_col not in df:
raise ValueError(f'"{ref_col}" not in df; specify ref_col')
if alt_col not in df:
raise ValueError(f'"{alt_col}" not in df; specify alt_col')
if flip_col in df:
raise ValueError(f'"{flip_col}" already in df; rename it or specify '
f'a different column name for flip_col')
if fall_back_to_old_IDs and rs_col not in df:
raise ValueError(f'You specified fall_back_to_old_IDs=True, but '
f'rs_col "{rs_col}" is not in df; specify it')
check_valid_dbSNP(dbSNP)
# Construct the rs number column piecewise by chromosome for efficiency
df = df.with_row_index()
rs_numbers = None
for df_chrom_ID in df[chrom_col].unique(maintain_order=True):
try:
chrom = standardize_chromosomes(df_chrom_ID)
except ValueError:
raise ValueError(f'df contains non-standard chromosome '
f'"{df_chrom_ID}"!')
if verbose:
print(f'Getting rs numbers for {chrom}...')
# Subset to chromosome; convert df to minimal representations
dbSNP_chrom = dbSNP.filter(pl.col.CHROM == chrom) \
.drop('CHROM') \
.rename({'BP': bp_col, 'REF': ref_col, 'ALT': alt_col,
'RSID': rs_col})
df_chrom = df.filter(pl.col(chrom_col) == df_chrom_ID) \
.select('index', bp_col, ref_col, alt_col) \
.pipe(get_minimal_representations, bp_col=bp_col, ref_col=ref_col,
alt_col=alt_col) \
.with_columns(pl.col(bp_col).cast(pl.Int32),
pl.col(ref_col, alt_col).cast(pl.Categorical))
# Allow matches without ref/alt flips...
matches_without_flips = df_chrom \
.join(dbSNP_chrom, on=[bp_col, ref_col, alt_col], how='left',
coalesce=True) \
.drop_nulls(rs_col)
# ...or with ref/alt flips, but only for SNVs (since for indels, e.g.
# "21:15847757:A:AG" and "21:15847757:AG:A" are different variants: the
# first is an insertion and the second is a deletion)
# Remove .cast(pl.String) once polars allows Categorical
# .str.len_bytes(): github.com/pola-rs/polars/issues/9773
matches_with_flips = df_chrom \
.filter(pl.col(ref_col).cast(pl.String).str.len_bytes() == 1,
pl.col(alt_col).cast(pl.String).str.len_bytes() == 1) \
.join(dbSNP_chrom, left_on=[bp_col, ref_col, alt_col],
right_on=[bp_col, alt_col, ref_col], how='left') \
.drop_nulls(rs_col)
# Ensure no variants match both with and without flips (theoretically
# possible since dbSNP has lots of edge cases, but we don't support it)
assert not matches_with_flips['index'].is_in(
matches_without_flips['index']).any()
# Merge matches with and without flips; in the rare case that a variant
# has multiple rs numbers, report all of them as a list column (if
# include_merged_IDs=True), or take only the lowest-numbered one (
# if include_merged_IDs=False)
chrom_rs_numbers = pl.concat([
matches_without_flips.with_columns(pl.lit(False).alias(flip_col)),
matches_with_flips.with_columns(pl.lit(True).alias(flip_col))]) \
.group_by('index') \
.agg(pl.col(rs_col).sort_by(pl.col(rs_col).str.slice(2).cast(int))
if include_merged_IDs else
pl.col(rs_col).sort_by(pl.col(rs_col).str.slice(2).cast(int))
.first(),
pl.first(flip_col))
rs_numbers = chrom_rs_numbers if rs_numbers is None else \
rs_numbers.extend(chrom_rs_numbers)
del dbSNP_chrom, df_chrom, matches_without_flips, matches_with_flips, \
chrom_rs_numbers
if rs_col in df:
df = df \
.lazy() \
.drop(rs_col) \
.join(rs_numbers.lazy(), on='index', how='left') \
.drop('index') \
.collect()
else:
df = df \
.join(rs_numbers, on='index', how='left') \
.drop('index')
# Print how many variants had rs numbers in dbSNP
if verbose:
num_mapped = len(df) - df[rs_col].null_count()
print(f'{num_mapped:,} of {len(df):,} variants '
f'({100 * num_mapped / len(df):.2f}%) had rs numbers in dbSNP')
# # Fall back to old IDs, if specified
if rs_col in df and fall_back_to_old_IDs:
df = df.with_columns(pl.col(rs_col).fill_null(df[rs_col]))
return df
def get_positions(df, dbSNP, *, rs_col='SNP', ref_col='REF', alt_col='ALT',
chrom_col='CHROM', bp_col='BP', flip_col='FLIP',
fall_back_to_old_positions=False):
"""
The reverse of get_rs_numbers(): given a DataFrame of variants with rs_col,
ref_col, and alt_col columns, adds columns chrom_col and bp_col to the
DataFrame giving the chromosome and base-pair positions of each variant.
Also adds a column flip_col, saying which variants needed to have their ref
and alt alleles flipped to match dbSNP; flipping is only attempted for
single-nucleotide variants.
Use this function to remap sumstats to a different genome build!
If chrom_col and bp_col are already columns of df, variants present in
dbSNP will have their chromosomes and base-pairs in chrom_col and bp_col
overwritten with those from dbSNP. Variants not present in dbSNP (including
multi-mapping variants; see load_dbSNP()) will have their chromosomes and
base-pair positions set to null, unless fall_back_to_old_positions=True, in
which case their original chromsomes/base-pair positions wil be retained.
If chrom_col and bp_col are not already columns of df, variants missing
from dbSNP will always have their chromosomes and base-pair positions set
to null, as there are no old chromosomes and base-pair positions to fall
back to.
Note: this function does not include defunct rs numbers that have been
merged into other rs numbers, since they don't appear in the dbSNP files
used in load_dbSNP(). Unlike for get_rs_numbers(), which also has this
behavior, here it is a limitation!
Args:
df: a DataFrame with rs_col, ref_col, and alt_col columns
dbSNP: a DataFrame returned by load_dbSNP() for the genome build you
want to get chrom/bp positions for
rs_col: the name of the rs number column in df
ref_col: the name of the reference allele column in df
alt_col: the name of the alternate allele column in df
chrom_col: the name of the chromosome column to be added to df
bp_col: the name of the base-pair column to be added to df
flip_col: the name of the flip column to be added to df. True where
alleles had to be flipped to match dbSNP, False where they
matched without flipping, null if the variant didn't match
dbSNP either with or without flipping
fall_back_to_old_positions: if True, variants missing from dbSNP will
retain their original chromosomes and
base-pair positions, instead of having them
set to null. Requires chrom_col and bp_col
to already be present in df.
Returns: df with three additional columns: chrom_col and bp_col, containing
the chromosomes and base-pair positions, and flip_col, containing
which variants were flipped.
"""
if df.is_empty():
raise ValueError(f'df is empty!')
if rs_col not in df:
raise ValueError(f'"{rs_col}" not in df; specify rs_col')
if ref_col not in df:
raise ValueError(f'"{ref_col}" not in df; specify ref_col')
if alt_col not in df:
raise ValueError(f'"{alt_col}" not in df; specify alt_col')
if chrom_col in df and bp_col not in df:
raise ValueError(f'chrom_col "{chrom_col}" is present in df but '
f'bp_col "{bp_col}" is not; either both must be '
f'present, or neither')
if chrom_col not in df and bp_col in df:
raise ValueError(f'bp_col "{bp_col}" is present in df but chrom_col '
f'"{chrom_col}" is not; either both must be present, '
f'or neither')
if flip_col in df:
raise ValueError(f'"{flip_col}" already in df; rename it or specify '
f'a different column name for flip_col')
check_valid_dbSNP(dbSNP)
# If fall_back_to_old_positions=True, save chromosomes and base-pair
# positions for later; otherwise, drop them so that they don't mess up the
# join with dbSNP
if fall_back_to_old_positions:
if chrom_col not in df:
raise ValueError(f'You specified fall_back_to_old_positions=True, '
f'but chrom_col "{chrom_col}" and bp_col '
f'"{bp_col}" are not in df; specify them')
df = df.rename({chrom_col: f'_GET_VARIANT_POSITIONS_{chrom_col}',
bp_col: f'_GET_VARIANT_POSITIONS_{bp_col}'})
else:
df = df.drop(chrom_col, bp_col)
# Convert df's ref and alt to their minimal representations
df = df \
.pipe(get_minimal_representations, ref_col=ref_col, alt_col=alt_col,
bp_col=None) \
.with_columns(pl.col(ref_col, alt_col).cast(pl.Categorical))
# Unlike for get_rs_numbers(), there isn't much of an efficiency gain in
# constructing the chromosome and base-pair columns piecewise by chromosome
# because the variant's chromosome isn't known a priori, so it needs to be
# matched against the entirety of dbSNP. However, as an optimization,
# subset dbSNP to just the rsIDs in df.
dbSNP = dbSNP \
.rename({'CHROM': chrom_col, 'BP': b

fibromyalgia_utils.py at commit de5c2c2, no license · at the source

Overview

Authors: Isabel Kerrebijn^1,2, Gyda Bjornsdottir³, Keon Arbabi^1,2,4, Lea Urpa^5,6,7, Hele Haapaniemi⁷, Gudmar Thorleifsson³, Lilja Stefansdottir³, Stephan Frangakis⁸, Jesse Valliere^5,7, Lovemore Kunorozva^6,7,9,10, Erik Abner¹¹, Caleb Ji^1,12, Markus Kangur¹, Bitten Aagaard¹³, Henning Bliddal¹⁴, Søren Brunak^15,16, Mie T Bruun¹⁷, Maria Didriksen¹⁸, Christian Erikstrup^19,20, Sarah Finer²¹

and 38 other authors

Arni J Geirsson²², Daniel F Gudbjartsson³, Thomas F Hansen^23,24,25, David van Heel^26,27, Ingileif Jonsdottir³, Stacey Knight²⁸, Kirk U Knowlton²⁸, Christina Mikkelsen¹⁸, Lincoln D Nadauld²⁹, Thorunn A Olafsdottir^3,30, Sisse R Ostrowski^18,31, Ole B V Pedersen^31,32, Saedis Saevarsdottir^3,30, Astros T Skuladottir^3,30, Erik Sørensen¹⁸, Hreinn Stefansson³, Patrick Sulem³, Olafur A Sveinsson^30,33, Gudny E Thorlacius³, Unnur Thorsteinsdottir³, Henrik Ullum³⁴, Arnor Vikingsson²², Thomas M Werge^31,35, Chronic Pain Genomics Consortium, FinnGen, DBDS Genomic Consortium, Estonian Biobank Research Team, Genes & Health Research Team, Richa Saxena^6,9,10,36, Kari Stefansson³, Chad M Brummett^8,37,38, Bente Glintborg^31,39, Daniel J Clauw⁴⁰, Thorgeir E Thorgeirsson³, Frances M K Williams⁴¹, Nasa Sinnott-Armstrong^{42,43,44,45,46}, Hanna M Ollila^6,7,9,36, Michael Wainberg^1,2,12,47,48

48 affiliations

Lunenfeld-Tanenbaum Research Institute, Mount Sinai Hospital, Toronto, Ontario Canada
Institute of Medical Science, University of Toronto, Toronto, Ontario Canada
Amgen deCODE Genetics, Reykjavik, Iceland
Krembil Centre for Neuroinformatics, Centre for Addiction and Mental Health, Toronto, Ontario Canada
Stanley Center for Psychiatric Research, Broad Institute, Cambridge, MA USA
Center for Genomic Medicine, Massachusetts General Hospital, Boston, MA USA
Institute for Molecular Medicine Finland (FIMM), HiLIFE, University of Helsinki, Helsinki, Finland
Department of Anesthesiology, University of Michigan Medical School, Ann Arbor, MI USA
Molecular and Population Genetics Program, Broad Institute, Cambridge, MA USA
Division of Sleep and Circadian Disorders, Brigham and Women’s Hospital, Boston, MA USA
Institute of Genomics, Estonian Genome Center, University of Tartu, Tartu, Estonia
Department of Computer Science, University of Toronto, Toronto, Ontario Canada
Department of Clinical Immunology, Aalborg University Hospital, Aalborg, Denmark
The Parker Institute, Copenhagen University Hospital Bispebjerg Frederiksberg, Copenhagen, Denmark
Department of Public Health, University of Copenhagen, Copenhagen, Denmark
Novo Nordisk Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
Clinical Immunology Research Unit, Department of Clinical Immunology, Odense University Hospital, Odense, Denmark
Department of Clinical Immunology, Copenhagen University Hospital, Rigshospitalet, Copenhagen, Denmark
Department of Clinical Immunology, Aarhus University Hospital, Aarhus, Denmark
Department of Clinical Medicine, Aarhus University, Aarhus, Denmark
Wolfson Institute of Population Health, Queen Mary University of London, London, UK
Department of Rheumatology, Landspitali University Hospital, Reykjavik, Iceland
Neurogenomics, Translational Research Centre, Copenhagen University Hospital, Glostrup, Denmark
Danish Multiple Sclerosis Center, Copenhagen University Hospital, Glostrup, Denmark
Danish Headache Center, Copenhagen University Hospital, Glostrup, Denmark
Blizard Institute, Queen Mary University of London, London, UK
Precision Healthcare University Research Institute, Queen Mary University of London, London, UK
Intermountain Medical Center, Intermountain Heart Institute, Salt Lake City, UT USA
Intermountain Healthcare, Saint George, UT USA
School of Health Sciences, Faculty of Medicine, University of Iceland, Reykjavik, Iceland
Department of Clinical Medicine, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark
Department of Clinical Immunology, Zealand University Hospital, Koege, Denmark
Department of Neurology, Landspitali University Hospital, Reykjavik, Iceland
Statens Serum Institut, Copenhagen, Denmark
Institute of Biological Psychiatry, Mental Health Services, Copenhagen University Hospital, Copenhagen, Denmark
Department of Anesthesiology, Mass General Brigham, Harvard Medical School, Boston, MA USA
Opioid Research Institute, Office for the Vice President for Research, University of Michigan, Ann Arbor, MI USA
Overdose Prevention Engagement Network, Institute for Healthcare Policy and Innovation, Michigan Medicine, Ann Arbor, MI USA
Copenhagen Center for Arthritis Research (COPECARE) and DANBIO, Center for Rheumatology and Spine Diseases, Rigshospitalet, Glostrup, Denmark
Chronic Pain and Fatigue Research Center, Department of Anesthesiology, University of Michigan, Ann Arbor, MI USA
Department of Twin Research and Genetic Epidemiology, School of Life Course Sciences, King’s College London, London, UK
Herbold Computational Biology Program, Fred Hutchinson Cancer Center, Seattle, WA USA
Department of Genome Sciences, University of Washington, Seattle, WA USA
Brotman Baty Institute, University of Washington, Seattle, WA USA
Center for Synthetic Biology, University of Washington, Seattle, WA USA
Center for One Health Research, University of Washington, Seattle, WA USA
Department of Psychiatry, University of Toronto, Toronto, Ontario Canada
Division of Biostatistics, Dalla Lana School of Public Health, University of Toronto, Toronto, Ontario Canada

Institutions: Mount Sinai Hospital (Canada); University of Toronto (Canada); Lunenfeld-Tanenbaum Research Institute (Canada); deCODE Genetics (Iceland) (Iceland); Centre for Addiction and Mental Health (Canada); Broad Institute (United States); University of Helsinki (Finland); Massachusetts General Hospital (United States); Institute for Molecular Medicine Finland (Finland); HiLIFE – Elämäntieteiden Instituutti; Stanley Center for Psychiatric Research; University of Michigan (United States); Michigan Medicine (United States); Brigham and Women's Hospital (United States); University of Tartu (Estonia); Aalborg University Hospital (Denmark); Bispebjerg Hospital (Denmark); Copenhagen University Hospital (Denmark); University of Copenhagen (Denmark); Odense University Hospital (Denmark); Rigshospitalet (Denmark); Aarhus University Hospital (Denmark); Aarhus University (Denmark); National University Hospital of Iceland (Iceland); Queen Mary University of London (United Kingdom); Intermountain Medical Center (United States); Danish Multiple Sclerosis Center (Denmark); Blizard Institute (United Kingdom); Intermountain Healthcare (United States); Zealand University Hospital (Denmark); Zealand University Hospital Køge (Denmark); University of Iceland (Iceland); Statens Serum Institut (Denmark); Mental Health Services (Denmark); Harvard University (United States); Mass General Brigham (United States); Glostrup Hospital (Denmark); King's College London (United Kingdom); University of Washington (United States); Fred Hutch Cancer Center (United States); Brotman Baty Institute (United States)

Journal: Nature medicine, volume 32, issue 8, pages 3060-3070

Dates: received 18 September 2025; accepted 28 May 2026; published online 28 July 2026; in print 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41591-026-04492-6 · PMID 42521817 · PMCID PMC13472937 · OpenAlex W4414345372

Open access: hybrid, a free copy (OpenAlex)

Preprint: doi:10.1101/2025.09.18.25335914

Status: code verified

Categories: genetics / omics (modality), human (organism), pain (population)

Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning

Keywords: Genome-wide association studies, Fibromyalgia

MeSH: Fibromyalgia*, Genetic Predisposition to Disease*, Female, Genome-Wide Association Study, Humans, Male, Polymorphism, Single Nucleotide, Stress Disorders, Post-Traumatic (* major topic)

Topic: Fibromyalgia and Chronic Fatigue Syndrome Research (Psychiatry and Mental health, Medicine), according to OpenAlex

Funding: Gouvernement du Canada | Instituts de Recherche en Santé du Canada | CIHR Skin Research Training Centre (Skin Research Training Centre) (FBD-199459, MHP-192163); European Commission (EC) (H2020-2020-848099); Novo Nordisk Fonden (Novo Nordisk Foundation) (NNF17OC0027864, NNF23OC0082015, NNF14CC0001, NNF17OC0027594); Oak Foundation (OFIL-24-074); Ministry of Education and Research | Estonian Research Competency Council (Research Competency Council) (PRG1291); NHGRI NIH HHS (RM1 HG010461); EC | Horizon 2020 Framework Programme (EU Framework Programme for Research and Innovation H2020) (894987, 101137201, 101137154); NIAMS NIH HHS (K08 AR082454); Natur og Univers, Det Frie Forskningsråd (Natural Sciences, Danish Council for Independent Research) (09-069412); U.S. Department of Health & Human Services | National Institutes of Health (NIH) (K08AR082454, RM1HG010461)

Citations: not cited yet (Europe PMC); 98 references in the paper

Abstract

Fibromyalgia is a common and debilitating chronic pain syndrome of poorly understood etiology. Here we conduct a multi-ancestry genome-wide association study meta-analysis across 2,563,755 individuals (54,629 cases and 2,509,126 controls) from 11 cohorts, identifying 26 risk loci for fibromyalgia. The strongest association was with a coding variant in HTT, the causal gene for Huntington’s disease. Gene prioritization implicated the HTT regulator GPR52, as well as diverse genes with neural roles, including DCC, DRD2/NCAM1, MDGA2 and CELF4. Fibromyalgia heritability was exclusively enriched within brain tissues and neural cell types. Fibromyalgia showed strong, positive genetic correlation with a wide range of chronic pain, psychiatric and somatic disorders, including genetic correlations above 0.7 with low back pain, post-traumatic stress disorder and irritable bowel syndrome. Despite large sex differences in fibromyalgia prevalence, the genetic architecture of fibromyalgia was nearly identical between males and females. This study provides robust genetic evidence defining fibromyalgia as a central nervous system disorder, thereby establishing a biological framework for its complex pathophysiology and extensive clinical comorbidities.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

FINNGEN/regenie-pipelines

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 320fd10df35be37d4b62977a03fbc86379cfa93c, 23 September 2026

Languages: Shell (5), R (2)

Size: 101 files, 7 scripts

Software Heritage: not archived

Found in: the text, “FinnGen”

Holds: README, license file, environment (requirements.txt, docker/Dockerfile)

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: data.table (2 files)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

9 files

scripts/build_docker.sh, Shell, 100 lines
scripts/conditional_release.sh, Shell, 104 lines
scripts/filter_hits_regions.sh, Shell, 92 lines
scripts/install.packages.R, R, 10 lines
scripts/qqplot.R, R, 134 lines
scripts/regenie_conditional.sh, Shell, 282 lines
scripts/return_bgen_chunks_limits.sh, Shell, 49 lines
LICENSE, License, 21 lines
README.md, Text, 232 lines

At the source: github.com/FINNGEN/regenie-pipelines

i-kerrebijn/fibromyalgia_GWAS_meta-analysis

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: de5c2c227a1a2176ca608105b25f1f27a5fdc9b7, 20 January 2026

Languages: Python (4)

Size: 5 files, 4 scripts

Software Heritage: not archived

Found in: “Code availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (1 file), NumPy (1 file), SciPy (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

4 files

LDSC_intercept.py, Python, 152 lines
fibromyalgia_utils.py, Python, 3,893 lines, 4 matches
meta_analysis.py, Python, 102 lines
munge_sumstats.py, Python, 395 lines

At the source: github.com/i-kerrebijn/fibromyalgia_GWAS_meta-analysis

Code availability

The code is available via GitHub at https://github.com/i-kerrebijn/fibromyalgia_GWAS_meta-analysis.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
11 scripts, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability Statement

GWAS summary statistics are available in the GWAS Catalog (accession codes GCST90838603 (https://www.ebi.ac.uk/gwas/search?query=GCST90838603) to GCST90838628 (https://www.ebi.ac.uk/gwas/search?query=GCST90838628)) and PRSs are available in the PGS catalog (score IDs PGS012560 and PGS012561). All summary statistics and PRSs are also available at https://paingenomics.org.

This study used data from the All of Us Research Program’s Controlled Tier Dataset Curated Data Repository version 8, available to authorized users on the Researcher Workbench.

The code is available via GitHub at https://github.com/i-kerrebijn/fibromyalgia_GWAS_meta-analysis.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

Publisher: n/a → Nature Portfolio

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 58 authors, 2 keywords, 8 MeSH terms, 10 funders, 97 references.

Cite

This paper

Kerrebijn, I., Bjornsdottir, G., Arbabi, K., Urpa, L., Haapaniemi, H., Thorleifsson, G., Stefansdottir, L., Frangakis, S., Valliere, J., Kunorozva, L., Abner, E., Ji, C., Kangur, M., Aagaard, B., Bliddal, H., Brunak, S., Bruun, M. T., Didriksen, M., Erikstrup, C., . . . Wainberg, M. (2026). The genetic architecture of fibromyalgia across 2.5 million individuals. Nature medicine, 32(8), 3060-3070. https://doi.org/10.1038/s41591-026-04492-6

BibTeX

@article{kerrebijn2026genetic, author = {Kerrebijn, Isabel and Bjornsdottir, Gyda and Arbabi, Keon and Urpa, Lea and Haapaniemi, Hele and Thorleifsson, Gudmar and Stefansdottir, Lilja and Frangakis, Stephan and Valliere, Jesse and Kunorozva, Lovemore and Abner, Erik and Ji, Caleb and Kangur, Markus and Aagaard, Bitten and Bliddal, Henning and Brunak, Søren and Bruun, Mie T and Didriksen, Maria and Erikstrup, Christian and Finer, Sarah and Geirsson, Arni J and Gudbjartsson, Daniel F and Hansen, Thomas F and van Heel, David and Jonsdottir, Ingileif and Knight, Stacey and Knowlton, Kirk U and Mikkelsen, Christina and Nadauld, Lincoln D and Olafsdottir, Thorunn A and Ostrowski, Sisse R and Pedersen, Ole B V and Saevarsdottir, Saedis and Skuladottir, Astros T and Sørensen, Erik and Stefansson, Hreinn and Sulem, Patrick and Sveinsson, Olafur A and Thorlacius, Gudny E and Thorsteinsdottir, Unnur and Ullum, Henrik and Vikingsson, Arnor and Werge, Thomas M and {Chronic Pain Genomics Consortium} and {FinnGen} and {DBDS Genomic Consortium} and {Estonian Biobank Research Team} and {Genes \& Health Research Team} and Saxena, Richa and Stefansson, Kari and Brummett, Chad M and Glintborg, Bente and Clauw, Daniel J and Thorgeirsson, Thorgeir E and Williams, Frances M K and Sinnott-Armstrong, Nasa and Ollila, Hanna M and Wainberg, Michael}, title = {{The genetic architecture of fibromyalgia across 2.5 million individuals}}, journal = {Nature medicine}, year = {2026}, month = jul, volume = {32}, number = {8}, pages = {3060--3070}, publisher = {Nature Portfolio}, issn = {1078-8956}, doi = {10.1038/s41591-026-04492-6}, url = {https://doi.org/10.1038/s41591-026-04492-6}, pmid = {42521817}, pmcid = {PMC13472937} }

RIS

TY - JOUR AU - Kerrebijn, Isabel AU - Bjornsdottir, Gyda AU - Arbabi, Keon AU - Urpa, Lea AU - Haapaniemi, Hele AU - Thorleifsson, Gudmar AU - Stefansdottir, Lilja AU - Frangakis, Stephan AU - Valliere, Jesse AU - Kunorozva, Lovemore AU - Abner, Erik AU - Ji, Caleb AU - Kangur, Markus AU - Aagaard, Bitten AU - Bliddal, Henning AU - Brunak, Søren AU - Bruun, Mie T AU - Didriksen, Maria AU - Erikstrup, Christian AU - Finer, Sarah AU - Geirsson, Arni J AU - Gudbjartsson, Daniel F AU - Hansen, Thomas F AU - van Heel, David AU - Jonsdottir, Ingileif AU - Knight, Stacey AU - Knowlton, Kirk U AU - Mikkelsen, Christina AU - Nadauld, Lincoln D AU - Olafsdottir, Thorunn A AU - Ostrowski, Sisse R AU - Pedersen, Ole B V AU - Saevarsdottir, Saedis AU - Skuladottir, Astros T AU - Sørensen, Erik AU - Stefansson, Hreinn AU - Sulem, Patrick AU - Sveinsson, Olafur A AU - Thorlacius, Gudny E AU - Thorsteinsdottir, Unnur AU - Ullum, Henrik AU - Vikingsson, Arnor AU - Werge, Thomas M AU - Chronic Pain Genomics Consortium AU - FinnGen AU - DBDS Genomic Consortium AU - Estonian Biobank Research Team AU - Genes & Health Research Team AU - Saxena, Richa AU - Stefansson, Kari AU - Brummett, Chad M AU - Glintborg, Bente AU - Clauw, Daniel J AU - Thorgeirsson, Thorgeir E AU - Williams, Frances M K AU - Sinnott-Armstrong, Nasa AU - Ollila, Hanna M AU - Wainberg, Michael TI - The genetic architecture of fibromyalgia across 2.5 million individuals T2 - Nature medicine J2 - Nat Med PY - 2026 DA - 2026/07/28 VL - 32 IS - 8 SP - 3060 EP - 3070 SN - 1078-8956 PB - Nature Portfolio DO - 10.1038/s41591-026-04492-6 UR - https://doi.org/10.1038/s41591-026-04492-6 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41591-026-04492-6", "type": "article-journal", "title": "The genetic architecture of fibromyalgia across 2.5 million individuals", "container-title": "Nature medicine", "author": [ { "family": "Kerrebijn", "given": "Isabel" }, { "family": "Bjornsdottir", "given": "Gyda" }, { "family": "Arbabi", "given": "Keon" }, { "family": "Urpa", "given": "Lea" }, { "family": "Haapaniemi", "given": "Hele" }, { "family": "Thorleifsson", "given": "Gudmar" }, { "family": "Stefansdottir", "given": "Lilja" }, { "family": "Frangakis", "given": "Stephan" }, { "family": "Valliere", "given": "Jesse" }, { "family": "Kunorozva", "given": "Lovemore" }, { "family": "Abner", "given": "Erik" }, { "family": "Ji", "given": "Caleb" }, { "family": "Kangur", "given": "Markus" }, { "family": "Aagaard", "given": "Bitten" }, { "family": "Bliddal", "given": "Henning" }, { "family": "Brunak", "given": "Søren" }, { "family": "Bruun", "given": "Mie T" }, { "family": "Didriksen", "given": "Maria" }, { "family": "Erikstrup", "given": "Christian" }, { "family": "Finer", "given": "Sarah" }, { "family": "Geirsson", "given": "Arni J" }, { "family": "Gudbjartsson", "given": "Daniel F" }, { "family": "Hansen", "given": "Thomas F" }, { "family": "van Heel", "given": "David" }, { "family": "Jonsdottir", "given": "Ingileif" }, { "family": "Knight", "given": "Stacey" }, { "family": "Knowlton", "given": "Kirk U" }, { "family": "Mikkelsen", "given": "Christina" }, { "family": "Nadauld", "given": "Lincoln D" }, { "family": "Olafsdottir", "given": "Thorunn A" }, { "family": "Ostrowski", "given": "Sisse R" }, { "family": "Pedersen", "given": "Ole B V" }, { "family": "Saevarsdottir", "given": "Saedis" }, { "family": "Skuladottir", "given": "Astros T" }, { "family": "Sørensen", "given": "Erik" }, { "family": "Stefansson", "given": "Hreinn" }, { "family": "Sulem", "given": "Patrick" }, { "family": "Sveinsson", "given": "Olafur A" }, { "family": "Thorlacius", "given": "Gudny E" }, { "family": "Thorsteinsdottir", "given": "Unnur" }, { "family": "Ullum", "given": "Henrik" }, { "family": "Vikingsson", "given": "Arnor" }, { "family": "Werge", "given": "Thomas M" }, { "literal": "Chronic Pain Genomics Consortium" }, { "literal": "FinnGen" }, { "literal": "DBDS Genomic Consortium" }, { "literal": "Estonian Biobank Research Team" }, { "literal": "Genes & Health Research Team" }, { "family": "Saxena", "given": "Richa" }, { "family": "Stefansson", "given": "Kari" }, { "family": "Brummett", "given": "Chad M" }, { "family": "Glintborg", "given": "Bente" }, { "family": "Clauw", "given": "Daniel J" }, { "family": "Thorgeirsson", "given": "Thorgeir E" }, { "family": "Williams", "given": "Frances M K" }, { "family": "Sinnott-Armstrong", "given": "Nasa" }, { "family": "Ollila", "given": "Hanna M" }, { "family": "Wainberg", "given": "Michael" } ], "container-title-short": "Nat Med", "volume": "32", "issue": "8", "page": "3060-3070", "DOI": "10.1038/s41591-026-04492-6", "PMID": "42521817", "PMCID": "PMC13472937", "ISSN": "1078-8956", "publisher": "Nature Portfolio", "URL": "https://doi.org/10.1038/s41591-026-04492-6", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 28 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1073/pnas.2609814123 [code]: Genetic architectures of brain-related traits are shaped by strong selective constraints.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: data.table, SciPy, Matplotlib, 1 other tool, genetics / omics, 7 references
[2] doi:10.1038/s41467-026-73428-y [code]: Regional heterogeneity in phenotypic and genetic associations between bone and brain in humans.
Journal: Nature communications
In common: data.table, SciPy, NumPy, genetics / omics, 7 references
[3] doi:10.2147/jpr.s619982 [code]: Spatially Contextualized Integrative Genomics Highlights Neuronal and Glial Regulatory Programs in Low Back Pain.
Journal: Journal of pain research
In common: data.table, SciPy, Matplotlib, 1 other tool, pain, genetics / omics, 5 references
[4] doi:10.1038/s41467-026-72164-7 [code]: Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders.
Journal: Nature communications
In common: data.table, SciPy, NumPy, genetics / omics, 5 references
[5] doi:10.1038/s41562-026-02486-5 [code]: Genome-wide association studies of infant and toddler temperament in European and multi-ancestry populations.
Journal: Nature human behaviour
In common: data.table, SciPy, Matplotlib, 1 other tool, genetics / omics, 4 references
[6] doi:10.1038/s41398-026-04137-9 [code]: An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression.
Journal: Translational psychiatry
In common: data.table, SciPy, Matplotlib, 1 other tool, genetics / omics, 3 references
[7] doi:10.1038/s41467-026-76676-0 [code]: Determinants of functional burden pleiotropy and gene dosage responses across human traits.
Journal: Nature communications
In common: data.table, SciPy, Matplotlib, 1 other tool, genetics / omics, 3 references
[8] doi:10.1038/s41467-026-73996-z [code]: Genetic architecture of white matter microstructure captured by unsupervised deep representation learning of fractional anisotropy maps.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, genetics / omics, 4 references
[9] doi:10.1073/pnas.2516601123 [code]: Unveiling the glymphatic system's role in brain aging: A comprehensive biomarker and modifiable intervention target.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: data.table, SciPy, Matplotlib, 1 other tool, 3 references
[10] doi:10.1186/s12967-026-08266-z [code]: Single-cell multi-omic integration analysis prioritizes druggable genes and reveals cell-type-specific causal effects in glioblastomagenesis.
Journal: Journal of translational medicine
In common: data.table, genetics / omics, 4 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

The genetic architecture of fibromyalgia across 2.5 million individuals.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

FINNGEN/regenie-pipelines

i-kerrebijn/fibromyalgia_GWAS_meta-analysis

Code availability

Tracing map

Data

Data Availability Statement

Versions

Version 2, 28 September 2026

Version 1, 27 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity