专栏首页科技记者简单的snptest要不要学

简单的snptest要不要学

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就得到了,前路漫漫,依然有好多不懂的,持续学习!

本文分享自微信公众号 - 科技记者(kejijizhe),作者:zd200572

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-03-08

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 宏转录组学习笔记(一)

    前面提到,已经有两家公司通过宏转录组(Metatranscriptomics)测序检测肠道微生物,面向消费者提供检测服务。对宏转录组充满了好奇,有这样的比方说,...

    用户1075469
  • 使用R语言读取PUBMED存入MYSQL数据库

    最近,在科研狗网站看到了一个有趣的项目,使用R语言读取pubmed存入mysql数据库,之前报名没有报上,还是决心要跟着做一下,无奈R语言水平比较渣渣,只能复制...

    用户1075469
  • 分析粪便微生物移植后患者高通量单分子实时测序数据的工作流程

    有许多基于测序的方法来了解复杂宏基因组,从全样本鸟枪法测序到靶向扩增。虽然靶向方法在低测序深度提供有价值的数据,但它们受引物设计和PCR限制。全样本鸟枪法通常使...

    用户1075469
  • 偏爱MySQL,Nifty使用4个Web Server支撑5400万个用户网站

    【编者按】Nifty运营网站已经有很长一段时间,而在基于HTML5的WYSIWYG网页制作平台推出后,用户在该公司建立的网站已超过5400万个,同时其中大部分网...

    CSDN技术头条
  • 【学习】数据分析与数据挖掘类的职位必备技能

    大数据催生数据分析师 薪酬比同等级职位高20% 随着大数据在国内的发展,大数据相关人才却出现了供不应求的状况,大数据分析师更是被媒体称为“未来最具发...

    小莹莹
  • 猎聘:分析70万在线职位后,告诉你数据分析师前景

    大数据文摘
  • 云计算对BAT来说是一次飞跃

    云时代的来临,对BAT来说是一次飞跃。百度,阿里、腾讯借助其产业优势以云计算为核心技术发展在行业中占尽先机。那么,云时代只属于BAT吗? 不容置否,云时代的到来...

    静一
  • 你还在用Python做数据分析么?

    被大数据分析算法刷屏的各种推荐,刷个抖音,被频繁的推荐可能认识的人,其中就包括分手一年多的前女友;淘宝闲逛,推送的都是你妈妈搜索过的中老年大码女装;微博浑水,你...

    叫我龙总
  • Hibternate框架笔记

    Hibernate框架 配置 配置文件: 1 <?xml version="1.0" encoding="UTF-8"?> 2 <!DOCTYPE hibe...

    二十三年蝉
  • 【程序源代码】工作流Activiti不错的学习资料总结

    在第一家公司工作的时候主要任务就是开发OA系统,当然基本都是有工作流的支持,不过当时使用的工作流引擎是公司一些牛人开发的(据说是用一个开源的引擎修改的),名称叫...

    程序源代码

扫码关注云+社区

领取腾讯云代金券