用 FastQTL 进行 cis-eQTL 分析

生信菜鸟团

发布于 2020-05-20 16:35:28

6.3K1

发布于 2020-05-20 16:35:28

文章被收录于专栏：生信菜鸟团

上周给大家介绍了 Matrix eQTL 的用法，它利用高效的矩阵运算实现了在很短的时间内完成关联分析。在 eqtl 分析中，我们对每个基因都进行了大量检验，所以我们必须进行多重检验校正。最简单的方案就是用 Bonferroni 法校正 P 值。然而由于不同基因组区域的特异性以及不同位点的等位基因频率和 LD，Bonferroni 方法通常都会过于严格，导致许多假阴性。为了解决这个问题，一般的我们可以分析每种表型的数千个置换数据集，以得到这些关联的零分布。接着就可以得到这些观察值来自零分布的可能性，从而得到一个调整后的 P 值。

但实际上，在这种情况下执行置换我们需要一种能够在合理时间内进行大量计算的方法。尽管 Matrix eQTL 已在多个大规模研究中使用，它的一个主要缺点在于没有高效的内置置换方案，这会导致我们使用了非最佳的多重检验校正方法。

今天给大家介绍的 FastQTL，它用快速高效的置换方法（通过 beta 分布来进行置换检验）改进了 Matrix eQTL。

软件首页：http://fastqtl.sourceforge.net/

软件特点

•运行速度快•用法灵活，支持离散性和连续性协变量•用法简单•对集群运算友好

软件安装

FastQTL 支持 Linux 和 MacOS，官网可下载二进制文件和源码以及示例数据：http://fastqtl.sourceforge.net/

也可直接用 docker 安装：

docker pull biocontainers/fastqtl:v2.184dfsg-6b1-deb_cv1

输入文件格式

FastQTL 需要三个输入文件，分别为基因型，基因表达矩阵以及协变量。

基因型

•VCF 格式：

##fileformat=VCFv4.1
#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT UNR1 UNR2    UNR3    UNR4
chr7    123    SNP1    A    G    100    PASS    INFO    GT:DS    0/0:0.001 0/0:0.000    0/1:0.999    1/1:1.999
chr7    456    SNP2    T    C    100    PASS    INFO    GT:DS    0/0:0.001 0/0:0.000    0/1:1.100    0/0:0.100
chr7    789    SNP3    A    T    100    PASS    INFO    GT:DS    1/1:2.000 0/1:1.001    0/0:0.010    0/1:0.890

需要建立索引：

bgzip genotypes.vcf && tabix -p vcf genotypes.vcf.gz

表型

•BED 格式：

#Chr    start    end    ID    UNR1    UNR2    UNR3    UNR4
chr1    173863 173864    ENSG123    -0.50    0.82    -0.71    0.83
chr1    685395    685396    ENSG456    -1.13    1.18    -0.03    0.11
chr1    700304    700305    ENSG789    -1.18    1.32    -0.36    1.26

需要建立索引：

bgzip phenotypes.bed && tabix -p bed phenotypes.bed.gz

协变量

id    UNR1    UNR2    UNR3    UNR4
PC1    -0.02    0.14    0.16    -0.02
PC2    0.01    0.11    0.10    0.01
PC3    0.03    0.05    0.08    0.07
BIN    1    0    0    1

两种分析方法

Nominal pass in cis

要对示例数据集执行简单的分类检验：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --out nominals.default.txt.gz

运行成功，屏幕上会输出如下代码：

一些参数

关联分析就是用线性回归对基因型与基因表达的进行关联，类似于 R 语言中的 lm 函数。该模型假定表型呈正态分布。如果表型量化值不是正态分布，则可以通过以下方式转为正态分布：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --out nominals.quantile.txt.gz --normal

可用 --window 参数修改 cis-window 大小，默认为 1 Mb：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --out nominals.window2Mb.txt.gz --window 2e6

可用 --seed 参数设置随机种子，以便重复实验结果：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --out nominals.seed.txt.gz --seed 123456789

可用 --cov 参数添加协变量文件：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --out nominals.cov.txt.gz --cov covariates.txt.gz

排除/包括数据

要从分析中排除样品，突变，表型或协变量，可以使用以下选项之一：

1.排除样本：--exclude-samples file.exc2.排除突变：--exclude-sites file.exc3.排除表型：--exclude-phenotypes file.exc4.排除协变量：--exclude-covariates file.exc

例如，我们要在分析示例数据集时忽略 3 个样本，首先需要创建一个文本文件，其中包含要排除的样本的 ID，这里命名为 file.exc：

UNR1
UNR2
UNR3

然后在运行命令时加入 --exclude-samples file.exc 即可：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --out nominals.sub.txt.gz --exclude-samples file.exc

和这四个数据排除参数类似，我们还可以用以下选项指定要包含在分析中的样本，突变，表型和协变量：--include-samples，--include-sites，--include-phenotypes 以及 --include-covariates。

并行运算

方法一

把整个分析拆成 30 个子任务运行：

for j in $(seq 1 30); do (nohup fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.chunk$j.txt.gz --chunk $j 30 & );done

--chunk [int] [int] 参数用于指定需要处理的数据块。需要两个参数（X1，X2），表示用 FastQTL 处理 X2 块中的 X1 块。

在这个简单的示例中，我们仅将数据分为 30 个块。但在实际的全基因组分析中一般需要将数据分成 200、500 甚至 1,000 个块。

方法二

先用 --commands 参数将 VCF 文件拆为多个区域并生成相应命令，再进行并行运算。

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out results --commands 10 commands.10.txt

commands.10.txt 中的内容如下：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:17517460-20748406 --region 22:17517460-20748406
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:36424473-39052635 --region 22:36424473-39052635
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:24407642-30163001 --region 22:24407642-30163001
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:42017123-45704850 --region 22:42017123-45704850
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:20792146-22307210 --region 22:20792146-22307210
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:39052641-39915701 --region 22:39052641-39915701
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:30163352-36044443 --region 22:30163352-36044443
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:45809500-51222092 --region 22:45809500-51222092
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:22337213-24322661 --region 22:22337213-24322661
fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --out nominals.22:39928860-42017101 --region 22:39928860-42017101

并行运算：

cat commands.10.txt | while read c; do (nohup $c &);done

输出结果

当并行运算完后，我们可以把所有文件合并到一起以简化下游分析：

zcat nominals.chunk*.txt.gz | gzip -c > nominals.all.chunks.txt.gz

输出文件包含 5 列和 N 行，分别对应于 N 个表型-突变配对。

ENSG00000237438.1 indel:4D_22_16518157 -999303 0.542909 -0.0510761
ENSG00000237438.1 snp_22_16519289 -998171 0.432764 0.124424
ENSG00000237438.1 snp_22_16520141 -997319 0.0945196 -0.251906
ENSG00000237438.1 snp_22_16520948 -996512 0.102846 -0.274157
ENSG00000237438.1 snp_22_16523696 -993764 0.0676318 -0.324492
ENSG00000237438.1 snp_22_16523730 -993730 0.0674215 -0.206578
ENSG00000237438.1 snp_22_16523848 -993612 0.489458 0.062382
ENSG00000237438.1 indel:4D_22_16524572 -992888 0.299065 -0.192933
ENSG00000237438.1 snp_22_16526596 -990864 0.769977 0.36254
ENSG00000237438.1 snp_22_16527116 -990344 0.0541876 -0.340604

这五列分别表示：

1.是基因 ID2.SNP3.突变与基因之间的距离4.p 值5.斜率

我们也可用 --threshold 仅输出小于 p 值的结果：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --threshold 0.001 --out nominals.threshold.txt.gz

Permutation pass in cis

用 1,000 次排列对示例数据集进行基于 permutation 的分析：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --permute 1000 --out permutations.default.txt.gz

运行成功，屏幕上会输出如下代码：

可用 --permute 参数设置 permutation 的数量：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --permute 10000 --out permutations.10K.txt.gz

我们还可以在示例数据集上使用自适应的 permutation 数量。例如，运行 100 至 100,000 次置换：

fastQTL --vcf genotypes.vcf.gz --bed phenotypes.bed.gz --region 22:17000000-18000000 --permute 100 100000 --out permutations.adaptive.txt.gz

其他参数以及并行运算的方法与之前的默认分类检验方法一致。

输出结果

结果文件包含 10 列和 M 行，其中 M 为基因的数量：

ENSG00000237438.1 7966 0.991684 2771.43 370.402 snp_22_17542810 25350 5.48728e-13 -0.547704 0.000999001 1.88501e-09
ENSG00000177663.8 8165 1.00051 4005.51 391.217 snp_22_17497295 -68549 0.000167489 -0.33122 0.371628 0.35903
ENSG00000183307.3 8279 1.07083 2223.39 354.37 snp_22_17587975 -14282 3.11324e-08 -0.436689 0.000999001 7.27221e-05
ENSG00000069998.8 8466 1.01933 3658.49 384.641 snp_22_17639837 -6340 5.01155e-11 -1.1834 0.000999001 5.88777e-08
ENSG00000093072.10 8531 0.962558 3824.52 383.029 snp_22_17699299 -39826 9.028e-05 0.267016 0.251748 0.249387

这些列分别表示：

1.基因 ID2.该基因范围内的突变数量3.Beta 分布 shape1 参数的最大似然估计4.Beta 分布 shape2 参数的最大似然估计5.dummy6.针对该基因发现的最佳的突变 ID（即具有最小的p值）7.突变与基因之间的距离8.p 值9.斜率10.用 direct method 得到的 permutation p-value11.通过 beta 近似获得的 permutation p-value。建议在下游分析中使用该值。（具体算法参见原文：https://academic.oup.com/bioinformatics/article/32/10/1479/1742545 ）

检查结果

用 R 检查是否正确估计了 p 值：

R> d = read.table("permutations.all.chunks.txt.gz", hea=F, stringsAsFactors=F)
R> colnames(d) = c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "ppval", "bpval")
R> plot(d$ppval, d$bpval, xlab="Direct method", ylab="Beta approximation", main="Check plot")
R> abline(0, 1, col="red")

从上图可知，结果和期望一致。

校正 P 值

这部分主要涉及从严格到宽松的 3 种校正方法。首先，载入数据：

R> d = read.table("permutations.all.chunks.txt.gz", hea=F, stringsAsFactors=F)
R> colnames(d) = c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "ppval", "bpval")

Bonferroni correction

R> d$bonferroni = p.adjust(d$bpval, method="bonferroni")
# 取 <= 0.05
R> write.table(d[which(d$bonferroni <= 0.05), c(1,6)], "permutations.all.chunks.bonferroni.txt", quote=F, row.names=F, col.names=T)

Benjamini & Hochberg correction

R> d$bh = p.adjust(d$bpval, method="fdr")
# 取 FDR <= 0.1
R> write.table(d[which(d$bh <= 0.10), c(1,6)], "permutations.all.chunks.benjamini.txt", quote=F, row.names=F, col.names=T)

Storey & Tibshirani correction

R> library(qvalue) R> d$st = qvalue(d$bpval)$qvalues
# 取 FDR <= 0.1
R> write.table(d[which(d$st <= 0.10), c(1,6)], "permutations.all.chunks.storey.txt", quote=F, row.names=F, col.names=T)

本期的内容就到这里啦 ~

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-05-10，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度