PLINK2的Association analysis文档

快乐是最好

已于 2024-12-24 10:46:05 修改

阅读量973

点赞数 25

文章标签：笔记

于 2024-12-11 15:42:03 首次发布

本文链接：https://blog.csdn.net/weixin_49113540/article/details/144394969

版权

PLINK2的关联分析部分

PLINK2的Association analysis
实操

PLINK2的Association analysis

1.PLINK2源文档地址
 2.自己参考的一个很全的教程

基本知识

三个注意事项

现在使用主成分校正人群分层的影响，已经成为基本流程。通常使用–pca。
PLINK的关联分析不会校正小家庭结构（small-scale family structure），通常需要使用-king-cutoff进行过滤。
–glm计算的结果在次要等位基因计数（minor allele count, MAC）很小时不准，可以在–glm之前用–mac 20过滤。后面讲firth回归的部分我没看懂。

基本语法

#<>是必选参数，[]是可选参数
--glm ['zs'] ['omit-ref'] [{sex | no-x-sex}] ['log10'] ['pheno-ids']
      [{genotypic | hethom | dominant | recessive | hetonly}] ['interaction']
      ['hide-covar'] ['skip-invalid-pheno'] ['allow-no-covars']
      ['qt-residualize'] [{intercept | cc-residualize | firth-residualize}]
      ['single-prec-cc'] [{no-firth | firth-fallback | firth}]
      ['cols='<col set desc.>] ['local-covar='<file>] ['local-psam='<file>]
      ['local-pos-cols='<key col #s> | 'local-pvar='<file>] ['local-haps']
      ['local-omit-last' | 'local-cats='<cat. ct> | 'local-cats0='<cat. ct>]
  (aliases: --linear, --logistic)
--ci <size>
--condition <variant ID> [{dominant | recessive}] ['multiallelic']
--condition-list <variant ID file> [{dominant | recessive}] ['multiallelic']
--parameters <number(s)/range(s)...>
--tests ['all'] [number(s)/range(s)...]
--vif <max VIF>
--max-corr <val>

技术细节

对于双等位基因位点，一般来说会以MAF较低的碱基作为效应碱基进行回归，如果想要固定以ALT碱基进行回归，需要使用omit-ref。这样这个等位基因会被列在A1列（A1的意思是当前的等位基因current allele）（这一点非常重要！！！！！）
X染色体有两点特殊之处
- 通常，性别（来自.fam/psam文件）被设置为协变量。如果不想要这样，使用no-x-sex参数，或者使用sex参数显式指定性别作为协变量。注意不要同时用-sex指定psam文件中性别作为协变量和用–covar指定性别作为协变量，否则重复的列会让回归失败。
  注意.fam/psam文件中男=1，女=2。
对于每次回归，–glm输出名为plink2.<表型名称>.glm.<回归类型>[.zst]的报告
- 如果有多个连续变量，他们都没有缺失值或者有缺失值的是相同样本，可以用1次–glm实现一起分析。进行单个连续变量回归分析时，PLINK2仅比PLINK1快几百倍，但是PLINK2针对多个连续变量有特殊优化，速度还可以再快十倍。
- 连续变量仅支持线性回归。
- 分类变量有三种回归方式：
  - no-firth：即使用logistic回归
  - firth-fallback：默认选项，先使用logistic回归，如果不收敛，改用firth回归。
  - firth：直接使用firth回归，如果常见变异不太值得这么做，计算量太大，研究显示MAC<400时，firth回归的结果变得相关（relevant），啥意思？
  - 三种选项生成的文件扩展名分别为logistic/logistic.hybrid/firth
- 对于每个变异进行回归时，输出文件为每个变异输出多行信息，每行为自变量及协变量的信息，如变异本身、性别、PC1-10的β值、p值等信息，整个表格会很大。如果不需要非截距协变量的信息，可以使用hide-covar减小文件。也可以使用intercept增加截距项的信息（即当所有协变量都为零时的预测值）
- 提供了牺牲精度换取速度的选项，用不上，不记了。
- 2020年5月以后的版本，因为认为一般需要校正前几个主成分的影响如果不校正协变量的影响，需要用allow-no-covars，否则会报错。
- log10 选项将p值转化为-log(10)的形式输出
- –ci <size> 会在输出里显示每个变量的置信区间，95%置信区间就是–ci 0.95
每个回归执行前都会进行多重共线性检验，如果存在多重共线性，这个回归会被跳过，结果输出为NA。
- 主要是方差膨胀因子（variance inflation factor, VIF）计算，如果VIF>50，检验不通过，用–vif可以修改阈值
- 自变量predictors之间的相关性也被检验，如果有相关系数>0.999，检验不通过，用–max-corr可以修改阈值
- ERRCODE列报告了不通过多重共线性检验的变量，还包含了另一些错误的信息。下面列举报错信息：
  - 略，后面再补充
- VIF检验在某些情况下会过分严格，特别是在分类协变量存在很多水平（借用R levels的概念）时。可以考虑提高VIF阈值来解决这个问题。
在回归之前，会对协变量进行检查。包括只纳入协变量的多重共线性检查，和一个covariate-scale检查（用来识别需要使用–covar-variance-standardize的情况）。如果检查未通过，PLINK2会报错，如果想只跳过受影响的回归，用skip-invalid-pheno。
如果PLINK2发现有些样本或协变量与所有回归无关（例如常数协变量，全是0只有一个样本非0的协变量），在回归前这些会被去掉。可以用pheno-ids让PLINK2报告剩余的样本。
如果表型在所有剩余样本中都是常数，程序会报错，可以用skip-invalid-pheno跳过。
有时候可以用某些特定的snps作为协变量。感觉用不上，不翻译了。
可以使用局部协变量（local covariates），即在不同的变异之间不同的协变量，感觉用不上，不翻译。
‘genotypic’ modifier的作用，没看懂
dominant/recessive 标签指定了编码A1等位基因的模型。hetonly标签把基因型列替换成了0…1…0的显性偏差（dominance-deviation，指杂合的情况与两种纯合都不同？）列
[{genotypic | hethom | dominant | recessive | hetonly}] 标签用来确定基因效应的模型
- dominant：显性模型，AA/Aa被编码为1，aa为0
- recessive：隐性模型，aa为1，AA/Aa为0
- genotypic：没看懂
- hethom：没看懂
- hetonly：只把Aa编码为1，其余编码为0
- 啥都不加：加性模型，AA/Aa/aa分别编码为2/1/0
interaction 标签增加基因型与协变量之间的交互项。具体略。也可以增加特定的交互项，而不是所有协变量与基因型。
最后的提示：考虑到线性回归的速度比logistic/firth回归快很多，如果有大规模数据需要试探性研究，可以把分类变量编码成连续变量。但这个结果只有提示意义，毕竟是不一样的回归。
–pfilter <threshold>参数可以设定p值的阈值（1e-6/5e-8）。只会输出小于这个阈值的列。
–xchr-model <mode number> 设定了X染色体得编码模式
- 0 跳过X染色体
- 1 男性X染色体编码为0/1，女性0/1/2。PLINK1默认参数
- 2 男女X染色体都编码成0/1/2。PLINK2默认参数（是因为X染色体同源区段的存在所以男的也可以有2？）

重编码GWAS结果以提交给GWAS Catalog

略

多重校正

--adjust-file <filename> ['zs'] ['gc'] ['cols='<col set descriptor>]
              ['log10'] ['input-log10'] ['test='<test name, case-sensitive>]
--adjust-chr-field <field name search order>
--adjust-pos-field <field name search order>
--adjust-id-field <field name search order>
--adjust-ref-field <field name search order>
--adjust-alt-field <field name search order>
--adjust-provref-field <field name search order>
--adjust-a1-field <field name search order>
--adjust-test-field <field name search order>
--adjust-p-field <field name search order>
--adjust ['zs'] ['gc'] ['log10'] ['cols='<col set descriptor>]
--lambda <value>

对于未过滤的PLINK报告，–adjust-file会提供基本的多重校正的信息（bonferroni/FDR/…），以p值升序排列，储存为plink2.adjusted[.zst].

默认会从数据中推断膨胀系数（inflation factor）/基因组控制系数（genomic-control, GC）λ或称并将其记录在log中，也可以使用–lambda显式指定λ。
$\lambda_{GC}=\frac{median(\chi^2_{observed})}{0.456}$
- λ的原理：大多数变异与表型不相关，因此观察到的卡方分布与预期卡方分布之间没有偏差，0.456是df=1的卡方分布的中位数。（如果假设疾病与大量snps有很小的相关性，那就不适用这个方法了）
- 计算出λ后通过下面的公式就可以计算校正 $\chi^2$ 值
  $\chi^2_{corrected}=\frac{\chi^2_{observed}}{\lambda_{GC}}$
gc标签可以用GC矫正过的p值进行多重校正，否则默认使用未校正过的P值。
log10标签使得P值以负对数的形式输出，input-log10则是指定输入文件的P值为负对数

实操

#设置工作目录
WD = <work_directory>
#这里的主成分（principal component, PC）是先用--pca计算得到的，然后和临床资料合并到一个文件里去了
plink2 \
    --pfile ${WD}/Step67 \             #输入质控过的P文件
    --pheno ${WD}/PCsam.txt \          #指定表型和协变量的文件，如果在分开的文件里，也可以单独指定，略
    --pheno-name name1,name2,name3 \   #表型的列名，可以一次跑多个表型，列名用逗号分隔
    --covar-name Age,SEX,PC1-PC10 \    #协变量名，PC1-PC10指定前10个PC作协变量
    --covar-variance-standardize \     #协变量归一化
    --glm log10 hide-covar omit-ref \           #P值负对数处理；不显示协变量的信息
    --pfilter 1e-6 \                   #只展示P值小于这个数的变异
    --ci 0.95 \                        #展示β值的95%置信区间
    --out ${WD}/Test