eQTL结果使用coloc共定位分析

共定位分析（Colocalization Analysis）之前有提到过,它是一种统计方法，旨在确定GWAS识别的与复杂性状或疾病相关的遗传变异位点是否与xQTL共享相同的因果变异。这种分析帮助揭示GWAS信号的潜在...

共定位分析（Colocalization Analysis）之前有提到过,它是一种统计方法，旨在确定GWAS识别的与复杂性状或疾病相关的遗传变异位点是否与xQTL共享相同的因果变异。这种分析帮助揭示GWAS信号的潜在分子机制，特别是当GWAS位点位于非编码区时，通过与xQTL的共定位推断其可能通过调控基因表达或其他分子表型（如剪接、甲基化）影响性状.这里我们详细介绍一下非模式数据如何进行共定位分析

—

数据准备

首先需要做一个标准的leafcutter分析,得到每个样本的内含子表达量信息既然是共定位分析,必须要有的肯定是两个分析出来的结果,主要是beta和se重要一点

gwas数据

chr   ps      af      beta           se             p_adj---   -----   -----   ------------   ------------   ------------1     17949   0.107   1.507454e-02   1.530044e-02   3.253326e-011     17956   0.089   3.360965e-02   1.557065e-02   3.171018e-021     17961   0.222   2.423437e-03   9.881948e-03   8.064456e-011     17966   0.119   2.326773e-02   1.576214e-02   1.409847e-011     17978   0.179   3.636033e-03   1.577490e-02   8.178703e-011     17990   0.143   1.364413e-02   1.336977e-02   3.083351e-011     17991   0.105   2.168526e-02   1.649959e-02   1.897904e-011     17992   0.063   1.032086e-02   1.625761e-02   5.260398e-011     17994   0.251   2.367801e-04   9.529350e-03   9.801938e-01

都是很重要的信息,尤其是beta和se,eqtl数据也类似

eqtl数据

variant_id    phenotype_id   beta          se      p-value
-----------   ------------   -----------   ----------   ----------------------
Chr01:17961   Pca01G000020   -0.48926565   0.1056443    1.4508091617399284e-05
Chr01:17994   Pca01G000020   -0.5147455    0.10087119   2.3444868519549754e-06
Chr01:18008   Pca01G000020   0.70167667    0.23522985   0.003822462008448789
Chr01:18036   Pca01G000020   -0.565934     0.11335434   3.617333592398793e-06
Chr01:18057   Pca01G000020   -0.62238723   0.16974024   0.00044987542666411926
Chr01:18083   Pca01G000020   -0.5732085    0.18021214   0.002117700496080771
Chr01:18181   Pca01G000020   -0.50499636   0.13173887   0.0002566866447425307
Chr01:18529   Pca01G000020   0.48160124    0.14975485   0.0019017360125251617
Chr01:18532   Pca01G000020   -0.6260348    0.1406607    2.849532661018963e-05

这里的染色体啊位置啊什么的都不一样,但是没关系,后续代码里统一就行,比gwas多了一个基因,为了防止出现假阳,后面还需要算ld,所以需要一个完整的vcf文件,用plink转换成bed,这里不多赘述

—

计算coloc.abf

大规模数据计算共定位,计算压力非常大,所以我们可以用coloc.abf来进行一下粗筛,找到有可能的位点,首先是根据显著gwas画locus,然后逐locus寻找邻近的显著eqtl,计算步骤较为简单,不过最好是并行,因为locus可能很多,其中的maf直接从gwas的af计算就行

maf  <- pmin(g_sub$af, 1 - g_sub$af)
d1 <- list(snp = common, beta = g_sub$beta, varbeta = g_sub$se^2, MAF = maf, N = args$gwas_num, type = "quant")
d2 <- list(snp = common, beta = e_sub$beta, varbeta = e_sub$se^2, MAF = maf, N = args$eqtl_num, type = "quant")
res <- tryCatch(coloc.abf(d1, d2), error = function(e) NULL)

得到的结果文件格式如下

chr   start      end        Gene           nsnps   PP.H0.abf              PP.H1.abf              PP.H2.abf            PP.H3.abf             PP.H4.abf
---   --------   --------   ------------   -----   --------------------   --------------------   ------------------   -------------------   -----------------
1     37799402   37999402   Pca01G033150   290     0.120332595044591      0.0941437539823218     0.153938265297287    0.119923921024389     0.511661464651411
1     37799402   37999402   XLOC_001723    239     0.123001896931881      0.10760919761577       0.0887707111989779   0.0770582098108113    0.60355998444256
18    11253932   11453932   Pca18G008610   26      0.456580929669092      0.0152949835026247     0.0679093283250741   0.00181649319184851   0.458398265311361
18    14727056   14927056   Pca18G011790   368     2.15479581406245e-06   9.29730283794106e-07   0.469421762178083    0.202213162045844     0.328361991249976
18    14727056   14927056   Pca18G012390   172     0.328424949214389      0.0678871657998478     0.218899422792352    0.044907782276696     0.339880679916715
2     10716098   10916098   Pca02G013730   30      0.226134341636627      0.0215949529174234     0.138380333824193    0.0126135085251071    0.60127686309665
3     21575172   21775172   XLOC_011782    31      0.274839598463053      0.152821998305832      0.0285545462025935   0.0153490545979302    0.528434802430591
6     15668407   15868407   Pca06G014670   124     0.00396210175906738    0.000177465280308773   0.647975478740978    0.0287040893489277    0.319180864870719

其中最重要的是最后一列H4的后验概率,越高说明gwas和eqtl的共享这个因果变异的概率越高,如果只关注单因果变异假设,那么这个结果就是直接可以用的了但是阈值一般推荐>0.8,我这里的都不够显著,但是其实不一定,因为大多数情况共定位并非单因果,而是多个独立因果驱动,这个我们下一章在介绍

然后给大家提供一下并行计算的代码,这里的需要最终数据格式,和你的不一样就需要微调,思路可以借鉴一下

library(future.apply)
library(coloc)
plan(multicore, workers = 8)
screening_results <- future_lapply(locus_tasks, function(locus_df) {
    chr_i <- locus_df$chr[1]
    start <- locus_df$start[1]
    end   <- locus_df$end[1]
    g_local <- gwas_full[chr == chr_i & pos >= start & pos <= end]
    if (nrow(g_local) < 10) return(data.table())
    res_list <- list()
    for (target_gene in unique(locus_df$gene)) {
      e_local <- eqtl_full[gene == target_gene & chr == chr_i & pos >= start & pos <= end]
      common  <- intersect(g_local$snp_id, e_local$snp_id)
      if (length(common) < 5) next
      g_sub <- g_local[match(common, snp_id)]
      e_sub <- e_local[match(common, snp_id)]
      maf  <- pmin(g_sub$af, 1 - g_sub$af)
      d1 <- list(snp = common, beta = g_sub$beta, varbeta = g_sub$se^2, MAF = maf, N = args$gwas_num, type = "quant")
      d2 <- list(snp = common, beta = e_sub$beta, varbeta = e_sub$se^2, MAF = maf, N = args$eqtl_num, type = "quant")
      res <- tryCatch(coloc.abf(d1, d2), error = function(e) NULL)

      if(!is.null(res)) {
      res_list[[target_gene]] <- as.list(c(
        chr = chr_i,
        start = start,
        end   = end,
        Gene  = target_gene,
        res$summary))
      }
    }

    return(rbindlist(res_list))
},future.seed = TRUE, future.globals = FALSE, future.packages = "coloc")
screening_df <- rbindlist(screening_results)
candidates <- screening_df[PP.H4.abf > 0.3]

更多生信课程：

发表于 22小时前
阅读 ( 23 )
分类：GWAS

eQTL结果使用coloc共定位分析

你可能感兴趣的文章

相关问题

0 条评论

作家榜 »