SAIGE用户手册笔记2

基于集合的分析

SAIGE-GENE(现在称为SAIGE-GENE+)采取两个步骤来执行基于集合的关联测试

  • BURDEN, SKAT, and SKAT-O
  • MAC < = 10的超罕见变体在测试中被归结为伪marker
  • 多个函数注释,例如仅 Lof、Lof+Missense、Lof+Missense+Synonymous
  • 多个最大 MAF 阀值,例如 0.0001、0.001 和 0.01
  • 对于每个基因,基于柯西分布组合对应于不同注释的多个p值和最大MAF截止值

步骤 1:拟合零逻辑/线性混合模型

  1. 当使用稀疏 GRM 拟合空模型(-useSparseGRMtoFitNULL=TRUE)且未建立方差比(–skipVarianceRatioEstimation=TRUE)时,基于集合的检验的步骤 1 与 SAIGE 中的单变量检验相同

    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        
        --skipVarianceRatioEstimation=TRUE 
        --outputPrefix=./output/example_binary_sparseGRM
  2. 当使用完整的 GRM 来拟合空模型时(GRM 是使用 PLINK 文件中的基因型动态构建的,–plinkfile=)

     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     
         --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_fullGRM 
         --nThreads=64    
         --useSparseGRMtoFitNULL=FALSE    
         --isCateVarianceRatio=TRUE      
         --outputPrefix_varRatio=./output/example_binary_cate 
         --useSparseGRMforVarRatio=TRUE  
         --IsOverwriteVarianceRatioFile=TRUE
    • 使用完整的 GRM 来拟合空模型(默认情况下,–useSparseGRMtoFitNULL=FALSE)
    • 使用完整的 GRM(使用 PLINK 文件中的基因型进行构造)和带有 –useSparseGRMforVarRatio=TRUE 的稀疏 GRM 来估计方差比 ** 对包含稀疏 GRM 的文件使用 –稀疏GRMFile 对包含稀疏 GRM 中样本 ID 的文件使用 –稀疏GRMSampleIDFile
    • 需要根据不同的次要等位基因计数类别估计多个方差比,其中 –isCateVarianceRatio=TRUE ** 默认情况下,10 <= MAC < 20 和 MAC >= 20 时,估计两个方差比。** 使用 –cateVarRatioMinMACVecExclude 和 –cateVarRatioMaxMACVecInclude 修改 MAC 类别 请注意,PLINK 文件需要包含至少 200 个 MAC 属于这些类别的变体 与 SAIGE 中用于 SAIGE 中单变量测试的步骤 1 不同,SAIGE 中仅估计单个方差比

输入文件(与 SAIGE 步骤 1 中的输入相同)

  • 与 SAIGE 在步骤 1 中用于单变量关联测试的输入相同

  • 对于稀有变体,PLINK 文件 –plinkFile= 特定于 SAIGE-GENE+,需要包含 MAC 低至 10 的稀有变体。

  1. (必填)表型文件(包含协变量(如果有),如性别和年龄)文件可以是空格,也可以是用标题以制表符分隔的。该文件必须包含一列用于样本 ID,一列用于表型。它可能包含协变量列。
image
less -S ./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt
  1. (可选)用于构建完整遗传关系矩阵(GRM)和估计方差比的基因型 文件 SAIGE 采用基因型的 PLINK 二进制文件,并假设文件前缀与 .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 文件。有关更多详细信息,请参阅步骤 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")

输出文件((与 SAIGE 步骤 1 中用于单变量关联测试的输出相同))

  • 与 SAIGE 在步骤 1 中用于单变量关联测试的输出相同

  • 特定于稀有变异的 SAIGE-GENE+,方差比文件包含多个方差比,而不是单个方差比

  1. 模型文件(步骤 2 的输入))

     ./output/example_binary.rda

     #load the model file in R
     load("./output/example_binary.rda")
     names(modglmm)
     modglmm$theta
    image
  2. 方差比文件(如果在步骤 1 中估计了类别方差比,则将在步骤 2 中输入)

     less -S ./output/example_binary_cate.varianceRatio.txt

  3. 随机选择的标记子集的关联结果文件(如果在步骤 1 中估计了方差比,则将生成此文件。它是一个中间文件,后续步骤不需要它。)

     less -S ./output/example_binary_30markers.SAIGE.results.txt

第2步:执行基于区域或基因的关联测试

  • 这些命令与单变体联合测试的步骤 2 相同,不同之处在于
    • 指定了一个组文件 (–groupFile),其中包含要测试的每个集合的遗传标记 ID、注释和权重(如果有)
  • 允许每组(基因或区域)使用多个掩码
    • 例如,singletons和doubletons为 1,2
    • 默认情况下,不应用此选项
    • 最大 MAC 将转换为最大 MAF,并在分析中合并到 –maxMAF_in_groupTest
    • 例如 0.0001,0.001,0.01(默认值)
    • 例如,”lof,missense:lof,missense:lof:synonymous”仅用于测试lof,missense +lof和missense+lof+synonymous
    • 使用 –annotation_in_groupTest 列出用逗号分隔的不同批注。在每个批注组合中,批注由”:”分隔
    • 使用 –maxMAF_in_groupTest 表示以逗号分隔的不同最大MAF 截止值
    • 使用 –maxMAC_in_groupTest 表示用逗号分隔的不同最大MAC 截止值
    • 在示例中,每个集合将应用 9 个掩码,并且将基于柯西组合合并9个p值
  • 默认情况下,将执行 SKAT-O 测试(同时输出 BURDEN 和 SKAT 测试结果)。使用 –r.corr=1 仅执行 BURDEN 检验
    • 如果执行 SKAT-O 测试(–r.corr=0),则还会输出单变量联合测试结果
    • 如果仅执行 BURDEN 检验(–r.corr=1),则默认情况下不执行单变量联合分析。使用 –is_single_in_groupTest=TRUE 输出单变量 关联分析结果。
  • 使用 –is_output_markerList_in_groupTest=TRUE 输出每个测试中包含的标记列表。
  • 默认情况下,程序首先检查组文件中是否提供了每个标记权重。如果不是,程序将基于 Beta 分布中的 MAF 与 paraemters weights.beta 计算权重。使用 –is_no_weight_in_groupTest=TRUE 可不在分析中使用任何权重。
  • 使用 –minMAF、–minMAC 和 –minInfo 指定的单个标记的截止值也应用于基于区域/集的分析中
  • 与单变量关联检验相同,可以执行基于条件分析的汇总统计信息(–条件)
  1. 在步骤 1 中,如果使用稀疏 GRM 来拟合空模型,并且没有估计方差比,则在步骤 2 中,使用与输入相同的稀疏 GRM(–稀疏 GRMFile、–稀疏 GRMSampleIDFile) 作为输入

     Rscript step2_SPAtests.R        
         --bgenFile=./input/genotype_100markers.bgen    
         --bgenFileIndex=./input/genotype_100markers.bgen.bgi 
         --SAIGEOutputFile=./output/genotype_100markers_bgen_groupTest_out_sparseGRMforStep1.txt 
         --chrom=1 
         --AlleleOrder=ref-first 
         --minMAF=0 
         --minMAC=0.5 
         --sampleFile=./input/samplelist.txt 
         --GMMATmodelFile=./output/example_binary_sparseGRM.rda 
         --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
         --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     
         --groupFile=./input/group_new_chrposa1a2.txt    
         --annotation_in_groupTest=lof,missense:lof,missense:lof:synonymous        
         --maxMAF_in_groupTest=0.0001,0.001,0.01
  2. 如果在步骤 1 中使用了完整的 GRM 来拟合空模型,并且使用完整和稀疏的 GRM 估计方差比,则在步骤 2 中,稀疏 GRM(–稀疏 GRMFile、–稀疏GRMSampleIDFile)和方差比 (–varianceRatioFile) 将用作输入。使用 –is_output_markerList_in_groupTest=TRUE 输出用于每个测试的标记。

     Rscript step2_SPAtests.R        
         --bgenFile=./input/genotype_100markers.bgen    
         --bgenFileIndex=./input/genotype_100markers.bgen.bgi 
         --SAIGEOutputFile=./output/genotype_100markers_bgen_groupTest_out.txt 
         --chrom=1 
         --LOCO=TRUE    
         --AlleleOrder=ref-first 
         --minMAF=0 
         --minMAC=0.5 
         --sampleFile=./input/samplelist.txt 
         --GMMATmodelFile=./output/example_binary_fullGRM.rda 
         --varianceRatioFile=./output/example_binary_cate.varianceRatio.txt 
         --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
         --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt  
         --groupFile=./input/group_new_chrposa1a2.txt 
         --annotation_in_groupTest="lof,missense:lof,missense:lof:synonymous"        
         --maxMAF_in_groupTest=0.0001,0.001,0.01 
         --is_output_markerList_in_groupTest=TRUE
  3. 仅执行 –r.corr=1 的 BURDEN 检验。使用 –minGroupMAC_in_BurdenTest表示负担测试中测试”负载标记”的最小 MAC。

     Rscript step2_SPAtests.R        
         --bgenFile=./input/genotype_100markers.bgen    
         --bgenFileIndex=./input/genotype_100markers.bgen.bgi 
         --SAIGEOutputFile=./output/genotype_100markers_bgen_groupTest_out_onlyBURDEN.txt 
         --chrom=1 
         --LOCO=TRUE    
         --AlleleOrder=ref-first 
         --minMAF=0 
         --minMAC=0.5 
         --sampleFile=./input/samplelist.txt 
         --GMMATmodelFile=./output/example_binary_fullGRM.rda 
         --varianceRatioFile=./output/example_binary_cate.varianceRatio.txt      
         --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
         --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt 
         --groupFile=./input/group_new_chrposa1a2.txt    
         --annotation_in_groupTest="lof,missense;lof,missense;lof;synonymous"        
         --maxMAF_in_groupTest=0.0001,0.001,0.01 
     --r.corr=1
  4. Use –condition= to perform conditioning analysis

     Rscript step2_SPAtests.R        
         --bgenFile=./input/genotype_100markers.bgen    
         --bgenFileIndex=./input/genotype_100markers.bgen.bgi 
         --SAIGEOutputFile=./output/genotype_100markers_bgen_groupTest_out_cond.txt 
         --chrom=1 
         --LOCO=TRUE    
         --AlleleOrder=ref-first 
         --minMAF=0 
         --minMAC=0.5 
         --sampleFile=./input/samplelist.txt 
         --GMMATmodelFile=./output/example_binary_fullGRM.rda 
         --varianceRatioFile=./output/example_binary_cate.varianceRatio.txt      
         --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
         --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     
         --groupFile=./input/group_new_chrposa1a2.txt    
         --annotation_in_groupTest=lof,missense:lof,missense:lof:synonymous        
         --maxMAF_in_groupTest=0.0001,0.001,0.01 
         --condition=1:30:A:C,1:79:A:C
  5. 使用PLINK文件作为基因型/剂量的输入(–bedFile=, –bimFile=, –famFile=, –AlleleOrder=)

     Rscript step2_SPAtests.R        
         --bedFile=./input/genotype_100markers.bed       
         --bimFile=./input/genotype_100markers.bim       
         --famFile=./input/genotype_100markers.fam       
         --SAIGEOutputFile=./output/genotype_100markers_plink_groupTest_out.txt 
         --chrom=1 
         --LOCO=TRUE    
         --AlleleOrder=alt-first 
         --minMAF=0 
         --minMAC=0.5 
         --sampleFile=./input/samplelist.txt 
         --GMMATmodelFile=./output/example_binary_fullGRM.rda 
         --varianceRatioFile=./output/example_binary_cate.varianceRatio.txt      
         --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
         --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     
         --groupFile=./input/group_new_chrposa1a2.txt    
         --annotation_in_groupTest=lof,missense:lof,missense:lof:synonymous        
         --maxMAF_in_groupTest=0.0001,0.001,0.01
    • –等位基因顺序可以是 alt-first 或 ref-first。必须正确指定,否则,PLINK 文件中的变体将与 groupFile 中的标记不匹配
  6. VCF 文件作为输入

    Rscript step2_SPAtests.R        
        --vcfFile=./input/genotype_100markers.vcf.gz    
        --vcfFileIndex=./input/genotype_100markers.vcf.gz.csi     
        --vcfField=GT   
        --SAIGEOutputFile=./output/genotype_100markers_vcf_groupTest_out.txt 
        --LOCO=FALSE    
        --chrom=1 
        --minMAF=0 
        --minMAC=0.5 
        --sampleFile=./input/samplelist.txt 
        --GMMATmodelFile=./output/example_binary_fullGRM.rda 
        --varianceRatioFile=./output/example_binary_cate.varianceRatio.txt      
        --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
        --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     
        --groupFile=./input/group_new_chrposa1a2.txt    
        --annotation_in_groupTest=lof,missense:lof,missense:lof:synonymous        
        --maxMAF_in_groupTest=0.0001,0.001,0.01
    • –vcfFileIndex 将 .csi 索引文件作为输入,可以使用 tabix –csi -p vcf ./input/genotype_100markers.vcf.gz
    • –vcfField 是 GT 或 DS
    • –必须为 VCF 文件指定 chrom。–chrom必须与 VCF 文件中的CHROM字符串相同

输入文件

详情请参考 SAIGE (Single-variant test) Step 2 输入文件。此外,基于集合的测试还需要额外的输入文件

  1. (必填。特定于基于集合的测试)组文件,其中包含每组(基因或区域)的标记ID,注释和/或权重。* 第一列包含设置的名称。* 组文件每组有 2 或 3 行。* 标记 ID 和注释是必需的,权重是可选的。* 第二列具有 var(指示该行用于标记 ID)、anno(指示该行用于注释)和权重(指示该行用于基于集的测试中使用的标记的权重)

    less -S ./input/group_new_chrposa1a2.txt

    less -S ./input/group_new_chrposa1a2_withWeights.txt

    • 带权重的分组文件
    • 如果权重未包含在组文件中,则默认情况下权重将计算为 beta(MAF,1,25)

    • 不带权重的组文件

输出文件

  • 包含基于区域或基因的关联测试结果的文件

      head -n 1  ./output/genotype_100markers_bgen_groupTest_out_cond.txt

    image
Region: set name
Group: annotation mask
max_MAF: maximum MAF cutoff
Pvalue: p value for SKAT-O test
Pvalue_Burden: p value for BURDEN test
Pvalue_SKAT: p value for SKAT test
BETA_Burden: effect size of BURDEN test
SE_Burden: standard error of BETA_Burden

#if --condition= is used for conditioning analysis, the conditional analysis results will be output
Pvalue_cond, Pvalue_Burden_cond, Pvalue_SKAT_cond, BETA_Burden_cond, SE_Burden_cond

MAC: minor allele count in the set
MAC_case: minor allele count in cases
MAC_control: minor allele count in controls
Number_rare: number of markers that are not ultra-rare with MAC > MACCutoff_to_CollapseUltraRare (=10 by default)
Number_ultra_rare: number of markers that are ultra-rare with MAC <= MACCutoff_to_CollapseUltraRare (=10 by default)

  • 包含基于集合的测试中单个标记的关联测试结果的文件 有关详细信息,请参阅 SAIGE(单变量测试)步骤 2 输出文件。 对于二元性状,通过设置 * –is_Firth_beta=TRUE 和 –pCutoffforFirth=0.01 * 注意:此选项目前仅适用于单变量协会检验和负担检验(当执行唯一的负荷检验 (–r.corr=1) 时),可以通过 Firth 的偏倚减少逻辑回归更准确地估计单个变异的效应大小。当执行 SKAT-O 检验 (–r.corr=0) 时,不应用 Firth 的偏倚减少逻辑回归。** 如果仅执行 Burden 检验 (–r.corr=1),请使用 –is_single_in_groupTest=TRUE 输出单变量 assoc 检验 ** 如果执行了 SKAT-O 检验 (–r.corr=0),则会自动输出单变异体联合检验结果。

      less -S ./output/genotype_100markers_bgen_groupTest_out_cond.txt.singleAssoc.txt

  • 包含基于区域/集的测试的标记列表的文件 (–is_output_markerList_in_groupTest=TRUE)

      less -S ./output/genotype_100markers_bgen_groupTest_out.txt.markerList.txt

示例 1

  • 使用稀疏 GRM 拟合空模型 (–useSparseGRMtoFitNULL=TRUE)并且没有方差比被估计(–skipVarianceRatioEstimation=TRUE)

  • a9 和 a10 是分类协变量,将针对不同级别重写 (–qCovarColList)

  • 需要在步骤 1 和步骤 2 中使用相同的稀疏 GRM 文件(对于方差比方法)

  • 为每组(基因或区域)使用多个掩码

      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        
          --skipVarianceRatioEstimation=TRUE      
          --outputPrefix=./output/example_binary_sparseGRM

      Rscript step2_SPAtests.R        
          --bgenFile=./input/genotype_100markers.bgen    
          --bgenFileIndex=./input/genotype_100markers.bgen.bgi 
          --SAIGEOutputFile=./output/genotype_100markers_bgen_groupTest_out_sparseGRMforStep1.txt 
          --chrom=1 
          --AlleleOrder=ref-first 
          --minMAF=0 
          --minMAC=0.5 
          --sampleFile=./input/samplelist.txt 
          --GMMATmodelFile=./output/example_binary_sparseGRM.rda 
          --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
          --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     
          --groupFile=./input/group_new_chrposa1a2.txt    
          --annotation_in_groupTest=lof,missense:lof,missense:lof:synonymous        
          --maxMAF_in_groupTest=0.0001,0.001,0.01
    • 应用了三个批注masks:仅 lof、missense+lof 和 missense+lof+同义词 (–annotation_in_groupTest)
    • 应用三个最大 MAF 截止值:0.0001,0.001,0.01
    • 每个集合的 3×3 检验的 p 值基于柯西组合进行组合

示例 2

  • 使用完整的 GRM 拟合空模型,该 GRM 将使用 plink 文件中的基因型动态计算 (*–plinkFile=*, –useSparseGRMtoFitNULL=FALSE)

  • a9 和 a10 是分类协变量,针对不同级别重编码(–qCovarColList)

  • 使用完整的 GRM(使用 PLINK 文件中的基因型进行构造)和sparseGRM 与 –useSparseGRMforVarRatio=TRUE *对包含稀疏 GRM 的文件使用 –sparseGRMFile 对包含稀疏 GRM 中样本 ID 的文件使用 –sparseGRMSampleIDFile

  • 需要根据不同的次要等位基因计数类别估计多个方差比,其中 –isCateVarianceRatio=TRUE ** 默认情况下,10 <= MAC < 20 和 MAC >= 20 时,估计两个方差比。使用 –cateVarRatioMinMACVecExclude 和 –cateVarRatioMaxMACVecInclude 修改 MAC 类别 请注意,PLINK 文件需要包含至少 1000 个 MAC 属于这些类别的变体 ** 与 SAIGE 中用于 SAIGE 中单变量测试的步骤 1 不同,SAIGE 中仅估计单个方差比

  • 需要在步骤 1 和步骤 2 中使用相同的稀疏 GRM 文件(对于方差比方法)

  • 输出测试的标记列表 – is_output_markerList_in_groupTest=TRUE

      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     
          --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_fullGRM 
          --nThreads=64    
          --useSparseGRMtoFitNULL=FALSE    
          --isCateVarianceRatio=TRUE      
          --outputPrefix_varRatio=./output/example_binary_cate 
          --useSparseGRMforVarRatio=TRUE  
          --IsOverwriteVarianceRatioFile=TRUE

      Rscript step2_SPAtests.R        
          --bgenFile=./input/genotype_100markers.bgen    
          --bgenFileIndex=./input/genotype_100markers.bgen.bgi 
          --SAIGEOutputFile=./output/genotype_100markers_bgen_groupTest_out.txt 
          --chrom=1 
          --LOCO=TRUE    
          --AlleleOrder=ref-first 
          --minMAF=0 
          --minMAC=0.5 
          --sampleFile=./input/samplelist.txt 
          --GMMATmodelFile=./output/example_binary_fullGRM.rda 
          --varianceRatioFile=./output/example_binary_cate.varianceRatio.txt      
          --sparseGRMFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx   
          --sparseGRMSampleIDFile=output/sparseGRM_relatednessCutoff_0.125_1000_randomMarkersUsed.sparseGRM.mtx.sampleIDs.txt     
          --groupFile=./input/group_new_chrposa1a2.txt    
          --annotation_in_groupTest=lof,missense:lof,missense:lof:synonymous        
          --maxMAF_in_groupTest=0.0001,0.001,0.01 
      --is_output_markerList_in_groupTest=TRUE

创建稀疏GRM

步骤 0:创建稀疏 GRM

SAIGE 和 SAIGE-GENE 可以采用稀疏 GRM 来拟合零模型和关联分析

  • 这种稀疏的GRM只需要为每个数据集创建一次,例如生物库,并且只要所有测试的样品都在稀疏GRM中,就可以用于所有不同的表型。
  • 可以使用多个程序来生成稀疏 GRM
  1. SAIGE 提供了一个用于创建稀疏 GRM 的脚本 *程序将输出一个以 sampleID 结尾的文件.txt其中包含稀疏 GRM 的示例 ID,以及一个以 .sparseGRM.mtx 结尾的文件,其中包含稀疏 GRM

     #For help information
     Rscript createSparseGRM.R --help

     Rscript createSparseGRM.R       
         --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly 
         --nThreads=4  
         --outputPrefix=./output/sparseGRM       
         --numRandomMarkerforSparseKin=2000      
         --relatednessCutoff=0.125
    • 然后,这两个文件可以直接用于后续步骤
  2. GCTA[1]

     gcta64 
         --bfile ./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly 
         --out ./output/sparseGRM 
         --make-grm-part 3 1 
         --maf 0.01 
         --geno 0.15 
         --thread-num 2
  3. KING[2]

同时运行单变异和集合的关联分析

如果 SAIGE 以前曾用于单变异关联测试。步骤 1 中的空模型拟合结果可重复用于基于区域/集的分析 首先,使用完整和稀疏的 GRM 估计分类方差比 使用 –sparseGRMFile 和 –sparseGRMSampleIDFile 作为稀疏 GRM 和 –useSparseGRMforVarRatio=TRUE 使用 –isCateVarianceRatio=TRUE 来构造分类方差比 使用 –skipModelFitting=FASLE 以避免重新拟合空模型。使用 –outputPrefix=prefix 时,prgoram 将从 prefix.rda 读取模型。如果尝试避免覆盖以前的方差比文件,请使用 –outputPrefix_varRatio为新的方差比结果指定单独的文件前缀,否则 –IsOverwriteVarianceRatioFile=TRUE 可用于覆盖以前的文件。注意:在使用 –plinkFile 指定的 plink 文件中,至少需要添加约 200 个 10<= MAC < 20 的标记,这将用于估计较低 MAF 类别的方差比。

  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     
      --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_fullGRM 
      --nThreads=64    
      --useSparseGRMtoFitNULL=FALSE    
      --isCateVarianceRatio=TRUE      
  --skipModelFitting=TRUE 
      --outputPrefix_varRatio=./output/example_binary_cate 
      --useSparseGRMforVarRatio=TRUE

如果之前未运行任何 SAIGE 作业,请按照步骤 1 和 2 进行基于集的测试。使用 –is_single_in_groupTest=TRUE 同时输出单变量联合测试

参考资料

[1]

GCTA: https://yanglab.westlake.edu.cn/software/gcta/#MakingaGRM

[2]

KING: https://www.kingrelatedness.com/manual.shtml





欢迎微信交流






小编微信






扫描二维码

获取更多精彩

公众号




本篇文章来源于微信公众号: 微因

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注