簡單的snptest要不要學

  • 2020 年 3 月 18 日
  • 筆記

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, score 和 ml. 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就得到了，前路漫漫，依然有好多不懂的，持續學習！

VirMach 便宜 VPS

QNews

QNews