This reproducible R Markdown analysis was created with workflowr (version 1.6.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20190823) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
Here we describe an end-to-end RSS-NET analysis of inflammatory bowel disease (IBD) GWAS summary statistics (Liu et al, 2015) and a gene regulatory network inferred for natural killer (NK) cells. This example illustrates the actual data analyses performed in Zhu et al (2019).
To reproduce results of this example, please use scripts in the directory script_dir, and follow the step-by-step guide below. Before running any script in script_dir, please install RSS-NET.
Since a real genome-wide analysis is conducted here, this example is more complicated than the previous simulation example. It is advisable to go through the previous simulation example before diving into this real data example.
Note that the working directory here is assumed to be wdtba. Please modify scripts accordingly if a different directory is used.
Step-by-step illustration
Download input data files
1. ibd2015_sumstat.mat: processed GWAS summary statistics and LD matrix estimates
This file is large (43G) because it has a LD matrix of 1.1 million common SNPs. Please contact me (xiangzhu[at]stanford.edu) if you have trouble accessing this file.
$ md5sum ibd2015_sumstat.matad1763079ee7e46b21722f74e037a230 ibd2015_sumstat.mat$ du -sh ibd2015_sumstat.mat 43G ibd2015_sumstat.mat
Let’s look at the contents of ibd2015_sumstat.mat.
GWAS summary statistics and LD estimates are stored as cell arrays. RSS-NET only uses the following variables:
betahat{j,1}, single-SNP effect size estimates of all SNPs on chromosome j;
se{j,1}, standard errors of betahat{j, 1};
chr{j,1} and pos{j, 1}, physical positions of these SNPs (GRCh37 build);
SiRiS{j,1}, a sparse matrix, defined as repmat((1./se),1,p) .* R .* repmat((1./se)',p,1), where R is the estimated LD matrix of these p SNPs.
2. ibd2015_snp2gene.mat: physical distance between SNPs and genes
This file contains the physical distance between each GWAS SNP and each protein-coding gene, within 1 Mb. This file corresponds to \({\bf G}_j\) in the RSS-NET model.
$ md5sum ibd2015_snp2gene.mat 7832838e2675e4cf3b85f471fed95554 ibd2015_snp2gene.mat$ du -sh ibd2015_snp2gene.mat 224M ibd2015_snp2gene.mat
In this example, there are 18334 SNPs and 1081481 genes.
The SNP-to-gene distance information is captured by a three-column matrix [colid rowid val]. For example, the distance between gene 1 and SNP 6 is 978947 bps.
For a given network, transcription factors (TFs) are stored in rowid and target genes (TGs) are stored in colid. In this example there are 3105 TGs and 376 TFs. Among these TFs and TGs, there are 92399 edges. The edge weights range from 0.61 to 1. These TF-to-TG connections and edge weights correspond to \(\{{\bf T}_g,v_{gt}\}\) in the RSS-NET model.
This file contains binary annotations whether a SNP is “near” the given network, that is, within 100 kb of any network element (TF, TG or associated regulatory elements).
$ md5sum ibd2015_Primary_Natural_Killer_cells_from_peripheral_blood_snp2net.matd96cd9b32759f954cc37680dc6aeafd8 ibd2015_Primary_Natural_Killer_cells_from_peripheral_blood_snp2net.mat$ du -sh ibd2015_Primary_Natural_Killer_cells_from_peripheral_blood_snp2net.mat21M ibd2015_Primary_Natural_Killer_cells_from_peripheral_blood_snp2net.mat
In this example, there are 1081481 GWAS SNPs and 382443 of them are near the NK cell network (i.e. val=1).
This file contains the SNP-to-gene cis regulation scores derived from context-matching cis eQTL studies. This file corresponds to \((c_{jg}-1)\) in the RSS-NET model.
$ md5sum ibd2015_NK_snp2gene_cis.matdedad8e25773fad69576dbce0f7d9f93 ibd2015_NK_snp2gene_cis.mat$ du -sh ibd2015_NK_snp2gene_cis.mat165M ibd2015_NK_snp2gene_cis.mat
In this example, there are 10790012 SNP-gene pairs with cis regulation scores available, consiting of 829280 SNPs and 18230 genes. The cis regulation scores (val) range from 0 to 0.76.
$ pwd/Users/xiangzhu/GitHub/rss-net/examples/ibd2015_nkcell$ tree -f.├── ./analysis_template.m├── ./ibd2015_nkcell.m└── ./ibd2015_nkcell.sbatch0 directories, 3 files
3. Submit job arrays
For a given GWAS and a given regulatory network, all RSS-NET analysis tasks are almost identical and they only differs in hyper-parameter values. To exploit this, we run one RSS-NET analysis as a job array with multiple tasks that run in parallel.
To this end, we write a simple sbatch script ibd2015_nkcell.sbatch, and submit it to a cluster with Slurm available.
$ sbatch ibd2015_nkcell.sbatch
After the submission, multiple jobs should run in different nodes simultaneously.
$ squeue -u xiangzhu JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON) 62554249_107 owners ibd2015_ xiangzhu R 0:32 1 sh02-17n12 62554249_108 owners ibd2015_ xiangzhu R 0:32 1 sh02-17n12 62554249_109 owners ibd2015_ xiangzhu R 0:32 1 sh01-28n08 62554249_110 owners ibd2015_ xiangzhu R 0:32 1 sh01-17n18 62554249_111 owners ibd2015_ xiangzhu R 0:32 1 sh01-26n33 62554249_112 owners ibd2015_ xiangzhu R 0:32 1 sh01-27n30
Start at Mar 2, 2020, 3:04 PM.
End at Mar 2, 2020, 11:37 PM.
Job ID: 62554249
Array Job ID: 62554249_125
Cluster: sherlock
User/Group: xiangzhu/whwong
State: COMPLETED (exit code 0)
Nodes: 1
Cores per node: 8
CPU Utilized: 1-03:08:31
CPU Efficiency: 74.54% of 1-12:24:40 core-walltime
Job Wall-clock time: 04:33:05
Memory Utilized: 26.58 GB
Memory Efficiency: 85.06% of 31.25 GB