STAARpipeline

README.md

@@ -1,7 +1,7 @@
 # STAARpipeline-Tutorial
 This is a tutorial for (1) automatically functionally annotating the variants of whole-genome/whole-exome sequencing (WGS/WES) studies and integrating the functional annotations with the genotype data using **FAVORannotator** and (2) performing single-/multi-trait association analysis of WGS/WES studies, summarizing and visualization results using **STAARpipeline** and **STAARpipelineSummary**. The software prerequisites, dependencies and installation can be found in <a href="https://github.com/xihaoli/STAARpipeline">**STAARpipeline**</a> and <a href="https://github.com/xihaoli/STAARpipelineSummary">**STAARpipelineSummary**</a> packages.
 
-**FAVORannotator**, **STAARpipeline** and **STAARpipelineSummary** are implemented as a collection of apps. Please see the apps
+**FAVORannotator**, **STAARpipeline** and **STAARpipelineSummary** are implemented as a collection of apps. Please see the following apps
 
 **favorannotator** (<a href="https://platform.sb.biodatacatalyst.nhlbi.nih.gov/public/apps/admin/sbg-public-data/favorannotator-1-0-0">**BDC-Seven Bridges**</a>, <a href="https://github.com/xihaoli/favorannotator-rap">**RAP-DNAnexus**</a>) <br>
 **staarpipeline** (<a href="https://platform.sb.biodatacatalyst.nhlbi.nih.gov/public/apps/admin/sbg-public-data/staarpipeline-0-9-6">**BDC-Seven Bridges**</a>, <a href="https://github.com/xihaoli/staarpipeline-rap">**RAP-DNAnexus**</a>) <br>
@@ -65,13 +65,13 @@ R package **FastSparseGRM** provides functions and a pipeline to efficiently cal
 * `jobs_num.Rdata`: a data frame containing the number of jobs for association analysis, including individual analysis, sliding window analysis and dynamic window analysis (SCANG-STAAR).
 
 ### Step 1: Fit STAAR or MultiSTAAR null model
-#### Script: <a href="STAARpipeline_Null_Model.r">**STAARpipeline_Null_Model.r**</a> or <a href="STAARpipeline_Null_Model_GENESIS.r">**STAARpipeline_Null_Model_GENESIS.r**</a> or <a href="STAARpipeline_Null_Model_Multi.r">**STAARpipeline_Null_Model.r**</a>
+#### Script: <a href="STAARpipeline_Null_Model.r">**STAARpipeline_Null_Model.r**</a> or <a href="STAARpipeline_Null_Model_GENESIS.r">**STAARpipeline_Null_Model_GENESIS.r**</a> or <a href="STAARpipeline_Null_Model_Multi.r">**STAARpipeline_Null_Model_Multi.r**</a>
 * `STAARpipeline_Null_Model.r` fits the STAAR null model using the STAARpipeline package. <br>
 * `STAARpipeline_Null_Model_GENESIS.r` fits the null model using the GENESIS package and convert it to the STAAR null model using the STAARpipeline package. <br>
 * `STAARpipeline_Null_Model_Multi.r` fits the MultiSTAAR null model using the STAARpipeline package. <br>
 #### Input: Phenotype data and (sparse) genetic relatedness matrix. For more details, please see the R scripts.
 #### Output: a Rdata file of the STAAR or MultiSTAAR null model.
-Note: Once the STAAR or MultiSTAAR null model is fit, all the remaining steps of STAARpipeline and STAARpipelineSummary share the same scripts (the information of single-trait or multi-trait analysis being considered is automatically retrieved from the null model object).
+Note: Once the <a href="https://github.com/xihaoli/STAAR">STAAR</a> or <a href="https://github.com/xihaoli/MultiSTAAR">MultiSTAAR</a> null model is fit, all the remaining steps of STAARpipeline and STAARpipelineSummary share the same scripts (the information of single-trait or multi-trait analysis being considered is automatically retrieved from the null model object).
 
 ### Step 2: Individual (single-variant) analysis
 #### Script: <a href="STAARpipeline_Individual_Analysis.r">**STAARpipeline_Individual_Analysis.r**</a>

STAARpipeline_Null_Model.r

@@ -29,9 +29,9 @@ output_name <- "obj_nullmodel.Rdata"
 #           Main Function 
 ###########################################################
 ## fit null model
-obj_nullmodel <- fit_nullmodel(LDLadj.norm~age+age2+sex+PC1+PC2+PC3+PC4+PC5+PC6+PC7+PC8+PC9+PC10+as.factor(study_race),
+obj_nullmodel <- fit_nullmodel(LDLadj.norm~age+age2+sex+PC1+PC2+PC3+PC4+PC5+PC6+PC7+PC8+PC9+PC10+as.factor(study_ethnicity),
                                data=phenotype,kins=sgrm,use_sparse=TRUE,kins_cutoff=0.022,id="sample.id",
-                               groups="study_race",family=gaussian(link="identity"),verbose=TRUE)
+                               groups="study_ethnicity",family=gaussian(link="identity"),verbose=TRUE)
 
 save(obj_nullmodel,file=paste0(output_path,output_name))
 

STAARpipeline_Null_Model_GENESIS.r

@@ -29,8 +29,8 @@ output_name <- "obj_nullmodel_GENESIS.Rdata"
 ## fit null model using GENESIS 
 data_GENESIS <- as(phenotype,"AnnotatedDataFrame") # Make AnnotatedDataFrame (specifically required by GENESIS)
 obj_nullmodel_GENESIS <- fitNullModel(data_GENESIS,outcome="LDLadj.norm",
-                                      covars=c("age","age2","sex","PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","study_race"),
-                                      cov.mat=sgrm,group.var="study_race",AIREML.tol=1e-4,verbose=TRUE)
+                                      covars=c("age","age2","sex","PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","study_ethnicity"),
+                                      cov.mat=sgrm,group.var="study_ethnicity",AIREML.tol=1e-4,verbose=TRUE)
 
 ## convert GENESIS null model to STAAR null model
 obj_nullmodel <- genesis2staar_nullmodel(obj_nullmodel_GENESIS)

STAARpipeline_Null_Model_Multi.r

@@ -0,0 +1,38 @@
+###########################################################
+# fit STAAR null model
+# Xihao Li, Zilin Li
+# Initiate date: 11/04/2021
+# Current date: 06/13/2024
+###########################################################
+rm(list=ls())
+gc()
+
+library(gdsfmt)
+library(SeqArray)
+library(SeqVarTools)
+library(STAAR)
+library(MultiSTAAR)
+library(STAARpipeline)
+
+###########################################################
+#           User Input
+###########################################################
+## Phenotype file
+phenotype <- read.csv("/path_to_the_file/pheno.csv")
+## (sparse) GRM file
+sgrm <- get(load("/path_to_the_file/sGRM.Rdata"))
+## file directory for the output file 
+output_path <- "/path_to_the_output_file/"
+## output file name
+output_name <- "obj_nullmodel.Rdata"
+
+###########################################################
+#           Main Function 
+###########################################################
+## fit null model
+obj_nullmodel <- fit_nullmodel(cbind(LDLadj.norm,HDLadj.norm,TGadj.norm)~age+age2+sex+PC1+PC2+PC3+PC4+PC5+PC6+PC7+PC8+PC9+PC10+PC11+as.factor(study_ethnicity),
+                               data=phenotype,kins=sgrm,use_sparse=TRUE,kins_cutoff=0.022,id="sample.id",
+                               family=gaussian(link="identity"),verbose=TRUE)
+
+save(obj_nullmodel,file=paste0(output_path,output_name))
+