snptest学习笔记

在阅读肺癌PRS分型那篇文章时发现作者使用这个软件进行了多元逻辑回归,值得学习一下。在微信上能找到一篇教程来自陈文燕[bio生物信息]公众号的对性染色体进行关联分析(转载在了下一篇文章里),也可以学习下一些知识。

For the GWAS datasets, we calculated per-allele ORs and SEs using logistic regression analysis with the SNPTEST software (version 2.5.4) based on a probabilistic dosage model.

使用基于概率剂量模型的SNPTEST软件(版本2.5.4)使用Logistic回归分析计算每个等位基因的ORs和Se.

根据这句话,尝试做点学习,记录点笔记。

Logit笔记

先补充学习了Logit的相关知识,手记笔记:

Logistic回归适用于二值响应型变量(0和1),模型假设响应亦是Y服从二项分布。通过一系列连续和/类别型 预测变量来预测二值型结果时使用。《R语言实战》(一本初看没有多少头绪,慢慢才有些感觉的书。最近还读了《R语言轻松入门与提高/达人迷》觉得这本书浅显易懂,虽然有点老,推荐阅读。

软件下载和准备

下载地址地这个网站 https://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.html

#我是黑苹果,下载的mac版,虽然两天死机一次,仍在坚持用,因为实在受不了bug10的咖喱风
 wget -c http://www.well.ox.ac.uk/~gav/resources/snptest_v2.5.4-beta3_MacOSX_x86_64.tgz
 #ubuntu,动态的比较推荐,适应更多平台
 wget -c http://www.well.ox.ac.uk/~gav/resources/snptest_v2.5.4-beta3_linux_x86_64_dynamic.tgz
 #Cent
 wget -c http://www.well.ox.ac.uk/~gav/resources/snptest_v2.5.4-beta3_CentOS6.6_x86_64_dynamic.tgz
 ##解压安装
 tar zxvf snptest_v2.5.4-beta3_MacOSX_x86_64.tgz
 #测试
snptest_v2.5.4-beta3_MacOSX_x86_64/snptest_v2.5.4-beta3 -help
Welcome to SNPTEST v2.5.3 (revision 8a0c171055ff6ac4c69079501f52ac93c563ed20)
© University of Oxford 2008-2015
https://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.html
Read LICENCE file for conditions of use.

Usage: snptest_v2.5.3 <options>
 #好多信息,此处省略n百字符

使用

1.先确定哪个选项完成这个功能

官网的文档中,有几个功能作了介绍,第一个Computing summary statistics,应该是汇总统计,主要功能介绍为:计算每个基因型的计数,频率,SNP缺失率和OR值。显然不是这个功能,顺便用搜索logistic 确认了下,没有。第二个是Frequentist Association Tests,频率关联测试,没找到相关内容。第三个,Bayesian Tests (Bayes Factors),贝叶斯测试,找到了相关内容。后面两个选项也是有相关内容的,但是在这里又看到了这句话的引文,来自2007年的一篇nature,A new multipoint method for genome-wide association studies by imputation of genotypes,发现这篇文章中涉及的两个方法是 frequentist和bayesian,所以断定应该是后者了。

简单看下这篇引文,会不会有什么收获,基本上是关于detect causal variants that have not been directly genotyped的,翻译一下就是检测尚未直接进行基因分型的因果变异。所推断的基因型有时(适当地)是不确定的,所以我们发现有必要开发关联性检验(包括频度检验和贝叶斯检验),以考虑到我们推定的基因型的不确定性。贝叶斯因子在某种程度上类似于频率P值,它们的使用开始出现在文献中,作为经典关联检验的一种更强大和更容易解释的选择。

还是来自引文的内容:使用贝叶斯因子比频率测试统计量或P值有几个优点。正确解释P值需要了解所用测试的威力。非正式地,小的P值可能在NULL下偶然出现或来自真正的关联。在没有对可能的效应大小的研究的威力了解的情况下,评估其中的哪一种可能是困难的(例如,对于动力不足的研究,我们预计最重要的P值会偶然出现)。贝叶斯因子的计算,就像幂计算一样,需要关于效应大小的假设,但贝叶斯因子本身具有自然的解释,作为根据数据改变我们先前的关联概率的因子。贝叶斯因子可以在给定的SNP下通过不同的关联模型自然地组合。例如,我们可以跨加性模型、显性模型、隐性模型和一般模型求贝叶斯因子的平均值,以避免必须指定在一个位点使用的单个模型。可以使用类似的思想来组合区域内跨SNP的贝叶斯因子。根据最近关于贝叶斯方法获得的能力的证据,我们重点研究了基于贝叶斯因子的测试统计,并在使用的两组测试统计中对方法进行了比较,以便将结果集中在每种方法预测因果变量的能力上,而不是集中在不同测试统计数据的能力差异上。我们分别使用非共轭和共轭先验(补充方法)对2型糖尿病数据集和模拟数据集进行分析,以反映我们对适合这些数据集的遗传效应大小的信念。

然后,这里有带了个参考文献,顺藤摸瓜,找到了,一个大牛的教程,但是是nature!!! 厉害吧。简单扫了一眼,好像看不懂,这充分告诉我们,学好数学是多么地重要,回归今天的主题吧!

2.用示例数据学习下Bayesian Tests (Bayes Factors)

-bayesian选项,此选项控制您希望在每个SNP上测试与模型没有关联的SNP。 这五个不同的模型被编码为1 =Additive,2 =显性,3 =隐性,4 =常规和5 =杂合子。 使用此选项时,输出文件将为每个测试包含一列,其中包含该测试的log10贝叶斯因子以及模型参数(β值)及其标准误差的后验均值。 SNPTEST将allele_A编码为0,将allele_B编码为1,这定义了beta的含义以及se的含义。 例如,当使用加性模型时,β估计对数几率的增加,这可以归因于allele_B的每个副本。 贝叶斯因子将始终以每个SNP计算。

-method选项还用于控制贝叶斯模型拟合的方式,但并非所有选项都有效。

  • 如果表型是二进制,那么有效的选项只有 threshold, expected, scoreml. score选项使用单个牛顿-拉夫森迭代来估计后验模式,而ml选项使用多次迭代。
  • 如果是数量表型,那么有效的选项只有threshold和expected
贝叶斯模型的先验

下表说明了所使用的Logistic回归的线性预测值、模型参数使用的先验的形式、SNPTEST中使用的默认先验以及可用于更改先验的命令行选项。

Model Linear Predictor Priors Default Coding Command line option
Additive log(pi/(1-pi)) = µ + ßGi µ~N(a0, a12) ß~N(b0, b12) a0=0, a1=1 b0=0, b1=0.2 Gi is the additive coding of the SNP i.e. AA -> 0, AB ->1, BB -> 2. -prior_add a0 a1 b0 b1
Dominant log(pi/(1-pi)) = µ + ßDi µ~N(a0, a12) ß~N(b0, b12) a0=0, a1=1 b0=0, b1=0.5 Di is the dominant coding of the SNP i.e. AA -> 0, AB -> 1, BB -> 1. -prior_dom a0 a1 b0 b1
Recessive log(pi/(1-pi)) = µ + ßRi µ~N(a0, a12) ß~N(b0, b12) a0=0, a1=1 b0=0, b1=0.5 Ri is the recessive coding of the SNP i.e. AA -> 0, AB -> 0, BB -> 1. -prior_rec a0 a1 b0 b1
General log(pi/(1-pi)) = µ + ßGi + qHi µ~N(a0, a12) ß~N(b0, b12) q~N(c0, c12) a0=0, a1=1 b0=0, b1=0.2 c0=0, c1=0.5 Gi is the additive coding of the SNP i.e. AA -> 0, AB ->1, BB -> 2. Hi is the heterozygote coding of the SNP i.e. AA -> 0, AB ->1, BB -> 0. -prior_gen a0 a1 b0 b1 c0 c1
Heterozygote log(pi/(1-pi)) = µ + ßHi µ~N(a0, a12) ß~N(b0, b12) a0=0, a1=1 b0=0, b1=0.5 Hi is the heterozygote coding of the SNP i.e. AA -> 0, AB ->1, BB -> 0. -prior_het a0 a1 b0 b1
T分布先验

在SNPTEST v2中,有一个新的选项来指定在遗传效应上使用t分布先验。t分布的较粗尾部允许更灵活地指定关于遗传效应大小的信念。此选项由以下两个选项控制。

-t_prior 详细说明了t-分布先验在遗传效应上的应用。该选项有效地修改了上表中描述的先验,即t-分布的均值和方差由上表中给出的选项指定,但是正态分布被t-分布代替。注:t分布仅用于遗传效应,即上述模型中的参数?和q。例如,-Bayesian add-t_previous将指定线性预测器log(pi/(1-pi))=µ+?Gi,并且先验将是µ~N(a0,a12)和?t(b0,b12,df=3)。
**-t_df ** t分布的自由度参数。默认值为3。当此参数设置得非常大时,先验收敛到正态分布先验。
例子-Bayesian病例-对照test

使用Bayesian addictive 模型对二值型表型bin1计算,使用默认参数

snptest_v2.5.4-beta3_MacOSX_x86_64/snptest_v2.5.4-beta3  \
-data snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort1.gen snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort1.sample snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort2.gen snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort2.sample \
-o snptest_v2.5.4-beta3_MacOSX_x86_64//example/ex.out \
-bayesian 1 \
-method score \
-pheno bin1
#输出了一些信息
#这里显示2.5.3版本了
Welcome to SNPTEST v2.5.3 (revision 8a0c171055ff6ac4c69079501f52ac93c563ed20)
© University of Oxford 2008-2015
https://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.html
Read LICENCE file for conditions of use.

==============

Data Files :
 -gen files : snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort1.gen snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort2.gen
 -sample files : snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort1.sample snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort2.sample

Tests :
 -bayesian : 1
 -method score

Inspecting data (this may take some time)...
Sample and exclusions summary :#样本量500-500
 - Number of individuals in : (cohort 1)   (cohort 2)
                              500          500


Reading sample files :
Summary of covariates and phenotypes
 # discrete variables : 4 离散变量
  cov1 : type = D (Discrete covariate)
  cov2 : type = D (Discrete covariate)
  sex : type = D (Discrete covariate)
  bin3 : type = D (Discrete covariate)
 # continuous variables : 2 连续变量
  cov3 : type = C (Continuous covariate)
  cov4 : type = C (Continuous covariate)
 # phenotypes : 4 表型
  pheno1 : type = P (Continuous phenotype)
  pheno2 : type = P (Continuous phenotype)
  bin1 : type = B (Binary phenotype)
  bin2 : type = B (Binary phenotype)
Covariate summary :
  (no covariates in use.)
Phenotype summary :
  bin1    : missing  levels
            2        0(499) 1(499)

You have specified the following model: 模型
  bin1 ~ (baseline) + (genotype)

Phenotype being used : bin1

Data Summaries :
 -number of SNPs = (unknown)

Data with missing genotype data threshold and exclusion list applied :
 snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort1.gen : 500
 snptest_v2.5.4-beta3_MacOSX_x86_64//example/cohort2.gen : 500

--------------------------------------------------------------------------

SinglePhenotypeTest #检验
--------------------------------------------------------------------------

Analyzing Data :
 scanning... read chunk [1 of (unknown)]... done.
 scanning... read chunk [2 of (unknown)]... done.
 scanning... no more data.

finito

来看看结果文件,都有什么内容:

less snptest_v2.5.4-beta3_MacOSX_x86_64/example/ex.out
#有点不好分清列,还是用excel,数据-分列-字符-空格

看看表头,好多信息,分了n行才列出来,第一行分别是染色体号,位置,等位基因,最大分型率

chromosome position alleleA alleleB index average_maximum_posterior_call info cohort_1_AA
INSERTION_1 RSID_1 NA 0 A AGTGCTA 1 0.992234

第二行是两个队列的相关信息。

cohort_1_AB cohort_1_BB cohort_1_NULL cohort_2_AA cohort_2_AB cohort_2_BB cohort_2_NULL
0.970846 0 0 0 499 176.815 250.974

第三行是总体基因型汇总信息

all_AA all_AB all_BB all_NULL all_total cases_AA cases_AB cases_BB
71.2104 0 176.815 250.974 71.2104 499 998 0

第四行是样本对照统计信息

cases_NULL cases_total controls_AA controls_AB controls_BB controls_NULL controls_total
0 499 499 176.815 250.974 71.2104 0

第五行是接着上面的,还有数据缺失情况和OR值信息。

all_maf cases_maf controls_maf missing_data_proportion het_OR het_OR_lower
499 0.394184 0 0.394184 0.25 NA

第六行还是上一行的OR值信息

het_OR_upper hom_OR hom_OR_lower hom_OR_upper all_OR all_OR_lower
NA NA NA NA NA NA

第七行是本次分析得到的目标参数,log10(bf)贝叶斯因子?,beta值,标准差。

all_OR_upper bayesian_add_log10_bf bayesian_add_beta_1 bayesian_add_se_1 comment
NA -0.0174105 2.02E-17 0.192135 NA

大概到这里,这句话要获得的信息OR值和se就得到了,前路漫漫,依然有好多不懂的,持续学习!

发表评论