Purpose of the analysis

This workflow aims to use GWAS summary statistics data as input and implements the causal TWAS algorithm. The “causal TWAS” method is developed to identify causal genes in GWAS studies. We improve the current TWAS analysis mainly by reducing false positives. We start from z scores for each SNP in a GWAS study and expression model trained by the Fusion/TWAS software. In the end, we will provide a PIP (the posterior inclusion probality for being a causal gene) for each gene with an expression model.

Below is the workflow for running ctwas for summary statistics data. We provided example input files with the package and you can use the following code as a toy pipeline. You can prepare your own input files following the examples.

Prepare Input

Our package is called ctwas.

library(ctwas)

The example inputs are given as below:

  • z scores, a data frame:
# z scores
data(z_snp)
head(z_snp)
#>       id A1 A2          z
#> 1 rs1001  G  A  -6.644254
#> 2 rs1002  C  T  14.369485
#> 3 rs1003  G  A -14.180414
#> 4 rs1004  C  T  -6.739014
#> 5 rs1005  G  A -14.049735
#> 6 rs1006  C  T  -2.176475

Here, A1 is effect allele. A2 is the other allele. z indicates the z score from GWAS.

  • LD reference data.

There are two ways to provide LD information for the functions in ctwas. One is to provide LD genotype data (in plink format), the other way is to provide genetic correlation matrix (R matrix) for LD indpendent regions.

If you want to provide LD genotype type, then you need to provide a character vector of .pgen or .bed files. One file for one chromosome, in the order of 1 to 22. Therefore, the length of this vector needs to be 22. If .pgen files are given, then .pvar and .psam are assumed to present in the same directory. If .bed files are given, then .bim and .fam files are assumed to present in the same directory. For example:

# genotype for LD reference
ld_pgenfs <- system.file("extdata/example_genotype_files", paste0("example_chr", 1:22, ".pgen"), package = "ctwas")
ld_pgenfs
#>  [1] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr1.pgen" 
#>  [2] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr2.pgen" 
#>  [3] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr3.pgen" 
#>  [4] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr4.pgen" 
#>  [5] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr5.pgen" 
#>  [6] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr6.pgen" 
#>  [7] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr7.pgen" 
#>  [8] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr8.pgen" 
#>  [9] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr9.pgen" 
#> [10] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr10.pgen"
#> [11] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr11.pgen"
#> [12] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr12.pgen"
#> [13] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr13.pgen"
#> [14] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr14.pgen"
#> [15] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr15.pgen"
#> [16] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr16.pgen"
#> [17] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr17.pgen"
#> [18] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr18.pgen"
#> [19] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr19.pgen"
#> [20] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr20.pgen"
#> [21] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr21.pgen"
#> [22] "/tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_genotype_files/example_chr22.pgen"

Or if you want to provide R matrices , you need to provide a directory with all .RDS and matching .Rvar files. For example:

ld_R_dir <- system.file("extdata/example_ld_R", package = "ctwas")
list.files(ld_R_dir)[1:10]
#>  [1] "example_chr1.R_snp.0_5.RDS"      "example_chr1.R_snp.0_5.Rvar"    
#>  [3] "example_chr1.R_snp.125_155.RDS"  "example_chr1.R_snp.125_155.Rvar"
#>  [5] "example_chr1.R_snp.155_185.RDS"  "example_chr1.R_snp.155_185.Rvar"
#>  [7] "example_chr1.R_snp.185_215.RDS"  "example_chr1.R_snp.185_215.Rvar"
#>  [9] "example_chr1.R_snp.215_245.RDS"  "example_chr1.R_snp.215_245.Rvar"

The .RDS file, is a R .RDS format file. It stores the LD correlation matrix for a region (a \(p \times p\) matrix, \(p\) is the number of SNPs in the region). We also require that for each .RDS file, a file with the same stem name but ended with the suffix .Rvar to present in the same directory. This .Rvar files gives SNPs information in the region and its order should match columns/rows in the .RDS file. An example, if the .RDS file is example_chr1.R_snp.0_5.RDS. the matching .Rvar file should be example_chr1.R_snp.0_5.Rvar in the same directory, and it looks like this:

head(read.table(system.file("extdata/example_ld_R","example_chr1.R_snp.0_5.Rvar",  package = "ctwas"), header = T))
#>   chrom     id pos alt ref
#> 1     1 rs1001   1   G   A
#> 2     1 rs1002   2   C   T
#> 3     1 rs1003   3   G   A
#> 4     1 rs1004   4   C   T

The program will automatically search for all .RDS and .Rvar file in the directory, so no other .RDS files should be seen in this directory unless you want to use it as a LD R matrix in the program. Each SNP should uniquely belong to one of these R matrices.

Note: LD reference should contain as many GWAS SNPs as possible as only the overlapping SNPs are included in the analysis. Also only SNPs with z scores, in LD reference and are eQTLs (i.e. SNPs shared in GWAS, LD reference and weight files) are used in imputing gene expression z scores. Thus we suggest you to impute z scores for eQTLs if many are not available. The program will harmonize data internally.

  • Weight files. We accept two formats for weight information. 1. FUSION/TWAS format. Please check out Fusion/TWAS for the format of weights. Below is an example.
weight.fusion <- system.file("extdata/example_fusion_weights", "Tissue", package = "ctwas")

If you want to provide weight in FUSION format then you just to provide the directory that contains all the .rda files like above. We assume a file with same name as the directory but has the suffix .pos is present in the same level as the directory, the program will search for this file automatically. For example, we have both the directory Tissue/ and Tissue.pos present under the extdata/example_fusion_weights folder.

  1. Predictdb format. Please see here for details: http://predictdb.org/. An example is provided with the package:
weight.predictdb <- system.file("extdata", "example_tissue.db", package = "ctwas") 

If you want to provide weight in predict.db format, just download these files from http://predictdb.org/ and provide the desired .db file.

  • Regions file (Optional).

If you want to use your own regions file, please provide a file like this:

regionsfile <- system.file("extdata", "example_regions.bed", package = "ctwas")
head(read.table(regionsfile, header = T))
#>   chr start stop
#> 1   1     0    5
#> 2   1     5   35
#> 3   1    35   65
#> 4   1    65   95
#> 5   1    95  125
#> 6   1   125  155

The regions are left closed and right open, i.g. [start, stop).

Otherwise, the package provided regions files generated by ldetect. You can specify the population and genome build that matches your GWAS data in ctwas_rss function, using the ld_regions and ld_regions_version arguments respectively. See below.

Steps when using LD genotype data

Step 1: impute expression z scores

In the following code, we compute expression z scores for each gene. We do this chromosome by chromosome. impute_expr_z will return z scores for each gene. It will also generate files containing the reference LD genotypes for eQTLs in the gene and the filename is returned. We will need these files in the ctwas_rss function in step 2.

outputdir <- "~/temp"
res <- impute_expr_z(z_snp = z_snp,  weight = weight.predictdb, ld_pgenfs = ld_pgenfs,
                           method = "lasso", outputdir = outputdir,
                           outname = "test_ss")
#> 2021-06-06 20:20:47 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:47 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:47 INFO::Reading weights for chromosome 1
#> 2021-06-06 20:20:47 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:47 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:47 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:47 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:47 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:47 INFO::Imputation done: number of genes with imputed expression: 1 for chr 1
#> 2021-06-06 20:20:47 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:47 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:47 INFO::Reading weights for chromosome 2
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 2
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 3
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 3
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 4
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 4
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 5
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 5
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 6
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 6
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 7
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 7
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 8
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 8
#> 2021-06-06 20:20:48 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:48 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:48 INFO::Reading weights for chromosome 9
#> 2021-06-06 20:20:48 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:48 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:48 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:48 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:48 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:48 INFO::Imputation done: number of genes with imputed expression: 1 for chr 9
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 10
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 10
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 11
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 11
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 12
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 12
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 13
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 13
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 14
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 14
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 15
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 15
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 16
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 16
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 17
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:49 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:49 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:49 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:49 INFO::Imputation done: number of genes with imputed expression: 1 for chr 17
#> 2021-06-06 20:20:49 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:49 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:49 INFO::Reading weights for chromosome 18
#> 2021-06-06 20:20:49 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:49 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:50 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:50 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:50 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:50 INFO::Imputation done: number of genes with imputed expression: 1 for chr 18
#> 2021-06-06 20:20:50 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:50 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:50 INFO::Reading weights for chromosome 19
#> 2021-06-06 20:20:50 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:50 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:50 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:50 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:50 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:50 INFO::Imputation done: number of genes with imputed expression: 1 for chr 19
#> 2021-06-06 20:20:50 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:50 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:50 INFO::Reading weights for chromosome 20
#> 2021-06-06 20:20:50 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:50 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:50 INFO::Start gene z score imputation ...
#> 2021-06-06 20:20:50 INFO::ld genotype is given, using genotypes to impute gene z score.
#> 2021-06-06 20:20:50 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:50 INFO::Imputation done: number of genes with imputed expression: 1 for chr 20
#> 2021-06-06 20:20:50 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:50 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:50 INFO::Reading weights for chromosome 21
#> 2021-06-06 20:20:50 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:50 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:50 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:50 INFO::Imputation done: number of genes with imputed expression: 0 for chr 21
#> 2021-06-06 20:20:50 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:50 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:20:50 INFO::Reading weights for chromosome 22
#> 2021-06-06 20:20:50 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:20:50 INFO::collecting gene weight information ...
#> 2021-06-06 20:20:50 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:20:50 INFO::Imputation done: number of genes with imputed expression: 0 for chr 22
z_gene <- res$z_gene
ld_exprfs <- res$ld_exprfs

head(z_gene)
#>          id          z
#> gene1 gene1 18.4624001
#> gene2 gene2 -1.7169590
#> gene3 gene3 -0.4940096
#> gene4 gene4 -0.5448988
#> gene5 gene5 -0.1813120
#> gene6 gene6 -0.2698443
ld_exprfs
#>  [1] "~/temp/test_ss_chr1.expr.gz"  "~/temp/test_ss_chr2.expr.gz" 
#>  [3] "~/temp/test_ss_chr3.expr.gz"  "~/temp/test_ss_chr4.expr.gz" 
#>  [5] "~/temp/test_ss_chr5.expr.gz"  "~/temp/test_ss_chr6.expr.gz" 
#>  [7] "~/temp/test_ss_chr7.expr.gz"  "~/temp/test_ss_chr8.expr.gz" 
#>  [9] "~/temp/test_ss_chr9.expr.gz"  "~/temp/test_ss_chr10.expr.gz"
#> [11] "~/temp/test_ss_chr11.expr.gz" "~/temp/test_ss_chr12.expr.gz"
#> [13] "~/temp/test_ss_chr13.expr.gz" "~/temp/test_ss_chr14.expr.gz"
#> [15] "~/temp/test_ss_chr15.expr.gz" "~/temp/test_ss_chr16.expr.gz"
#> [17] "~/temp/test_ss_chr17.expr.gz" "~/temp/test_ss_chr18.expr.gz"
#> [19] "~/temp/test_ss_chr19.expr.gz" "~/temp/test_ss_chr20.expr.gz"
#> [21] "~/temp/test_ss_chr21.expr.gz" "~/temp/test_ss_chr22.expr.gz"

The weight parameter from impute_expr_z function can take either FUSION TWAS format or predictdb format. In the above example, we provided .db format weights. you can also use weight = weight.fusion if you provide fusion format.

Step 2: run ctwas_rss.

In this step we will perform the causal TWAS algorithm, the algorithm will run susie iteratively for parameter estimation and lastly provide PIPs for all genes and SNPs included in the analysis. If you don’t want to define your own LD regions, then you can use the one defined by ldetect by simply specifying the population name using the ld_regions argument and specify the genome build that matches your data by using ld_region_version. Currently, only genome build b37 and b38 are provide by the package. If you need other versions, please download the regions file from the package source, liftover and use the ld_regions_custom to provide your own customed regions. One feature of the ctwas function is that it allows parallel computing. You can specify number of cores to use by the ncore argument.

Run ctwas_rss using genotype data for LD with the custom regions.

pars <- ctwas_rss(z_gene = z_gene, z_snp = z_snp, ld_exprfs = ld_exprfs, ld_pgenfs = ld_pgenfs, ld_regions_custom = regionsfile, thin = 0.9, max_snp_region = 20, outputdir = outputdir, outname = "test_ss", ncore = 1)
#> 2021-06-06 20:20:50 INFO::ctwas started ...
#> 2021-06-06 20:20:50 INFO::LD region file: /tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_regions.bed
#> 2021-06-06 20:20:50 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:20:51 INFO::No. LD regions: 374
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr1: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr1
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr1 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr2: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr2
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr2 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr3: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr3
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr3 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr4: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr4
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr4 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr5: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr5
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr5 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr6: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr6
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr6 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr7: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr7
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr7 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr8: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr8
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr8 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr9: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr9
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr9 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr10: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr10
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr10 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr11: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr11
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr11 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr12: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr12
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr12 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr13: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr13
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr13 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr14: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr14
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr14 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr15: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr15
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr15 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr16: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr16 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr17: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr17 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr18: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr18
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr18 after merging: 16
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr19: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr19
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr19 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr20: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr20
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr20 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr21: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr21
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr21 after merging: 17
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr22: 17
#> 2021-06-06 20:20:51 INFO::Merge regions for chr22
#> 2021-06-06 20:20:51 INFO::No. regions with at least one SNP/gene for chr22 after merging: 17
#> 2021-06-06 20:20:51 INFO::Trim regions with SNPs more than Inf
#> 2021-06-06 20:20:51 INFO::Run susie iteratively, getting rough estimate ...
#> 2021-06-06 20:20:51 INFO::run iteration 1
#> 2021-06-06 20:21:01 INFO::After iteration 1, gene prior 0.15396123828879:, SNP prior:0.0375465027367227
#> 2021-06-06 20:21:01 INFO::run iteration 2
#> 2021-06-06 20:21:10 INFO::After iteration 2, gene prior 0.0882628368144239:, SNP prior:0.036038024148467
#> 2021-06-06 20:21:10 INFO::run iteration 3
#> 2021-06-06 20:21:20 INFO::After iteration 3, gene prior 0.0666064574425756:, SNP prior:0.0353425325788066
#> 2021-06-06 20:21:20 INFO::Blocks are filtered: 16 blocks left
#> 2021-06-06 20:21:20 INFO::Run susie iteratively, getting accurate estimate ...
#> 2021-06-06 20:21:20 INFO::run iteration 1
#> 2021-06-06 20:21:23 INFO::After iteration 1, gene prior 0.120355586178647:, SNP prior:0.0313256408158838
#> 2021-06-06 20:21:23 INFO::run iteration 2
#> 2021-06-06 20:21:27 INFO::After iteration 2, gene prior 0.12667132575201:, SNP prior:0.0298158963892055
#> 2021-06-06 20:21:27 INFO::run iteration 3
#> 2021-06-06 20:21:30 INFO::After iteration 3, gene prior 0.127968853532249:, SNP prior:0.0290029393354405
#> 2021-06-06 20:21:30 INFO::run iteration 4
#> 2021-06-06 20:21:33 INFO::After iteration 4, gene prior 0.128271093620608:, SNP prior:0.0285009525352293
#> 2021-06-06 20:21:33 INFO::run iteration 5
#> 2021-06-06 20:21:36 INFO::After iteration 5, gene prior 0.128362212919622:, SNP prior:0.0281705166835134
#> 2021-06-06 20:21:36 INFO::run iteration 6
#> 2021-06-06 20:21:40 INFO::After iteration 6, gene prior 0.128404928383303:, SNP prior:0.0279473157995042
#> 2021-06-06 20:21:40 INFO::run iteration 7
#> 2021-06-06 20:21:43 INFO::After iteration 7, gene prior 0.128433872355419:, SNP prior:0.0277968216319189
#> 2021-06-06 20:21:43 INFO::run iteration 8
#> 2021-06-06 20:21:46 INFO::After iteration 8, gene prior 0.128456950378025:, SNP prior:0.0276985170674317
#> 2021-06-06 20:21:46 INFO::run iteration 9
#> 2021-06-06 20:21:50 INFO::After iteration 9, gene prior 0.128476380690168:, SNP prior:0.0276393026884226
#> 2021-06-06 20:21:50 INFO::run iteration 10
#> 2021-06-06 20:21:53 INFO::After iteration 10, gene prior 0.128493081469183:, SNP prior:0.0276103551493495
#> 2021-06-06 20:21:53 INFO::run iteration 11
#> 2021-06-06 20:21:56 INFO::After iteration 11, gene prior 0.128507602618013:, SNP prior:0.0276054787861449
#> 2021-06-06 20:21:56 INFO::run iteration 12
#> 2021-06-06 20:21:59 INFO::After iteration 12, gene prior 0.128520339836496:, SNP prior:0.027620175268022
#> 2021-06-06 20:21:59 INFO::run iteration 13
#> 2021-06-06 20:22:02 INFO::After iteration 13, gene prior 0.128531596834339:, SNP prior:0.0276510876535266
#> 2021-06-06 20:22:02 INFO::run iteration 14
#> 2021-06-06 20:22:05 INFO::After iteration 14, gene prior 0.128541612372474:, SNP prior:0.0276956524931885
#> 2021-06-06 20:22:05 INFO::run iteration 15
#> 2021-06-06 20:22:09 INFO::After iteration 15, gene prior 0.128550576663043:, SNP prior:0.027751873640597
#> 2021-06-06 20:22:09 INFO::run iteration 16
#> 2021-06-06 20:22:12 INFO::After iteration 16, gene prior 0.128558642873839:, SNP prior:0.0278181703671663
#> 2021-06-06 20:22:12 INFO::run iteration 17
#> 2021-06-06 20:22:15 INFO::After iteration 17, gene prior 0.128565935576112:, SNP prior:0.0278932725078092
#> 2021-06-06 20:22:15 INFO::run iteration 18
#> 2021-06-06 20:22:18 INFO::After iteration 18, gene prior 0.128572557045647:, SNP prior:0.0279761463137103
#> 2021-06-06 20:22:18 INFO::run iteration 19
#> 2021-06-06 20:22:21 INFO::After iteration 19, gene prior 0.128578592006233:, SNP prior:0.0280659409031903
#> 2021-06-06 20:22:21 INFO::run iteration 20
#> 2021-06-06 20:22:25 INFO::After iteration 20, gene prior 0.128584111231672:, SNP prior:0.0281619488616015
#> 2021-06-06 20:22:25 INFO::run iteration 21
#> 2021-06-06 20:22:28 INFO::After iteration 21, gene prior 0.1285891743058:, SNP prior:0.0282635767670036
#> 2021-06-06 20:22:28 INFO::run iteration 22
#> 2021-06-06 20:22:31 INFO::After iteration 22, gene prior 0.12859383175696:, SNP prior:0.0283703228109487
#> 2021-06-06 20:22:31 INFO::run iteration 23
#> 2021-06-06 20:22:34 INFO::After iteration 23, gene prior 0.128598126723934:, SNP prior:0.0284817595772505
#> 2021-06-06 20:22:34 INFO::run iteration 24
#> 2021-06-06 20:22:37 INFO::After iteration 24, gene prior 0.128602096267853:, SNP prior:0.0285975206280642
#> 2021-06-06 20:22:37 INFO::run iteration 25
#> 2021-06-06 20:22:41 INFO::After iteration 25, gene prior 0.128605772414218:, SNP prior:0.0287172899394402
#> 2021-06-06 20:22:41 INFO::run iteration 26
#> 2021-06-06 20:22:44 INFO::After iteration 26, gene prior 0.128609182987344:, SNP prior:0.0288407934965567
#> 2021-06-06 20:22:44 INFO::run iteration 27
#> 2021-06-06 20:22:47 INFO::After iteration 27, gene prior 0.128612352283772:, SNP prior:0.02896779254484
#> 2021-06-06 20:22:47 INFO::run iteration 28
#> 2021-06-06 20:22:51 INFO::After iteration 28, gene prior 0.128615301619717:, SNP prior:0.0290980781242551
#> 2021-06-06 20:22:51 INFO::run iteration 29
#> 2021-06-06 20:22:54 INFO::After iteration 29, gene prior 0.128618049779188:, SNP prior:0.0292314666077282
#> 2021-06-06 20:22:54 INFO::run iteration 30
#> 2021-06-06 20:22:57 INFO::After iteration 30, gene prior 0.128620613383212:, SNP prior:0.0293677960324969
#> 2021-06-06 20:22:57 INFO::Run susie for all regions.
#> 2021-06-06 20:22:57 INFO::run iteration 1
#> 2021-06-06 20:23:08 INFO::After iteration 1, gene prior 0.113965838175324:, SNP prior:0.138489926675174
#> 2021-06-06 20:23:08 INFO::No. LD regions: 374
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr1: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr1
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr1 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr2: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr2
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr2 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr3: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr3
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr3 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr4: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr4
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr4 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr5: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr5
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr5 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr6: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr6
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr6 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr7: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr7
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr7 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr8: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr8
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr8 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr9: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr9
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr9 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr10: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr10
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr10 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr11: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr11
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr11 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr12: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr12
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr12 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr13: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr13
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr13 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr14: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr14
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr14 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr15: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr15
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr15 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr16: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr16 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr17: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr17 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr18: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr18
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr18 after merging: 16
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr19: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr19
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr19 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr20: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr20
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr20 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr21: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr21
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr21 after merging: 17
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr22: 17
#> 2021-06-06 20:23:08 INFO::Merge regions for chr22
#> 2021-06-06 20:23:08 INFO::No. regions with at least one SNP/gene for chr22 after merging: 17
#> 2021-06-06 20:23:08 INFO::Trim regions with SNPs more than 20
#> 2021-06-06 20:23:09 INFO::Number of regions that contains strong gene signals: 1
#> 2021-06-06 20:23:09 INFO::Rerun susie for regions with strong gene signals using full SNPs.
#> 2021-06-06 20:23:09 INFO::run iteration 1
#> 2021-06-06 20:23:11 INFO::After iteration 1, gene prior 1:, SNP prior:0.122064022562353

Steps when using LD R (correlation) matrices

Step 1: impute expression z scores

In the following code, we compute expression z scores for each gene. We do this chromosome by chromosome. impute_expr_z will return z scores for each gene. It will also generate files containing the reference LD genotypes for eQTLs in the gene and the filename is returned. We will need these files in the ctwas_rss function in step 2.

When providing LD R matrices:

outputdir <- "~/temp"

res <- impute_expr_z(z_snp = z_snp,  weight = weight.predictdb, ld_R_dir = ld_R_dir,
                           method = "lasso", outputdir = outputdir,
                           outname = "test_ss")
#> 2021-06-06 20:23:12 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:12 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:12 INFO::Reading weights for chromosome 1
#> 2021-06-06 20:23:12 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:12 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:12 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:12 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:12 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:12 INFO::Imputation done: number of genes with imputed expression: 1 for chr 1
#> 2021-06-06 20:23:12 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:12 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:12 INFO::Reading weights for chromosome 2
#> 2021-06-06 20:23:12 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:12 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:12 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:12 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:12 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:12 INFO::Imputation done: number of genes with imputed expression: 1 for chr 2
#> 2021-06-06 20:23:12 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:12 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:12 INFO::Reading weights for chromosome 3
#> 2021-06-06 20:23:12 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:12 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:12 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:12 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 3
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 4
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 4
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 5
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 5
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 6
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 6
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 7
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 7
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 8
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 8
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 9
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:13 INFO::Imputation done: number of genes with imputed expression: 1 for chr 9
#> 2021-06-06 20:23:13 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:13 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:13 INFO::Reading weights for chromosome 10
#> 2021-06-06 20:23:13 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:13 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:13 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:13 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:13 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 10
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 11
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 11
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 12
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 12
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 13
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 13
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 14
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 14
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 15
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 15
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 16
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:14 INFO::Imputation done: number of genes with imputed expression: 1 for chr 16
#> 2021-06-06 20:23:14 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:14 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:14 INFO::Reading weights for chromosome 17
#> 2021-06-06 20:23:14 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:14 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:14 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:14 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:14 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:15 INFO::Imputation done: number of genes with imputed expression: 1 for chr 17
#> 2021-06-06 20:23:15 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:15 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:15 INFO::Reading weights for chromosome 18
#> 2021-06-06 20:23:15 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:15 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:15 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:15 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:15 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:15 INFO::Imputation done: number of genes with imputed expression: 1 for chr 18
#> 2021-06-06 20:23:15 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:15 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:15 INFO::Reading weights for chromosome 19
#> 2021-06-06 20:23:15 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:15 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:15 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:15 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:15 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:15 INFO::Imputation done: number of genes with imputed expression: 1 for chr 19
#> 2021-06-06 20:23:15 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:15 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:15 INFO::Reading weights for chromosome 20
#> 2021-06-06 20:23:15 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:15 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:15 INFO::Start gene z score imputation ...
#> 2021-06-06 20:23:15 INFO::Using given LD matrices to impute gene z score.
#> 2021-06-06 20:23:15 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:15 INFO::Imputation done: number of genes with imputed expression: 1 for chr 20
#> 2021-06-06 20:23:15 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:15 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:15 INFO::Reading weights for chromosome 21
#> 2021-06-06 20:23:15 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:15 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:15 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:15 INFO::Imputation done: number of genes with imputed expression: 0 for chr 21
#> 2021-06-06 20:23:15 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:15 INFO::will also flip weights to match LD reference for each gene
#> 2021-06-06 20:23:15 INFO::Reading weights for chromosome 22
#> 2021-06-06 20:23:15 INFO::number of genes with weights provided: 20
#> 2021-06-06 20:23:15 INFO::collecting gene weight information ...
#> 2021-06-06 20:23:15 INFO::Imputation done, writing results to output...
#> 2021-06-06 20:23:15 INFO::Imputation done: number of genes with imputed expression: 0 for chr 22
z_gene <- res$z_gene
ld_exprfs <- res$ld_exprfs

head(z_gene)
#>          id          z
#> gene1 gene1 18.4624001
#> gene2 gene2 -1.7169590
#> gene3 gene3 -0.4940096
#> gene4 gene4 -0.5310785
#> gene5 gene5 -0.1813120
#> gene6 gene6 -0.2698443
ld_exprfs
#>  [1] "~/temp/test_ss_chr1.expr.gz"  "~/temp/test_ss_chr2.expr.gz" 
#>  [3] "~/temp/test_ss_chr3.expr.gz"  "~/temp/test_ss_chr4.expr.gz" 
#>  [5] "~/temp/test_ss_chr5.expr.gz"  "~/temp/test_ss_chr6.expr.gz" 
#>  [7] "~/temp/test_ss_chr7.expr.gz"  "~/temp/test_ss_chr8.expr.gz" 
#>  [9] "~/temp/test_ss_chr9.expr.gz"  "~/temp/test_ss_chr10.expr.gz"
#> [11] "~/temp/test_ss_chr11.expr.gz" "~/temp/test_ss_chr12.expr.gz"
#> [13] "~/temp/test_ss_chr13.expr.gz" "~/temp/test_ss_chr14.expr.gz"
#> [15] "~/temp/test_ss_chr15.expr.gz" "~/temp/test_ss_chr16.expr.gz"
#> [17] "~/temp/test_ss_chr17.expr.gz" "~/temp/test_ss_chr18.expr.gz"
#> [19] "~/temp/test_ss_chr19.expr.gz" "~/temp/test_ss_chr20.expr.gz"
#> [21] "~/temp/test_ss_chr21.expr.gz" "~/temp/test_ss_chr22.expr.gz"

The weight parameter from impute_expr_z function can take either FUSION TWAS format or predictdb format.

Step 2: run ctwas_rss.

In this step we will perform the causal TWAS algorithm, the algorithm will run susie iteratively for parameter estimation and lastly provide PIPs for all genes and SNPs included in the analysis. If you don’t want to define your own LD regions, then you can use the one defined by ldetect by simply specifying the population name using the ld_regions argument and specify the genome build that matches your data by using ld_region_version. Currently, only genome build b37 and b38 are provide by the package. If you need other versions, please download the regions file from the package source, liftover and use the ld_regions_custom to provide your own customed regions. . One feature of the ctwas function is that it allows parallel computing. You can specify number of cores to use by the ncore argument.

When providing LD R matrices, use the following command:

pars <- ctwas_rss(z_gene = z_gene, z_snp = z_snp, ld_exprfs = ld_exprfs, ld_R_dir = ld_R_dir, ld_regions_custom = regionsfile, thin = 0.9, max_snp_region = 20, outputdir = outputdir, outname = "test_ss", ncore = 1)
#> 2021-06-06 20:23:15 INFO::ctwas started ...
#> 2021-06-06 20:23:17 INFO::LD region file: /tmp/RtmpzZjoXP/temp_libpath6eb4bb559ef/ctwas/extdata/example_regions.bed
#> 2021-06-06 20:23:17 INFO::flipping z scores to match LD reference
#> 2021-06-06 20:23:17 INFO::No. LD regions: 374
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr1: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr1
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr1 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr2: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr2
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr2 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr3: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr3
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr3 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr4: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr4
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr4 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr5: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr5
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr5 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr6: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr6
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr6 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr7: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr7
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr7 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr8: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr8
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr8 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr9: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr9
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr9 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr10: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr10
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr10 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr11: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr11
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr11 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr12: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr12
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr12 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr13: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr13
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr13 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr14: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr14
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr14 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr15: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr15
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr15 after merging: 16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr16: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr16
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr16 after merging: 17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr17: 17
#> 2021-06-06 20:23:17 INFO::Merge regions for chr17
#> 2021-06-06 20:23:17 INFO::No. regions with at least one SNP/gene for chr17 after merging: 16
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr18: 17
#> 2021-06-06 20:23:18 INFO::Merge regions for chr18
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr18 after merging: 16
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr19: 17
#> 2021-06-06 20:23:18 INFO::Merge regions for chr19
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr19 after merging: 17
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr20: 17
#> 2021-06-06 20:23:18 INFO::Merge regions for chr20
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr20 after merging: 17
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr21: 17
#> 2021-06-06 20:23:18 INFO::Merge regions for chr21
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr21 after merging: 17
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr22: 17
#> 2021-06-06 20:23:18 INFO::Merge regions for chr22
#> 2021-06-06 20:23:18 INFO::No. regions with at least one SNP/gene for chr22 after merging: 17
#> 2021-06-06 20:23:18 INFO::Trim regions with SNPs more than Inf
#> 2021-06-06 20:23:18 INFO::Adding R matrix info, as genotype is not given
#> 2021-06-06 20:23:18 INFO::Adding R matrix info for chrom 1
#> 2021-06-06 20:23:18 INFO::Adding R matrix info for chrom 2
#> 2021-06-06 20:23:18 INFO::Adding R matrix info for chrom 3
#> 2021-06-06 20:23:18 INFO::Adding R matrix info for chrom 4
#> 2021-06-06 20:23:19 INFO::Adding R matrix info for chrom 5
#> 2021-06-06 20:23:19 INFO::Adding R matrix info for chrom 6
#> 2021-06-06 20:23:19 INFO::Adding R matrix info for chrom 7
#> 2021-06-06 20:23:19 INFO::Adding R matrix info for chrom 8
#> 2021-06-06 20:23:19 INFO::Adding R matrix info for chrom 9
#> 2021-06-06 20:23:20 INFO::Adding R matrix info for chrom 10
#> 2021-06-06 20:23:20 INFO::Adding R matrix info for chrom 11
#> 2021-06-06 20:23:20 INFO::Adding R matrix info for chrom 12
#> 2021-06-06 20:23:20 INFO::Adding R matrix info for chrom 13
#> 2021-06-06 20:23:21 INFO::Adding R matrix info for chrom 14
#> 2021-06-06 20:23:21 INFO::Adding R matrix info for chrom 15
#> 2021-06-06 20:23:21 INFO::Adding R matrix info for chrom 16
#> 2021-06-06 20:23:21 INFO::Adding R matrix info for chrom 17
#> 2021-06-06 20:23:21 INFO::Adding R matrix info for chrom 18
#> 2021-06-06 20:23:21 INFO::Adding R matrix info for chrom 19
#> 2021-06-06 20:23:22 INFO::Adding R matrix info for chrom 20
#> 2021-06-06 20:23:22 INFO::Adding R matrix info for chrom 21
#> 2021-06-06 20:23:22 INFO::Adding R matrix info for chrom 22
#> 2021-06-06 20:23:22 INFO::Run susie iteratively, getting rough estimate ...
#> 2021-06-06 20:23:22 INFO::run iteration 1
#> 2021-06-06 20:23:28 INFO::After iteration 1, gene prior 0.182071436616161:, SNP prior:0.03746440484442
#> 2021-06-06 20:23:28 INFO::run iteration 2
#> 2021-06-06 20:23:33 INFO::After iteration 2, gene prior 0.105026103339378:, SNP prior:0.0358910262612177
#> 2021-06-06 20:23:33 INFO::run iteration 3
#> 2021-06-06 20:23:39 INFO::After iteration 3, gene prior 0.0754309047368718:, SNP prior:0.035204891354971
#> 2021-06-06 20:23:39 INFO::Blocks are filtered: 11 blocks left
#> 2021-06-06 20:23:39 INFO::Run susie iteratively, getting accurate estimate ...
#> 2021-06-06 20:23:39 INFO::run iteration 1
#> 2021-06-06 20:23:41 INFO::After iteration 1, gene prior 0.122501316616366:, SNP prior:0.0245617877007893
#> 2021-06-06 20:23:41 INFO::run iteration 2
#> 2021-06-06 20:23:43 INFO::After iteration 2, gene prior 0.127432210427403:, SNP prior:0.0213106631766159
#> 2021-06-06 20:23:43 INFO::run iteration 3
#> 2021-06-06 20:23:45 INFO::After iteration 3, gene prior 0.128448878753661:, SNP prior:0.0195333397842195
#> 2021-06-06 20:23:45 INFO::run iteration 4
#> 2021-06-06 20:23:46 INFO::After iteration 4, gene prior 0.128678110322031:, SNP prior:0.018348890367087
#> 2021-06-06 20:23:46 INFO::run iteration 5
#> 2021-06-06 20:23:48 INFO::After iteration 5, gene prior 0.128738567338109:, SNP prior:0.0174748856892675
#> 2021-06-06 20:23:48 INFO::run iteration 6
#> 2021-06-06 20:23:50 INFO::After iteration 6, gene prior 0.128761361933475:, SNP prior:0.0167891912481062
#> 2021-06-06 20:23:50 INFO::run iteration 7
#> 2021-06-06 20:23:52 INFO::After iteration 7, gene prior 0.128774649289401:, SNP prior:0.0162287910244394
#> 2021-06-06 20:23:52 INFO::run iteration 8
#> 2021-06-06 20:23:54 INFO::After iteration 8, gene prior 0.128784674884784:, SNP prior:0.0157571957780964
#> 2021-06-06 20:23:54 INFO::run iteration 9
#> 2021-06-06 20:23:56 INFO::After iteration 9, gene prior 0.128792991292574:, SNP prior:0.015351509831133
#> 2021-06-06 20:23:56 INFO::run iteration 10
#> 2021-06-06 20:23:57 INFO::After iteration 10, gene prior 0.12880011521214:, SNP prior:0.0149964922625565
#> 2021-06-06 20:23:57 INFO::run iteration 11
#> 2021-06-06 20:23:59 INFO::After iteration 11, gene prior 0.128806306784461:, SNP prior:0.0146815251890626
#> 2021-06-06 20:23:59 INFO::run iteration 12
#> 2021-06-06 20:24:01 INFO::After iteration 12, gene prior 0.128811740650231:, SNP prior:0.0143989381832964
#> 2021-06-06 20:24:01 INFO::run iteration 13
#> 2021-06-06 20:24:03 INFO::After iteration 13, gene prior 0.128816548431023:, SNP prior:0.0141430232770228
#> 2021-06-06 20:24:03 INFO::run iteration 14
#> 2021-06-06 20:24:05 INFO::After iteration 14, gene prior 0.128820833207632:, SNP prior:0.013909426533009
#> 2021-06-06 20:24:05 INFO::run iteration 15
#> 2021-06-06 20:24:07 INFO::After iteration 15, gene prior 0.128824676994723:, SNP prior:0.0136947565678126
#> 2021-06-06 20:24:07 INFO::run iteration 16
#> 2021-06-06 20:24:09 INFO::After iteration 16, gene prior 0.12882814573029:, SNP prior:0.0134963238712331
#> 2021-06-06 20:24:09 INFO::run iteration 17
#> 2021-06-06 20:24:10 INFO::After iteration 17, gene prior 0.128831292922053:, SNP prior:0.0133119620572124
#> 2021-06-06 20:24:10 INFO::run iteration 18
#> 2021-06-06 20:24:12 INFO::After iteration 18, gene prior 0.12883416238073:, SNP prior:0.0131399021469813
#> 2021-06-06 20:24:12 INFO::run iteration 19
#> 2021-06-06 20:24:14 INFO::After iteration 19, gene prior 0.128836790287537:, SNP prior:0.0129786821688169
#> 2021-06-06 20:24:14 INFO::run iteration 20
#> 2021-06-06 20:24:16 INFO::After iteration 20, gene prior 0.128839206768737:, SNP prior:0.0128270808703129
#> 2021-06-06 20:24:16 INFO::run iteration 21
#> 2021-06-06 20:24:18 INFO::After iteration 21, gene prior 0.128841437103527:, SNP prior:0.0126840682603186
#> 2021-06-06 20:24:18 INFO::run iteration 22
#> 2021-06-06 20:24:20 INFO::After iteration 22, gene prior 0.128843502658103:, SNP prior:0.0125487681301962
#> 2021-06-06 20:24:20 INFO::run iteration 23
#> 2021-06-06 20:24:22 INFO::After iteration 23, gene prior 0.128845421613988:, SNP prior:0.0124204292532766
#> 2021-06-06 20:24:22 INFO::run iteration 24
#> 2021-06-06 20:24:24 INFO::After iteration 24, gene prior 0.128847209540713:, SNP prior:0.0122984029716046
#> 2021-06-06 20:24:24 INFO::run iteration 25
#> 2021-06-06 20:24:26 INFO::After iteration 25, gene prior 0.128848879849829:, SNP prior:0.0121821255519112
#> 2021-06-06 20:24:26 INFO::run iteration 26
#> 2021-06-06 20:24:27 INFO::After iteration 26, gene prior 0.128850444157677:, SNP prior:0.0120711041495756
#> 2021-06-06 20:24:27 INFO::run iteration 27
#> 2021-06-06 20:24:29 INFO::After iteration 27, gene prior 0.128851912577487:, SNP prior:0.0119649055349574
#> 2021-06-06 20:24:29 INFO::run iteration 28
#> 2021-06-06 20:24:31 INFO::After iteration 28, gene prior 0.128853293956248:, SNP prior:0.0118631469580429
#> 2021-06-06 20:24:31 INFO::run iteration 29
#> 2021-06-06 20:24:33 INFO::After iteration 29, gene prior 0.128854596068133:, SNP prior:0.0117654886851768
#> 2021-06-06 20:24:33 INFO::run iteration 30
#> 2021-06-06 20:24:35 INFO::After iteration 30, gene prior 0.128855825773462:, SNP prior:0.0116716278555945
#> 2021-06-06 20:24:35 INFO::Run susie for all regions.
#> 2021-06-06 20:24:35 INFO::run iteration 1
#> 2021-06-06 20:24:41 INFO::After iteration 1, gene prior 0.119773952871303:, SNP prior:0.0594530674263792
#> 2021-06-06 20:24:41 INFO::No. LD regions: 374
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr1: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr1
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr1 after merging: 17
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr2: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr2
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr2 after merging: 17
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr3: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr3
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr3 after merging: 17
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr4: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr4
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr4 after merging: 16
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr5: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr5
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr5 after merging: 17
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr6: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr6
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr6 after merging: 17
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr7: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr7
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr7 after merging: 16
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr8: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr8
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr8 after merging: 16
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr9: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr9
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr9 after merging: 17
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr10: 17
#> 2021-06-06 20:24:41 INFO::Merge regions for chr10
#> 2021-06-06 20:24:41 INFO::No. regions with at least one SNP/gene for chr10 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr11: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr11
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr11 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr12: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr12
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr12 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr13: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr13
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr13 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr14: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr14
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr14 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr15: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr15
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr15 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr16: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr16 after merging: 17
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr17: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr17
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr17 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr18: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr18
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr18 after merging: 16
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr19: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr19
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr19 after merging: 17
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr20: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr20
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr20 after merging: 17
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr21: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr21
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr21 after merging: 17
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr22: 17
#> 2021-06-06 20:24:42 INFO::Merge regions for chr22
#> 2021-06-06 20:24:42 INFO::No. regions with at least one SNP/gene for chr22 after merging: 17
#> 2021-06-06 20:24:42 INFO::Trim regions with SNPs more than 20
#> 2021-06-06 20:24:42 INFO::Adding R matrix info, as genotype is not given
#> 2021-06-06 20:24:42 INFO::Adding R matrix info for chrom 1
#> 2021-06-06 20:24:43 INFO::Adding R matrix info for chrom 2
#> 2021-06-06 20:24:43 INFO::Adding R matrix info for chrom 3
#> 2021-06-06 20:24:43 INFO::Adding R matrix info for chrom 4
#> 2021-06-06 20:24:43 INFO::Adding R matrix info for chrom 5
#> 2021-06-06 20:24:43 INFO::Adding R matrix info for chrom 6
#> 2021-06-06 20:24:43 INFO::Adding R matrix info for chrom 7
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 8
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 9
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 10
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 11
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 12
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 13
#> 2021-06-06 20:24:44 INFO::Adding R matrix info for chrom 14
#> 2021-06-06 20:24:45 INFO::Adding R matrix info for chrom 15
#> 2021-06-06 20:24:45 INFO::Adding R matrix info for chrom 16
#> 2021-06-06 20:24:45 INFO::Adding R matrix info for chrom 17
#> 2021-06-06 20:24:45 INFO::Adding R matrix info for chrom 18
#> 2021-06-06 20:24:45 INFO::Adding R matrix info for chrom 19
#> 2021-06-06 20:24:46 INFO::Adding R matrix info for chrom 20
#> 2021-06-06 20:24:46 INFO::Adding R matrix info for chrom 21
#> 2021-06-06 20:24:46 INFO::Adding R matrix info for chrom 22
#> 2021-06-06 20:24:47 INFO::Number of regions that contains strong gene signals: 1
#> 2021-06-06 20:24:47 INFO::Rerun susie for regions with strong gene signals using full SNPs.
#> 2021-06-06 20:24:47 INFO::run iteration 1
#> 2021-06-06 20:24:48 INFO::After iteration 1, gene prior 1:, SNP prior:0.0500167713726725

Output from ctwas

ctwas_rss returns the estimated parameters (in the order of gene, SNP):

pars
#> $group_prior
#> [1] 0.12885583 0.01050447
#> 
#> $group_prior_var
#> [1] 293.2541393   0.2139139

PIP results are given in outname.susieIrss.txt. This file contains PIP for each gene and SNP, please check susie_pip column.

head(data.table::fread("~/temp/test_ss.susieIrss.txt"))[, c("chrom", "id", "pos", "type", "susie_pip")]
#>    chrom     id pos type   susie_pip
#> 1:     1 rs1005   5  SNP 0.757344613
#> 2:     1 rs1006   6  SNP 0.002614718
#> 3:     1 rs1007   7  SNP 0.774952213
#> 4:     1 rs1008   8  SNP 0.002434619
#> 5:     1 rs1009   9  SNP 0.748231244
#> 6:     1 rs1010  10  SNP 0.686410868

You will also notice a few auxilary/intermediate files produced by ctwas_rss. These files can be useful for diagnoiss..Rd files contains parameters estimation updates with each iteration.

sessionInfo()
#> R version 3.6.1 (2019-07-05)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Scientific Linux 7.4 (Nitrogen)
#> 
#> Matrix products: default
#> BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ctwas_0.1.28
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.5        compiler_3.6.1    pillar_1.5.1     
#>  [4] iterators_1.0.10  tools_3.6.1       bit_1.1-14       
#>  [7] digest_0.6.20     RSQLite_2.2.7     debugme_1.1.0    
#> [10] evaluate_0.14     memoise_1.1.0     lifecycle_1.0.0  
#> [13] tibble_3.1.0      lattice_0.20-38   pkgconfig_2.0.2  
#> [16] rlang_0.4.10      Matrix_1.2-18     foreach_1.4.4    
#> [19] DBI_1.1.0         parallel_3.6.1    yaml_2.2.0       
#> [22] pkgdown_1.6.1     xfun_0.8          stringr_1.4.0    
#> [25] dplyr_1.0.5       knitr_1.23        pgenlibr_0.2     
#> [28] desc_1.2.0        generics_0.0.2    fs_1.3.1         
#> [31] vctrs_0.3.7       bit64_0.9-7       tidyselect_1.1.0 
#> [34] rprojroot_1.3-2   grid_3.6.1        glue_1.4.2       
#> [37] data.table_1.13.2 R6_2.4.0          fansi_0.4.0      
#> [40] rmarkdown_1.13    blob_1.2.0        purrr_0.3.4      
#> [43] magrittr_1.5      backports_1.1.4   codetools_0.2-16 
#> [46] htmltools_0.3.6   ellipsis_0.2.0.1  assertthat_0.2.1 
#> [49] logging_0.10-108  utf8_1.1.4        stringi_1.4.3    
#> [52] doParallel_1.0.14 crayon_1.3.4