获取数量性状位点或药靶数据

Xing Abao Lv3

安装程序

To run mtag, you will need to have Python 2.7 installed with the following packages:

  • numpy (>=1.13.1)
  • scipy
  • pandas (>=0.18.1)
  • argparse
  • bitarray (for ldsc)
  • joblib

(Note: if you already have the Python 3 version of the Anaconda distribution installed, then you will need to create and activate a Python 2.7 environment to run mtag. See here for details.)

mtag may be downloaded by cloning this github repository:

1
2
git clone https://github.com/omeed-maghzian/mtag.git
cd mtag

该包可以直接从 GitHub 上安装,使用remotes包。

1
remotes::install_github("alexploner/cfdr.pleio")

cfdr.pleio实现的算法需要一个参考数据集,该数据集描述了已知人类基因组中相当大比例变异体的局部 LD 结构。有关如何从头构建此参考数据的描述,请参见其他文档。在本文章中,我们将简单下载一个预先计算的参考数据集,下载地址:https://zenodo.org/record/5750318/files/genref.zip

请注意,该通用参考数据集对于cfdr.pleio的工作是必需的。由于其占用约 3 GB 的空间,通常最好将其存储在安全的地方,与使用cfdr.pleio的任何特定项目分开(因为可能不想重复下载)。这里,我们使用一个用户主目录下的建议目录名称,但您可以选择任何适合您的名称:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# https://zenodo.org/record/5750318/files/genref.zip
# 建议手动用浏览器下载,需要科学上网

# 解压缩文件并检查其内容:
REF_DIR = "E:/QTLMR/REF_DIR"
zip_file = paste0(REF_DIR, '/', 'genref.zip')
unzip(zip_file, exdir = REF_DIR)

dir(REF_DIR)
# [1] "all_chr_perVariant.rds" "chr01_LDpairs.rds" "chr02_LDpairs.rds" "chr03_LDpairs.rds" "chr04_LDpairs.rds"
# [6] "chr05_LDpairs.rds" "chr06_LDpairs.rds" "chr07_LDpairs.rds" "chr08_LDpairs.rds" "chr09_LDpairs.rds"
# [11] "chr10_LDpairs.rds" "chr11_LDpairs.rds" "chr12_LDpairs.rds" "chr13_LDpairs.rds" "chr14_LDpairs.rds"
# [16] "chr15_LDpairs.rds" "chr16_LDpairs.rds" "chr17_LDpairs.rds" "chr18_LDpairs.rds" "chr19_LDpairs.rds"
# [21] "chr20_LDpairs.rds" "chr21_LDpairs.rds" "chr22_LDpairs.rds" "genref.zip"

应该看到 23 个 .rda 文件,其中一个包含所有参考变异体的列表,其余 22 个描述非性染色体之间的 LD 结构。

截至目前已经安装了运行 condFdr 或 conjFdr 分析所需的一切

简单示例:BMI 和身高

使用 GIANT 发布的数据,我们希望利用我们的 cfdr.pleio 来回答两个问题:

  1. 在身高的条件下,有多少 SNP 与 BMI 相关,条件假阳性发现率condFdr < 0.001或更低,与常规全基因组显著性阈值p < 5E-8相比?
  2. 我们可以预期有多少 SNP 与身高和 BMI 都相关,联合假阳性发现率阈值为conjFdr < 0.01

加载 R 包

1
2
suppressWarnings(suppressMessages(library(data.table)))
suppressWarnings(suppressMessages(library(cfdr.pleio)))

数据下载

1
2
3
4
5
# BMI
https://portals.broadinstitute.org/collaboration/giant/images/1/15/SNP_gwas_mc_merge_nogc.tbl.uniq.gz

# Height
https://portals.broadinstitute.org/collaboration/giant/images/0/01/GIANT_HEIGHT_Wood_et_al_2014_publicrelease_HapMapCeuFreq.txt.gz

选择这两个数据是因为其大小适中,具有强烈的遗传信号,并且可以直接下载

读取数据

接下来,我们将数据读入 R,并做一些准备,以便将其输入cfdr.pleio,此时它对预期的列名仍然很严格。

这里使用data.table进行导入,主要是由于该包不仅具有快速从压缩文本文件读取的功能,cfdr.pleio也使用其数据格式进行内部存储,因此这相当高效。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
BMI = fread("E:/QTLMR/GWAS/SNP_gwas_mc_merge_nogc.tbl.uniq.gz")

BMI
# SNP A1 A2 Freq1.Hapmap b se p N
# <char> <char> <char> <num> <num> <num> <num> <num>
# 1: rs1000000 G A 0.6333 0.0001 0.0044 0.98190 231410.0
# 2: rs10000010 T C 0.5750 -0.0029 0.0030 0.33740 322079.0
# 3: rs10000012 G C 0.1917 -0.0095 0.0054 0.07853 233933.0
# 4: rs10000013 A C 0.8333 -0.0095 0.0044 0.03084 233886.0
# 5: rs10000017 C T 0.7667 -0.0034 0.0046 0.45980 233146.0
# ---
# 2554633: rs9999992 A G 0.0500 -0.0055 0.0124 0.65740 172167.0
# 2554634: rs9999993 T A 0.4583 -0.0063 0.0037 0.08862 234013.0
# 2554635: rs9999996 C A 0.1917 -0.0035 0.0053 0.50900 233967.0
# 2554636: rs9999997 A G 0.5500 -0.0055 0.0037 0.13720 233380.0
# 2554637: rs9999998 T C NA -0.0019 0.0085 0.82310 77362.5
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Height = fread("E:/QTLMR/GWAS/GIANT_HEIGHT_Wood_et_al_2014_publicrelease_HapMapCeuFreq.txt.gz")

Height
# MarkerName Allele1 Allele2 Freq.Allele1.HapMapCEU b SE p N
# <char> <char> <char> <num> <num> <num> <num> <int>
# 1: rs4747841 A G 0.551 -0.0011 0.0029 7.0e-01 253213
# 2: rs4749917 T C 0.436 0.0011 0.0029 7.0e-01 253213
# 3: rs737656 A G 0.367 -0.0062 0.0030 4.2e-02 253116
# 4: rs737657 A G 0.358 -0.0062 0.0030 4.1e-02 252156
# 5: rs7086391 T C 0.120 -0.0087 0.0038 2.4e-02 248425
# ---
# 2550854: rs4445756 T G NA -0.0070 0.0083 4.0e-01 67223
# 2550855: rs4299144 A G NA 0.0037 0.0058 5.3e-01 67210
# 2550856: rs2841648 A C NA 0.0140 0.0030 7.3e-06 249974
# 2550857: rs2468460 A G NA 0.0017 0.0065 7.9e-01 56490
# 2550858: rs2933064 C G NA -0.0073 0.0062 2.3e-01 55329

从上可以看到这两个数据集具有相同的基本格式,尽管列名略有不同,且 SNP 的数量也大致相同。

修正列名

有效的cfdr.pleio性状数据集需要有三列,具体名称以及内容如下:

  1. SNP,变异位点的 rs 标识符;
  2. BETA,基于底层回归模型的 SNP 的估计效应值;
  3. PVAL,效应估计对应的 P 值;
  4. 所有其他变量将被忽略。
1
2
3
4
5
6
7
8
9
10
11
colnames(BMI)[c(1, 5, 7)] = c("SNP", "BETA", "PVAL")

head(BMI)
# SNP A1 A2 Freq1.Hapmap BETA se PVAL N
# <char> <char> <char> <num> <num> <num> <num> <num>
# 1: rs1000000 G A 0.6333 0.0001 0.0044 0.98190 231410
# 2: rs10000010 T C 0.5750 -0.0029 0.0030 0.33740 322079
# 3: rs10000012 G C 0.1917 -0.0095 0.0054 0.07853 233933
# 4: rs10000013 A C 0.8333 -0.0095 0.0044 0.03084 233886
# 5: rs10000017 C T 0.7667 -0.0034 0.0046 0.45980 233146
# 6: rs10000023 G T 0.4083 0.0024 0.0038 0.52770 233860
1
2
3
4
5
6
7
8
9
10
11
colnames(Height)[c(1, 5, 7)] = c("SNP", "BETA", "PVAL")

head(Height)
# SNP Allele1 Allele2 Freq.Allele1.HapMapCEU BETA SE PVAL N
# <char> <char> <char> <num> <num> <num> <num> <int>
# 1: rs4747841 A G 0.551 -0.0011 0.0029 7.0e-01 253213
# 2: rs4749917 T C 0.436 0.0011 0.0029 7.0e-01 253213
# 3: rs737656 A G 0.367 -0.0062 0.0030 4.2e-02 253116
# 4: rs737657 A G 0.358 -0.0062 0.0030 4.1e-02 252156
# 5: rs7086391 T C 0.120 -0.0087 0.0038 2.4e-02 248425
# 6: rs878177 T C 0.300 0.0140 0.0031 8.2e-06 251271

此时,我们可以对 SNP 进行一些过滤,例如基于质量分数、主要等位基因频率或样本大小等等。不过这个示例数据作者在分析过程中已经经过了严格的过滤,未提提供太多的其他信息,因此这里跳过这一步。

是否需要进行过滤将根据项目具体情况进行判断。

查看数据

有多少 SNP 在 BMI 上达到了传统的全基因组显著性(与问题 1 相关)

1
2
3
4
5
6
7
8
9
table(BMI$PVAL < 5E-8)
#
# FALSE TRUE
# 2552777 1860

proportions(table(BMI$PVAL < 5E-8))
#
# FALSE TRUE
# 0.9992719122 0.0007280878

一共发现有 1860 个 SNP,未进行连锁不平衡分析,对应于 Locke 等人在文章中报告的 77 个良好分离的位点。

相比之下,发现与身高相关的 SNP 数量惊人地达到了 26593 个,这占所有 SNP 的 1% 左右。

1
2
3
4
5
6
7
8
9
table(Height$PVAL < 5E-8)
#
# FALSE TRUE
# 2524265 26593

proportions(table(Height$PVAL < 5E-8))
#
# FALSE TRUE
# 0.98957488 0.01042512

运行分析

现在我们已准备好开始实际分析,创建一个新的分析对象:

1
BMI_Height = cfdr_pleio$new()

该对象最初是空的,后续的分析将通过调用适当的方法进行,这与大多数其他 R 包略有不同,因为 cfdr.pleio 基于 R6 类系统。

初始化数据

第一步是用必要的信息初始化空对象以运行分析:

1. 保存两个感兴趣性状(本例中为 BMI 和身高)汇总统计信息的数据对象,

2. 下载的通用参考数据的位置(在本例中,这是上面存储为 REF_DIR 的目录),

3. 针对特定数据的参考数据的工作版本的存储位置。

后两者之间的区别很重要:我们之前下载的通用参考数据是 cfdr.pleio 设置的一部分,并且与应用 cfdr.pleio 的项目无关。特定参考数据是在分析的下一步生成的,本质上是与项目数据匹配的参考数据的简化版本;出于内存效率的考虑,此简化版本也保存到磁盘,但通常应该作为项目文件夹的一部分进行存储。对于我们的示例:

1
2
3
4
5
6
7
BMI_Height$init_data(
trait1 = BMI, trait2 = Height,
trait_names = c("BMI", "Height"),
refdat = refdata_location(REF_DIR),
local_refdat_path = "E:/tmp/cond_conjFDR/BMI_height",
verbose = TRUE
)

生成本地参考数据的过程相对耗时,请耐心等待。

初始化修剪索引

第二步是初始化随机修剪索引。这样做的动机在 Andreassen 等人的论文中提及,但基本上,条件 FDR 应从近似独立的变异体中估计。这是通过从两个性状共享的变异体集合中选择随机变异体,并系统性地丢弃与所选变异体在指定 LD 范围内的变异体来完成的;这个过程生成了一个近似独立的变异体子集。为了最小化随机选择的影响,使用不同的随机选择重复进行相当多的次数。对于我们的示例,我选择了n = 50次迭代:

1
BMI_Height$initialize_pruning_index(n_iter = 50, seed = 154226, verbose = TRUE)

计算 FDR

在第三步也是最后一步,我们将指定的数据和随机修剪索引结合在一起,计算条件和联合 FDR。在我们的示例中,我们特别希望知道 BMI 性状的条件 FDR,条件是身高性状。这可以通过指定哪个性状是 FDR 性状(即主要的关注性状),哪个是次要的条件性状来完成;因为 BMI 是第一个性状,我们这样做:

1
BMI_Height$calculate_cond_fdr(fdr_trait = 1, verbose = TRUE)

为了计算联合 FDR,我们需要将条件反转,重新进行计算,即也需要计算身高的条件 FDR,条件为 BMI:

1
BMI_Height$calculate_cond_fdr(fdr_trait = 2, verbose = TRUE)

此时,所有的条件 FDR 和联合 FDR 已被计算并存储在分析对象中,可以提取为 data.table

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
BMI_Height_res = BMI_Height$get_trait_results()

BMI_Height_res
# SNP BETA1 LOG10PVAL1 BETA2 LOG10PVAL2 cfdr12 cfdr21 conj_fdr
# <char> <num> <num> <num> <num> <num> <num> <num>
# 1: rs2153920 0.0135 1.96417017 -0.0081 1.26760624 0.2907734 0.5023832 0.5023832
# 2: rs6604648 -0.0010 0.09550084 0.0020 0.27572413 0.9944421 0.9808161 0.9944421
# 3: rs1277751 0.0012 0.08751224 0.0004 0.03621217 0.9971372 0.9981689 0.9981689
# 4: rs12670580 0.0073 0.49852993 -0.0041 0.31875876 0.9479804 0.9590807 0.9590807
# 5: rs12893676 -0.0004 0.02291688 0.0160 3.11350927 0.9899864 0.2864715 0.9899864
# ---
# 2458100: rs13259729 -0.0248 0.79263496 0.0099 0.32790214 0.8840471 0.9387551 0.9387551
# 2458101: rs13156977 -0.0078 0.85047299 -0.0052 0.67778071 0.8188316 0.8501499 0.8501499
# 2458102: rs13026991 -0.0031 0.38237070 -0.0010 0.13667714 0.9758350 0.9878107 0.9878107
# 2458103: rs1346263 0.0029 0.23336411 -0.0008 0.07058107 0.9912054 0.9950950 0.9950950
# 2458104: rs1323344 -0.0029 0.35124979 -0.0022 0.33724217 0.9714638 0.9658762 0.9714638

此对象返回了指定的两个性状的原始数据,以及在分析过程中计算的所有新 FDR 列(此处为三个,取其最大值)。

保存结果

此时,保存结果以便进一步处理是个好主意。主要结果是包含所有设置、种子等信息的完整分析对象,但也可以方便地只保存最终的 data.table 以便进一步处理:

1
2
3
DATA_DIR = "E:/tmp/cond_conjFDR/BMI_height"
saveRDS(BMI_Height, file = paste(DATA_DIR, "BMI_Height_AnalysisObject.rds", sep = "/"))
saveRDS(BMI_Height_res, file = paste(DATA_DIR, "BMI_Height_ResultsTable.rds", sep = "/"))

结果解读

问题 1

让我们看看在条件身高的情况下,BMI 的条件 FDR 小于 0.001 的变异体数量:

1
2
3
4
table(BMI_Height_res$cfdr12 < 0.001)
#
# FALSE TRUE
# 2455319 2785

与全基因组显著性相比,我们发现大约有 50% 的 SNP 在这个保守的 FDR 阈值下显著,即相比单纯用全基因组显著性发现的更多,提升了检测能力。这说明了,condFDR 可以利用”多效性”信息,发现更多可能与 BMI 相关的 SNP

我们还可以反转这个问题:在对 BMI 具有(无条件)全基因组显著性的变异中,与身高相关的最大条件 FDR 是多少?

1
2
3
summary(subset(BMI_Height_res, LOG10PVAL1 > -log10(5E-8))$cfdr12)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 5.444e-08 1.469e-07 2.583e-07 7.862e-06 4.263e-06 1.160e-04

在这组 SNP 中,BMI 的最大条件 FDR 约为 0.00012,这非常保守;我们可以放宽至 0.001 的阈值,增加约 50% 的功效,同时仍能对识别的 SNP 保持非常保守的 FDR 控制。

问题 2

我们发现有n = 892个 SNP 似乎在联合 FDR 水平为 0.01 或更低的情况下与 BMI 和身高相关:

1
2
3
4
table(BMI_Height_res$conj_fdr < 0.01)
#
# FALSE TRUE
# 2457212 892

让我们查看这些重叠的 SNP:我们为这一子集的 SNP 绘制两个性状的原始 p 值的散点图:

1
2
3
4
5
6
7
8
9
10
11
plot(
LOG10PVAL2 ~ LOG10PVAL1,
data = BMI_Height_res[conj_fdr < 0.01],
xlab = "-log10(p) for BMI",
ylab = "-log10(p) for height",
xlim = c(0,50),
ylim = c(0, 50),
pch = 19,
col = gray(0.4, 0.25)
)
abline(h = -log10(5E-8), v = -log10(5E-8), lty = 3)

每个点代表 892 个多效性 SNP,而虚线表示各性状原始 p 值的全基因组截断线。我们发现联合 FDR 找到的重叠远大于各性状全基因组显著变异体的简单交集(如右上象限所示);此外,SNP 的增益并不是对称的:对 BMI 的条件关联证据较弱的 SNP 更多地被纳入(左上象限),而身高则较少(右下象限),这并不令人惊讶,因为身高的关联相对广泛。

conjFDR 发现的变异体要多于两个性状各自全基因组显著 SNP 的交集,这说明联合分析更敏感,通过散点图,可以发现很多 SNP 在单个性状下可能不够显著,但是联合分析下被纳入了。

代码简洁版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
REF_DIR = "E:/QTLMR/REF_DIR"
zip_file = paste0(REF_DIR, '/', 'genref.zip')
unzip(zip_file, exdir = REF_DIR)
dir(REF_DIR)

suppressWarnings(suppressMessages(library(data.table)))
suppressWarnings(suppressMessages(library(cfdr.pleio)))

BMI = fread("E:/QTLMR/GWAS/SNP_gwas_mc_merge_nogc.tbl.uniq.gz")
BMI

Height = fread("E:/QTLMR/GWAS/GIANT_HEIGHT_Wood_et_al_2014_publicrelease_HapMapCeuFreq.txt.gz")
Height

colnames(BMI)[c(1, 5, 7)] = c("SNP", "BETA", "PVAL")
head(BMI)

colnames(Height)[c(1, 5, 7)] = c("SNP", "BETA", "PVAL")
head(Height)

table(BMI$PVAL < 5E-8)

proportions(table(BMI$PVAL < 5E-8))

table(Height$PVAL < 5E-8)
proportions(table(Height$PVAL < 5E-8))

BMI_Height = cfdr_pleio$new()

BMI_Height$init_data(
trait1 = BMI, trait2 = Height,
trait_names = c("BMI", "Height"),
refdat = refdata_location(REF_DIR),
local_refdat_path = "E:/tmp/cond_conjFDR/BMI_height",
verbose = TRUE
)

BMI_Height$initialize_pruning_index(n_iter = 50, seed = 154226, verbose = TRUE)

BMI_Height$calculate_cond_fdr(fdr_trait = 1, verbose = TRUE)

BMI_Height$calculate_cond_fdr(fdr_trait = 2, verbose = TRUE)

BMI_Height_res = BMI_Height$get_trait_results()
BMI_Height_res

DATA_DIR = "E:/tmp/cond_conjFDR"
saveRDS(BMI_Height, file = paste(DATA_DIR, "BMI_Height_AnalysisObject.rds", sep = "/"))
saveRDS(BMI_Height_res, file = paste(DATA_DIR, "BMI_Height_ResultsTable.rds", sep = "/"))

table(BMI_Height_res$cfdr12 < 0.001)

summary(subset(BMI_Height_res, LOG10PVAL1 > -log10(5E-8))$cfdr12)

table(BMI_Height_res$conj_fdr < 0.01)

png("C:/Users/Administrator/Desktop/cfdr.pleio.png", width = 7, height = 7, units = 'in', res = 300)
plot(
LOG10PVAL2 ~ LOG10PVAL1,
data = BMI_Height_res[conj_fdr < 0.01],
xlab = "-log10(p) for BMI",
ylab = "-log10(p) for height",
xlim = c(0,50),
ylim = c(0, 50),
pch = 19,
col = gray(0.4, 0.25)
)
abline(h = -log10(5E-8), v = -log10(5E-8), lty = 3)
dev.off()

封装函数

1
2
3
4
5
6
7
8
9
10
# cond.conjFDR.v01.R
dat = cond_conjFDR(
GWASfile = c("E:/QTLMR/GWAS/MTAG/BMI_MTAG.txt", "E:/QTLMR/GWAS/MTAG/Height_MTAG.txt"),
trait_names = c("BMI", "Height"),
refdata_location = "E:/QTLMR/REF_DIR",
n_iter = 50,
seed = 154226,
save_name = "res",
save_path = "E:/tmp/cond_conjFDR"
)
1
2
3
4
5
6
7
8
# 参数说明
GWASfile : 包含两个或多个 GWAS 数据文件的路径,可由`format_dat()`转换获得,MATG 或 METAL 格式数据
trait_names : 性状名称,默认为 trait1, trait2
refdata_location : 字符串,遗传参考数据的文件夹路径。下载地址:https//zenodo.org/record/5750318/files/genref.zip,需解压使用
n_iter : 整数,迭代的次数,默认 50;为了最小化随机选择的影响;
seed : 整数;设置随机种子数,默认 12306
save_name : 字符串,输出文件名
save_path : 字符串,存储结果文件的目录

要点总结

cfdr.pleio是对原始 MATLAB 版 pleioFDR 的 R 语言重写,流程几乎一致。

它通过借用多效性信息,提升了检测关联变异体的能力,特别适合多性状基因组分析。

condFDR适合发现”次要表型”在”主要表型”条件下的关联;conjFDR适合发现同时影响两个性状的多效性变异体。

分析结果不仅能帮助发现新位点,还能对多性状遗传结构有更深入的理解。

相关文献

Age Cell, 10.1111/acel.14271

神经退行性疾病与表观遗传衰老和人类寿命的因果关系和共同遗传病因

Causal associations and shared genetic etiology of neurodegenerative diseases with epigenetic aging and human longevity

在本研究中,研究者旨在检测多种神经退行性疾病(AD、PD、LBD、ALS 和多发性硬化症 MS)与四种表观遗传时钟(GrimAge、PhenoAge、IEAA 和 HannumAge)以及神经退行性疾病和多变量长寿相关表型(父母寿命、健康寿命和异常长寿)之间的因果关系和遗传病因学重叠。为了实现这一目标,主要利用 MR 和条件/联合错误发现率(cond/conjFDR)方法,使用大规模全基因组关联研究(GWAS)数据集。此外,我们鉴定了包括神经退行性疾病、表观遗传老化和多变量长寿相关表型的共享分子表型的多效性遗传变体、基因和生物学途径。

Nature Communications, 10.1038/s41467-024-52121-y

抑郁症和皮质下脑结构体积的共享遗传机制

Investigating the shared genetic architecture between depression and subcortical volumes

识别抑郁症和皮质下体积的共享遗传基础,使用双变量因果混合模型 MiXeR 方法,揭示两者间多基因架构的重叠,并通过条件/联合错误发现 cond/conjFDR 分析来识别两者间的共享基因位点。此外,研究还分析了这些共享基因的功能注释及其在不同发育阶段的表达模式,探索它们与认知能力和行为症状的关联。

  • Title: 获取数量性状位点或药靶数据
  • Author: Xing Abao
  • Created at : 2025-05-02 23:33:04
  • Updated at : 2025-05-04 12:47:38
  • Link: https://bioinformatics.vip/2025/05/02/post-GWAS/获取数量性状位点或药靶数据/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments