SAIGE用戶手冊筆記1

群體遺傳學(xué)領(lǐng)域，看起來華人（中國人）至少是占據(jù)前線的，至少很多軟件的一作就是中國人的名字，至少活主要是我們干的。當(dāng)然，也有很多華人（國人）大牛開發(fā)了史詩級別的軟件，向勤勞的華夏人致敬！這個軟件也是如此。

總覽

簡介

SAIGE是與Rcpp一起開發(fā)的R包，用于大規(guī)模數(shù)據(jù)集和生物庫中的全基因組關(guān)聯(lián)測試。方法

  • 基于廣義混合模型的樣本相關(guān)性
  • 允許使用全遺傳關(guān)系矩陣或稀疏遺傳關(guān)系矩陣（GRM）進行模型擬合
  • 適用于數(shù)量和二元分類性狀
  • 處理二元性狀的病例-對照不平衡
  • 對于大型數(shù)據(jù)集，計算效率高
  • 執(zhí)行單變異關(guān)聯(lián)測試
  • 通過 Firth 偏倚減少邏輯回歸提供效應(yīng)大小估計
  • 執(zhí)行條件關(guān)聯(lián)分析

SAIGE-GENE（現(xiàn)在稱為SAIGE-GENE+）是R包中的新方法擴展，用于基于集合的罕見變異分析。

  • 執(zhí)行 BURDEN、SKAT 和 SKAT-O 分析
  • 允許對多個次要等位基因頻率閥值和功能注釋進行分析
  • 允許在基于集合的分析中指定標(biāo)記的權(quán)重
  • 執(zhí)行條件分析以識別獨立于近 GWAS 信號的關(guān)聯(lián)

該軟件包采用以下格式的基因型文件輸入

  • PLINK （bed， bim， fam）， BGEN， VCF， BCF， SAV

引用

SAIGE:

  • Zhou W, Nielsen JB, Fritsche LG, Dey R, Gabrielsen ME, Wolford BN, LeFaive J, VandeHaar P, Gagliano SA, Gifford A, Bastarache LA, Wei WQ, Denny JC, Lin M, Hveem K, Kang HM, Abecasis GR, Willer CJ, Lee S. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nat Genet. 2018 Sep;50(9):1335-1341. doi: 10.1038/s41588-018-0184-y. Epub 2018 Aug 13. PMID: 30104761; PMCID: PMC6119127.

SAIGE-GENE:

  • Zhou W, Zhao Z, Nielsen JB, Fritsche LG, LeFaive J, Gagliano Taliun SA, Bi W, Gabrielsen ME, Daly MJ, Neale BM, Hveem K, Abecasis GR, Willer CJ, Lee S. Scalable generalized linear mixed model for region-based association tests in large biobanks and cohorts. Nat Genet. 2020 Jun;52(6):634-639. doi: 10.1038/s41588-020-0621-6. Epub 2020 May 18. PMID: 32424355; PMCID: PMC7871731.

SAIGE-GENE+:

  • Wei Zhou, Wenjian Bi, Zhangchen Zhao*, Kushal K. Dey, Karthik A. Jagadeesh, Konrad J. Karczewski, Mark J. Daly, Benjamin M. Neale, Seunggeun Lee. Set-based rare variant association tests for biobank scale sequencing data sets. medRxiv 2021.07.12.21260400; doi: https://doi.org/10.1101/2021.07.12.21260400

許可證

SAIGE在GPL許可證下分發(fā)。

聯(lián)系

如果您對SAIGE有任何疑問，請聯(lián)系 saige.genetics@gmail.com

框架圖

安裝

如何安裝和運行 SAIGE 和 SAIGE-GENE

安裝 SAIGE/SAIGE-GENE（當(dāng)前版本 1.0.0（更新于 2022 年 3 月 15 日））

依賴項列表：

  • R >= 3.6.1， gcc >= 5.4.0， cmake 3.14.1， cget， savvy

  • R 包： Rcpp （>= 1.0.7）， RcppArmadillo， RcppParallel， data.table， SPAtest （== 3.1.2）， RcppEigen， Matrix， methods， BH， optparse， SKAT， qlcMatrix， RhpcBLASctl， roxygen2， rversions， devtools， dplyr， dbplyr

  • Rscript ./extdata/install_packages.R 可用于安裝 R 包

  • 兩個用于讀取 VCF 文件的庫。將在 SAIGE 安裝過程中自動安裝

  • cget，savvy

日志：

v1.0.1（2022 年 3 月 18 日）：

修復(fù)錯誤： –AlleleOrder bgen 輸入單變異關(guān)聯(lián)分析輸出的標(biāo)題

改進：在 CCT 中將 PI 更改為M_PI.cpp在 Makevar 中添加標(biāo)志以安裝在 Centos 8 上

v1.0.0（2022 年 3 月 15 日）：第一個穩(wěn)定版本

v0.99.3（2022 年 3 月 12 日）：對于基于集合的分析，所有與 –condition 一起使用的標(biāo)記 ID 和組文件中都需要使用 chr:pos:ref:alt 格式

使用源代碼 https://github.com/saigegit/SAIGE 安裝 SAIGE

安裝依賴項

R >= 3.6.1， gcc >= 5.4.0， cmake 3.14.1，
從 github 下載 SAIGE

  repo_src_url=https://github.com/saigegit/SAIGE
  git clone --depth 1 -b $src_branch $repo_src_url

安裝依賴項：R 包

Rscript ./SAIGE/extdata/install_packages.R

安裝 SAIGE R 軟件包

將 SAIGE 安裝到存儲所有 R 庫的根目錄

–library=path_to_final_SAIGE_library 可用于指定安裝 SAIGE 的目錄
 R CMD INSTALL --library=path_to_final_SAIGE_library SAIGE

運行 SAIGE

如果未在根 R 庫路徑中安裝 SAIGE，請更改
library(SAIGE)
自

library(SAIGE, lib.loc="path_to_final_SAIGE_library")
在以下腳本中

 SAIGE/extdata/step1_fitNULLGLMM.R
 SAIGE/extdata/step2_SPAtests.R
 SAIGE/extdata/createSparseGRM.R

使用 conda 環(huán)境安裝 SAIGE

  1. 下載并安裝miniconda

  2. 使用

    • conda 環(huán)境文件位于 SAIGE 文件夾中：./conda_env/environment-RSAIGE.yml

    • 下載環(huán)境-RSAIGE.yml 后，運行以下命令

        conda env create -f environment-RSAIGE.yml

  3. 激活 CONDA 環(huán)境 RSAIGE

        conda activate RSAIGE
    FLAGPATH=`which python | sed 's|/bin/python$||'`
    export LDFLAGS="-L${FLAGPATH}/lib"
    export CPPFLAGS="-I${FLAGPATH}/include"

    請確保使用導(dǎo)出（最后兩個命令行）設(shè)置LDFLAGS和CPPFLAGS，以便在編譯SAIGE源代碼時可以正確鏈接庫。

注意：以下是創(chuàng)建 conda 環(huán)境文件的步驟

  1. 從源代碼安裝 SAIGE。

        src_branch=main
    repo_src_url=https://github.com/saigegit/SAIGE
    git clone --depth 1 -b $src_branch $repo_src_url
    Rscript ./SAIGE/extdata/install_packages.R
    R CMD INSTALL --library=path_to_final_SAIGE_library SAIGE

    當(dāng)在 R 中調(diào)用 SAIGE 時，設(shè)置 lib.loc=path_to_final_SAIGE_library

        library(SAIGE, lib.loc=path_to_final_SAIGE_library)

使用 Docker 映像運行 SAIGE

感謝 Juha Karjalainen 分享 Dockerfile。

Dockerfile 可以在 SAIGE 文件夾中找到：./docker/Dockerfile

可以拉取映像
docker pull wzhou88/saige:1.0.0
可以調(diào)用函數(shù)

docker run --rm -it  -v /data/haod/:/home wzhou88/saige:1.0.0
cd /home
step1_fitNULLGLMM.R --help
step2_SPAtests.R --help
createSparseGRM.R --help

### 從 bioconda 安裝 SAIGE （不是當(dāng)前版本）

警告：請不要將此 BIOCONDA 版本用于 BGEN 輸入。我們正在研究這個問題。

要從 conda 安裝 saige，只需使用最新版本的 R 和 saige 創(chuàng)建環(huán)境：

conda create -n saige -c conda-forge -c bioconda "r-base>=4.0" r-saige
conda activate saige

有關(guān) r-saige conda 軟件包和可用版本的更多信息，請參閱問題 #272。

單變異關(guān)聯(lián)分析

SAIGE采取兩個步驟來執(zhí)行單變異關(guān)聯(lián)分析

  • 我們建議對 MAC >= 20 的變體進行單變體關(guān)聯(lián)分析
  • 對于罕見的變異關(guān)聯(lián)，請使用 SAIGE-GENE+進行基于集合的關(guān)聯(lián)分析

步驟 1：擬合null(空)邏輯/線性混合模型

  • 對于二元性狀，將擬合一個空邏輯混合模型（–traitType=binary）。

  • 對于定量性狀，將擬合空線性混合模型（-traitType=quantitative），并且需要逆歸一化（–invNormalize=TRUE)

  • 默認情況下，協(xié)變量作為補償包含在模型中，請使用-isCovariateOffset=FALSE停用此功能。這將花費更多的計算時間與 -isCovariateOffset=FALSE

  • 示例腳本位于 extdata 文件夾中

#go to the folder
cd extdata
#check the help info for step 1
Rscript step1_fitNULLGLMM.R --help

使用完整 GRM 擬合空模型

對于二元性狀，將擬合一個空邏輯混合模型（–TRAITTYPE=BINARY）。

  • 使用 -qCovarColList 列出分類協(xié)變量。將為不同的級別創(chuàng)建虛擬變量（獨熱編碼）。
  • 使用 -nThreads表示要使用的 CPU 數(shù)量
    Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly_22chr  \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_binary \
        --covarColList=x1,x2,a9,a10 \
        --qCovarColList=a9,a10  \
        --sampleIDColinphenoFile=IID \
        --traitType=binary        \
        --outputPrefix=./output/example_binary \
        --nThreads=24   \
        --IsOverwriteVarianceRatioFile=TRUE

對于定量性狀，將擬合空線性混合模型（–TRAITTYPE=定量），并且需要逆歸一化（–INVNORMALIZE=TRUE)

    Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly_22chr  \
        --useSparseGRMtoFitNULL=FALSE    \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_quantitative \
        --covarColList=x1,x2,a9,a10 \
        --qCovarColList=a9,a10  \
        --sampleIDColinphenoFile=IID \
        --invNormalize=TRUE     \
        --traitType=quantitative        \
        --outputPrefix=./output/example_quantitative \
        --nThreads=24   \
        --IsOverwriteVarianceRatioFile=TRUE

使用稀疏 GRM 擬合空模型（將僅使用一個 CPU）

參考GWASLab的文章，這個方法不推薦，因為沒有LOCO。

  • Use –useSparseGRMtoFitNULL=TRUE

  • 對包含稀疏 GRM 的文件使用 –稀疏 GRMFile

  • 對包含稀疏 GRM 中樣本 ID 的文件使用 –稀疏 GRM 樣本 ID

  • 如果未使用 –plinkFile 指定 plink 文件，則不會估計方差比，請確保在步驟 2 中指定 –sparseSigmaFile（由步驟 1 輸出）以進行關(guān)聯(lián)測試

  • 要估計方差比并在步驟 2 中使用，請指定包含將用于使用 –plinkFile 進行方差比估計的標(biāo)記的 plink 文件。例如 –plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly_22chr 和 –skipVarianceRatioEstimation=FALSE 注意：這將加快步驟 2 的速度

  • 要僅包含樣本子集以擬合空模型，請使用 –SampleIDIncludeFile

  • 將僅使用一個 CPU，并且不會應(yīng)用 LOCO

    Rscript step1_fitNULLGLMM.R     \
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   \
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     \
        --useSparseGRMtoFitNULL=TRUE    \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_binary \
        --covarColList=x1,x2,a9,a10 \
        --qCovarColList=a9,a10  \
        --sampleIDColinphenoFile=IID \
        --traitType=binary        \
        --outputPrefix=./output/example_binary_sparseGRM

輸入文件

  1. （必填）表型文件（包含協(xié)變量（如果有），如性別和年齡）文件可以是空格，也可以是用標(biāo)題以制表符分隔的。該文件必須包含一列用于樣本 ID，一列用于表型。它可能包含協(xié)變量列。
less -S ./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt
輸入
  1. （可選）用于構(gòu)建完整遺傳關(guān)系矩陣（GRM）和估計方差比的基因型
    文件 SAIGE 采用基因型的 PLINK 二進制文件，并假設(shè)文件前綴與 .bed， .bim和 .fam的文件前綴相同。
     ./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly.bed
     ./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly.bim
     ./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly.fam
  1. （可選）稀疏 GRM 文件和稀疏 GRM 的示例 ID 文件。有關(guān)更多詳細信息，請參閱步驟 0。
   ./output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx
./output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt
   ##the sparse matrix can be read and viewed using the R library Matrix
   library(Matrix)
   sparseGRM = readMM("./output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx")

輸出文件

  1. 模型文件（步驟 2 的輸入）)

     ./output/example_quantitative.rda
 #load the model file in R
 load("./output/example_quantitative.rda")
 names(modglmm)
 modglmm$theta
    輸出

  2. 方差比文件（如果在步驟 1 中估計了方差比，則將在步驟 2 中輸入)

     less -S ./output/example_quantitative.varianceRatio.txt
 1.21860521240473

  3. 隨機選擇的標(biāo)記子集的關(guān)聯(lián)結(jié)果文件（如果在步驟 1 中估計了方差比，則將生成此文件。它是一個中間文件，后續(xù)步驟不需要它。)

     less -S ./output/example_quantitative_30markers.SAIGE.results.txt
markerIndex p.value p.value.NA var1 var2 Tv1 N AC AF
69970 0.214147816412411 0.173552665987772 1.14375062050603 0.953027027027038       
1.32853003974658 1000 74 0.037
96109 0.953143629473153 0.949713298740539 1.07324590784922 0.931653266331636     
-0.0608734702523558 1000 199 0.0995

步驟 2：執(zhí)行單變體關(guān)聯(lián)測試

  • 對于二元性狀，鞍點近似用于解釋病例--對照不平衡。
  • 可以使用要分析的遺傳變異的劑量/基因型的文件格式：PLINK，VCF，BGEN，SAV
  • 可以在步驟 2 中執(zhí)行基于條件分析的匯總統(tǒng)計信息（–condition）
  • 查詢和測試標(biāo)記子集
    • 變體 ID （chr:pos_ref/alt） 和染色體位置范圍（chr 開始 結(jié)束）都可以為 BGEN 輸入指定 （-idstoIncludeFile， -rangestoIncludeFile）
  • -markers_per_chunk可用于指定要測試的標(biāo)記數(shù)并將其輸出為一個塊。默認值 = 10000。請注意，小數(shù)字可能會減慢運行速度。要求此數(shù)字為 > = 1000。
  • 如果 LOCO=TRUE（默認情況下），則必須指定–chrom，因此只有基因型/劑量文件應(yīng)僅包含一條染色體
  • 對于 VCF/BCF/SAV 輸入，-vcfField=DS 用于劑量，-vcfField=GT 用于分析基因型
  • 默認情況下，缺失的基因型/劑量將被歸因為最佳猜測的 gentoypes/劑量（如 round（2xfreq），–impute_method=best_guess）。請注意，目前不支持刪除缺少基因型/劑量的樣本
  • –sampleFile 用于為 bgen 文件指定具有示例 ID 的文件。請不要在文件中包含表頭。
#check the help info for step 2
Rscript step2_SPAtests.R --help
  • 對于二元性狀，使用 –is_output_moreDetails=TRUE 輸出病例和對照組中的雜合子和純合子計數(shù)

如果在步驟 1 中估計了方差比 （–varianceRatioFile=）

  • Using –bgenFile， –bgenFileIndex， –AlleleOrder， –sampleFile for bgen input
Rscript step2_SPAtests.R        \
        --bgenFile=./input/genotype_100markers.bgen    \
        --bgenFileIndex=./input/genotype_100markers.bgen.bgi \
        --sampleFile=./input/samplelist.txt \
        --AlleleOrder=ref-first \
        --SAIGEOutputFile=./output/genotype_100markers_marker_bgen_Firth.txt    \
        --chrom=1   \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary.rda \
        --varianceRatioFile=./output/example_binary.varianceRatio.txt  \
        --is_Firth_beta=TRUE    \
        --pCutoffforFirth=0.1
  • 使用 –bedFile， –bimFile， –famFile， –AlleleOrder進行 PLINK 輸入
Rscript step2_SPAtests.R        \
        --bedFile=./input/genotype_100markers.bed       \
        --bimFile=./input/genotype_100markers.bim       \
        --famFile=./input/genotype_100markers.fam       \
        --AlleleOrder=alt-first \
        --SAIGEOutputFile=./output/genotype_100markers_marker_plink.txt \
        --chrom=1       \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary.rda \
        --varianceRatioFile=./output/example_binary.varianceRatio.txt   \
        --is_output_moreDetails=TRUE
  • 使用 –vcfFile、 –vcfFileIndex、 –vcfField 進行 VCF、BCF 和 SAV 輸入
Rscript step2_SPAtests.R        \
        --vcfFile=./input/genotype_100markers.vcf.gz    \
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     \
        --vcfField=GT   \
        --SAIGEOutputFile=./output/genotype_100markers_marker_vcf.txt \
        --chrom=1       \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary.rda \
        --varianceRatioFile=./output/example_binary.varianceRatio.txt   \
        --is_output_moreDetails=TRUE

如果在步驟 1 中未估計方差比，因為使用了稀疏 GRM 來擬合空模型（與步驟 1 中相同，–sparseGRMFile and –sparseGRMSampleIDFile 用于指定稀疏 GRM 和樣本 ID） 

Rscript step2_SPAtests.R        \
        --vcfFile=./input/genotype_100markers.vcf.gz    \
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     \
        --vcfField=GT   \
        --SAIGEOutputFile=./output/genotype_100markers_marker_vcf_step1withSparseGRM.txt \
        --chrom=1       \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary_sparseGRM.rda \
        --is_output_moreDetails=TRUE \
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   \
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt

對于二元性狀，可以通過Firth的偏差減少邏輯回歸更準(zhǔn)確地估計效應(yīng)大小，方法是設(shè)置

–is_Firth_beta=TRUE 和 –pCutoffforFirth=0.01

  • p 值為 <= pCutoffforFirth 的標(biāo)記的效應(yīng)大小將通過Firth 的偏倚減少邏輯回歸進行估計 
Rscript step2_SPAtests.R        \
        --vcfFile=./input/genotype_100markers.vcf.gz    \
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     \
        --vcfField=GT   \
        --SAIGEOutputFile=./output/genotype_100markers_marker_vcf_step1withSparseGRM.txt \
        --chrom=1       \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary_sparseGRM.rda \
        --is_output_moreDetails=TRUE    \
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   \
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt \
    --is_Firth_beta=TRUE    \
        --pCutoffforFirth=0.01

條件分析

  • -condition = 冒號分隔的變異id (chr:pos:ref:alt)
    比如：chr3:101651171:C:T,chr3:101651186:G:A
Rscript step2_SPAtests.R        \
        --vcfFile=./input/genotype_100markers.vcf.gz    \
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     \
        --vcfField=GT   \
        --SAIGEOutputFile=./output/genotype_100markers_marker_vcf_cond.txt \
        --chrom=1       \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary.rda \
        --varianceRatioFile=./output/example_binary.varianceRatio.txt   \
        --is_output_moreDetails=TRUE    \
        --condition=1:13:A:C,1:79:A:C

輸入文件

  1. （必需）劑量文件 SAIGE 支持不同的劑量格式： PLINK， VCF， BCF， BGEN 和 SAV.

    • plink
     ./input/genotype_100markers.bed
 ./input/genotype_100markers.bim
 ./input/genotype_100markers.fam
    • BGEN
     ./input/genotype_100markers.bgen
 ./input/genotype_100markers.bgen.bgi
    • 含有基因型的VCF
     #go to input
 cd ./input
 #index file can be generated using tabix
 tabix --csi -p vcf genotype_10markers.vcf.gz
 zcat genotype_10markers.vcf.gz | less -S
 genotype_10markers.vcf.gz.csi
    • 含劑量的vcf
     #index file can be generated using tabix
 tabix --csi -p vcf dosage_10markers.vcf.gz    
 zcat dosage_10markers.vcf.gz | less -S
 dosage_10markers.vcf.gz.csi
    • SAV
     dosage_10markers.sav
 dosage_10markers.sav.s1r

  2. （可選。僅適用于不包含樣品 ID 的 BGEN 文件） 示例文件包含一列用于與劑量文件中的示例順序相對應(yīng)的示例 ID。不包括標(biāo)頭。該選項最初用于不包含示例信息的 BGEN 文件。

     less -S ./input/samplelist.txt
1a1
1a2
1a3

  3. （必需。步驟 1） 步驟 1 中的模型文件輸出

     ./output/example_binary.rda

  4. （可選。步驟 1） 步驟 1 中的方差比文件中的輸出

     ./output/example_binary.varianceRatio.txt

  5. （可選。與步驟 1） 稀疏 GRM 文件中的輸入和稀疏 GRM 的示例 ID 文件相同。有關(guān)更多詳細信息，請參閱步驟 0。

     ./output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx
 ./output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt

 ##the sparse matrix can be read and viewed using the R library Matrix
 library(Matrix)
 sparseGRM = readMM("./output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx")

輸出文件

  1. 包含關(guān)聯(lián)分析結(jié)果的文件

     less -S ./output/genotype_100markers_marker_vcf.txt
  CHR     POS     MarkerID        Allele1 Allele2 AC_Allele2      AF_Allele2      MissingRate    
 BETA    SE      Tstat   var     p.value p.value.NA      Is.SPA  p.value.NA_c    Is.SPA.converge AF_case  AF_ctrl N_case  N_ctrl
1       4       rs4     A       C       20      0.01    0       -0.145553       0.774596        -0.242588       1.66667 0.850949        0.850949        false   0.0102041       0.00997783      98       902     0       2       0       18

  2. 輸出文件的索引。此文件包含關(guān)聯(lián)測試的進度。如果在步驟 2 中指定 –is_overwrite_output=FALSE，程序?qū)z查此索引并繼續(xù)未完成的分析，而不是從頭開始

     less -S ./output/genotype_100markers_marker_vcf.txt.index 
 This is the output index file for SAIGE package to record the end point in case users want to restart the analysis. Please do not modify this file.
  This is a Marker level analysis.
  nEachChunk = 10000
  Have completed the analysis of chunk 1
  Have completed the analyses of all chunks.

注意：

  • 關(guān)聯(lián)結(jié)果與等位基因2有關(guān)。
#check the header
head -n 1 output/genotype_100markers_marker_vcf_cond.txt

CHR: chromosome
POS: genome position 
SNPID: variant ID
Allele1: allele 1
Allele2: allele 2
AC_Allele2: allele count of allele 2
AF_Allele2: allele frequency of allele 2
MissingRate: missing rate (If the markers in the dosage/genotype input are imputed with is_imputed_data=TRUE, imputationInfo will be output instead of MissingRate)
imputationInfo: imputation info (output with is_imputed_data=TRUE). If not in dosage/genotype input file, will output 1
BETA: effect size of allele 2
SE: standard error of BETA
Tstat: score statistic of allele 2
var: estimated variance of score statistic
p.value: p value (with SPA applied for binary traits)
p.value.NA: p value when SPA is not applied (only for binary traits)
Is.SPA.converge: whether SPA is converged or not (only for binary traits)

#if --condition= is used for conditioning analysis, the conditional analysis results will be output  
BETA_c, SE_c, Tstat_c, var_c, p.value_c, p.value.NA_c 

AF_case: allele frequency of allele 2 in cases (only for binary traits)
AF_ctrl: allele frequency of allele 2 in controls (only for binary traits)
N_case: sample size in cases (only for binary traits)
N_ctrl: sample size in controls (only for binary traits)

if --is_output_moreDetails=TRUE
homN_Allele2_cases, hetN_Allele2_cases, homN_Allele2_ctrls, hetN_Allele2_ctrls will be output for sample sizes with different genotypes in cases and controls (only for binary traits)

示例 1

  • 二元性狀（–性狀類型=二元）
  • 使用完整的 GRM 擬合空模型，該 GRM 將使用 plink 文件中的基因型動態(tài)計算 （–plinkFile=）
  • 估計步驟 1 中的方差比，該比率將用作步驟 2 的輸入
  • 在步驟 1 中使用 2 個 CPU （–n 線程）
  • a9 和 a10 是分類協(xié)變量，將針對不同級別重寫 （–qCovarColList）
    Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly_22chr  \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_binary \
        --covarColList=x1,x2,a9,a10 \
        --qCovarColList=a9,a10  \
        --sampleIDColinphenoFile=IID \
        --traitType=binary        \
        --outputPrefix=./output/example_binary \
        --nThreads=2    \
        --IsOverwriteVarianceRatioFile=TRUE

    Rscript step2_SPAtests.R        \
        --vcfFile=./input/genotype_100markers.vcf.gz    \
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     \
        --vcfField=GT   \
        --SAIGEOutputFile=./output/genotype_100markers_marker_vcf.txt \
        --chrom=1       \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary.rda \
        --varianceRatioFile=./output/example_binary.varianceRatio.txt

示例 2

  • 二元性狀（–性狀類型=二元）
  • 使用稀疏 GRM 擬合空模型 （–useSparseGRMtoFitNULL=TRUE， –sparseGRMFile， –sparseGRMSampleIDFile）
  • 不要在步驟 1 中估計方差比。
  • 當(dāng)使用稀疏 GRM 擬合空模型且不會應(yīng)用 LOCO 時，僅使用一個 CPU
  • 在步驟 2 中將 PLINK 輸入用于基因型/劑量
  • p 值為 <= pCutoffforFirth 的標(biāo)記的效應(yīng)大小將通過 Firth 的偏倚減少邏輯回歸 –is_Firth_beta=TRUE 和 –pCutoffforFirth=0.01 來估計（請注意，此選項正在評估中）
    Rscript step1_fitNULLGLMM.R     \
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   \
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     \
        --useSparseGRMtoFitNULL=TRUE    \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_binary \
        --covarColList=x1,x2,a9,a10 \
        --qCovarColList=a9,a10  \
        --sampleIDColinphenoFile=IID \
        --traitType=binary        \
        --outputPrefix=./output/example_binary_sparseGRM

    Rscript step2_SPAtests.R        \
        --bedFile=./input/genotype_100markers.bed       \
        --bimFile=./input/genotype_100markers.bim       \
        --famFile=./input/genotype_100markers.fam       \
        --AlleleOrder=alt-first \
        --SAIGEOutputFile=./output/genotype_100markers_marker_plink_step1withSparseGRM_Firth.txt \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_binary_sparseGRM.rda \
        --is_output_moreDetails=TRUE    \
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   \
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt \
        --is_Firth_beta=TRUE    \
        --pCutoffforFirth=0.01

示例 3

  • 數(shù)量特征（–性狀類型=定量）和需要逆歸一化（–invNormalize=TRUE）
  • 使用稀疏 GRM 擬合空模型 （–useSparseGRMtoFitNULL=TRUE， –sparseGRMFile， –sparseGRMSampleIDFile）
  • 在步驟 1 中估計方差比 （–plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly_22chr）
  • 當(dāng)使用稀疏 GRM 擬合空模型且不會應(yīng)用 LOCO 時，僅使用一個 CPU
  • –需要為 VCF 輸入指定色度
    Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly_22chr  \
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   \
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     \
        --useSparseGRMtoFitNULL=TRUE    \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_quantitative \
        --covarColList=x1,x2,a9,a10 \
        --qCovarColList=a9,a10  \
        --sampleIDColinphenoFile=IID \
        --traitType=quantitative        \
    --invNormalize=TRUE \
        --outputPrefix=./output/example_quantitative_sparseGRM

    Rscript step2_SPAtests.R        \
        --vcfFile=./input/genotype_100markers.vcf.gz    \
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     \
        --vcfField=GT   \
    --chrom=1   \
        --SAIGEOutputFile=./output/genotype_100markers_marker_vcf_step1withSparseGRM.txt \
        --minMAF=0 \
        --minMAC=20 \
        --GMMATmodelFile=./output/example_quantitative_sparseGRM.rda \
        --is_output_moreDetails=TRUE    \
        --varianceRatioFile=./output/example_quantitative_sparseGRM.varianceRatio.txt

示例 4

  • 二元性狀（–性狀類型=二元）
  • 使用完整的 GRM 擬合空模型，該 GRM 將使用 plink 文件中的基因型動態(tài)計算 （–plinkFile=）
  • 估計步驟 1 中的方差比，該比率將用作步驟 2 的輸入
  • 覆蓋現(xiàn)有的步驟 1 輸出 – Is覆蓋方差RatioFile=TRUE
  • 步驟 2 中標(biāo)記的 p.值為 ~3.74 x 10^-7
Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly \
        --phenoFile=./input/Prev_0.1_nfam_1000.pheno_positive_pheno.txt \
        --phenoCol=y \
        --covarColList=x1,x2 \
        --sampleIDColinphenoFile=IID \
        --traitType=binary        \
        --outputPrefix=./output/example_binary_positive_signal \
        --nThreads=4    \
    --IsOverwriteVarianceRatioFile=TRUE

Rscript step2_SPAtests.R        \
        --vcfFile=./input/nfam_1000_MAF0.2_nMarker1_nseed200.vcf.gz \
        --vcfFileIndex=./input/nfam_1000_MAF0.2_nMarker1_nseed200.vcf.gz.tbi \
        --vcfField=GT \
        --chrom=1 \
        --minMAF=0.0001 \
        --minMAC=1 \
        --GMMATmodelFile=./output/example_binary_positive_signal.rda \
        --varianceRatioFile=./output/example_binary_positive_signal.varianceRatio.txt \
        --SAIGEOutputFile=./output/example_binary_positive_signal.assoc.step2.txt

參考

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容