<jupyter_start><jupyter_text>Enrichment analysis: TORUS and GREGORThis workflow document has several pipelines written in SoS Workflow language to perform data preparation and analysis for enrichment of GWAS signals in given annotation features. For enrichment analysis currently it implements the `torus` pipeline, based on some `snakemake` codes originally written by Jean at Xin He Lab, UChicago.<jupyter_code>%revisions -s -n 10
sos run gwas_enrichment.ipynb -h<jupyter_output>usage: sos run gwas_enrichment.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  zscore2bed
  get_variants
  merge_annotations
  range2var_annotation
  enrichment
  deepsea_apply
  gregor

Global Workflow Options:
  --cwd . (as path)
                        working directory
  --hg19 . (as path)
                        hg19 path
  --deepsea-model . (as path)
                        Deepsea model for use with `deepsee_apply`
  --annotation-dir . (as path)
                        Path to directory of annotation files
  --z-score . (as path)
                        Path to z-score file
  --single-annot . (as path)
                        Path to list of single annotati[...]<jupyter_text>Copy this GIT depository```git clone https://github.com/cumc/bioworkflows``` To run enrichment analysis (overview) Step 1 (start from `bed` format annotation)```sos run fine-mapping/gwas_enrichment.ipynb range2var_annotation --z-score ... --single-annot ... --multi-annot ...``` Step 2 (start from binary annotation for SNPs)```sos run fine-mapping/gwas_enrichment.ipynb enrichment --z-score ... --single-annot ... --multi-annot ...```  Software intallation `SoS`    SoS introduction [webpage](https://vatlab.github.io/sos-docs/)Installation:```pip install sos sos-notebook```SoS basic usage [webpage](https://vatlab.github.io/sos-docs/running.htmlcontent) `docker`Some steps in this pipeline uses docker images we prepare. If you do not have docker installed, you can install it via:- Run commands below:```curl -fsSL get.docker.com -o get-docker.shsudo sh get-docker.shsudo usermod -aG docker $USER```- Log out and log back in (no need to reboot computer).**In this pipeline, `bedops` and `torus` is executed via docker image we prepared and uploaded to dockerhub.com. The pipeline should automatically download the image and run the docker container instance.**  Reference data preparation hg19.fa- only for `deepsea` and `gregor` steps; if you only use enrichment analysis, you do not need to download it.- Genome Reference: hg19 (GRCh37) or hg38 (GRCh38). - Download [link](http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz)  Annotation files Annotation files downloadAlkes Price's lab has some annotation files in `bed` format available for download:https://data.broadinstitute.org/alkesgroup/LDSCORE/eg `baselineLD_v2.2_bedfiles.tgz`. `bed` format annotation input**If you use the type of annotation file below, you need to run step 1 and 2.**Annotation files are in `bed` format (for example: Promoter_UCSC.bed).        chr1	9873	16361        chr1	32610	36610        chr1	67090	71090        ... Binary annotations of interest**If you use the type of annotation file below, you only need to run step 2.**0/1 binary indicated. For example, a few lines in `Coding_UCSC_annotation.torus.gz`.    SNP     Coding_UCSC_d    9:140192349:G:C    0    9:140192576:C:T    0    9:140194020:G:A    0    9:140194735:G:T    1    9:140194793:A:G    1    9:140195019:C:A    1There are many annotations one can use. For this workflow you should prepare two lists:- For the ones that you'd like to use for one variable `torus`, put their file name prefix in `data/single_annotations.txt`- For the ones that you'd like to use for multi-variable `torus`, put their file name prefix in `data/multi_annotations.txt`If you want to remove one annotation just use `` to comment it out in the list files.<jupyter_code>%preview ../data/general_annotations.txt -l 10<jupyter_output><empty_output><jupyter_text>`torus` workflow data format GWAS z-scores specificationExample GWAS data (European)- [SCZ Sweden data](https://www.med.unc.edu/pgc/results-and-downloads/data-use-agreement-forms/SwedenSCZ_data_download_agreement)- [PTSD European Ancestry](https://www.med.unc.edu/pgc/results-and-downloads/data-use-agreement-forms/PTSD%20EA_data_download_agreement)Format `chr:pos:alt:ref   ld_label   z-score`. 3 columns in total: `chr:pos:alt:ref` (chromosome, position, alternative allele and reference allele), `ld_label` (LD chunk label) and `z-score` (summary statistics).```1:10583:A:G        1      0.116319...22:51228910:A:G    1703   -0.866894```where the 2nd column is the analysis block ID (see below). All GWAS summary statistics have to be converted to this format in order to use. LD block filesThe last column indicates LD chunk. LD chunk blocks depend on chromosome and position. There are 1703 LD chunks in European human genome.```            chr1    10583       1892607     1            chr1    1892607     3582736     2            ...            chr22   49824534    51243298    1703```- LD block reference paper [link](https://academic.oup.com/bioinformatics/article/32/2/283/1743626). SNP mapGWAS SNP map. Identical number of SNPs with GWAS.Format `chrom.pos  chr  pos`. 3 columns in total: `chrom.pos` (example: `chr1.729679`), `chr` and `pos`.```chr1.729679	1	729679chr1.731718	1	731718chr1.734349	1	734349...```  Minimal working example for `torus` workflow- GWAS z-score example [link](https://www.dropbox.com/s/y2q29a71bb2wmr1/SCZGWAS.zscore.gz?dl=0)- SNP map example [link](https://www.dropbox.com/s/oi4yr0au3oy3hi2/smap.gz?dl=0)- annotation example: `Promoter_UCSC.bed` download [link](https://www.dropbox.com/s/huo5k78wptxmrpe/Promoter_UCSC.bed?dl=0). 22,436 lines in total. - analysis block example: [link](https://www.dropbox.com/s/up78cnpkpoodark/ld_chunk.bed?dl=0). Used to identify LD chunk for each GWAS SNP.  Workflow**You can edit and change the following 5 directories and file pathes.**First edit and run the following 5 command to define bash variables.```work_dir=~/Documents/GWAS_ATACanno_dir=~/Documents/GWAS_ATAC/bed_annotationsz_file=~/Documents/GWAS_ATAC/SCZGWAS_zscore/SCZGWAS.zscore.gzsingle=~/GIT/fine-mapping/data/general_annotations.txtblk=~/Documents/GWAS_ATAC/blocks.txtsnps=~/Documents/GWAS_ATAC/SCZGWAS_zscore/SCZGWAS.snps.gz```Then run following commands.```sos run fine-mapping/gwas_enrichment.ipynb range2var_annotation --cwd $work_dir --annotation_dir $anno_dir --z-score $z_file --single-annot $singlesos run fine-mapping/gwas_enrichment.ipynb enrichment --cwd $work_dir --annotation_dir $anno_dir --z-score $z_file --single-annot $single --blocks $blk --snps $snps```<jupyter_code>[global]
# working directory
parameter: cwd = path()
# hg19 path
parameter: hg19 = path()
# Deepsea model for use with `deepsee_apply`
parameter: deepsea_model = path()
# Path to directory of annotation files
parameter: annotation_dir = path() # 2
# Path to z-score file
parameter: z_score = path() # 3
# Path to list of single annotations to use
parameter: single_annot = path() # 4
# Path to lists of multi-annotations to use
parameter: multi_annot = paths() # 5

def get_bed_name(directory, basename):
    filename = path(f"{directory}/{basename}.bed")
    if not filename.is_file():
        raise ValueError(f"Cannot find file ``{directory}/{basename}.bed``")
    return filename

try:
    single_anno = [get_bed_name(f"{annotation_dir:a}", x.split()[0]) for x in open(single_annot).readlines() if not x.startswith('#')]
    multi_anno = dict([(f'{y:bn}', [get_bed_name(f"{annotation_dir:a}", x.split()[0]) for x in open(y).readlines() if not x.startswith('#')]) for y in multi_annot])
except (FileNotFoundError, IsADirectoryError):
    single_anno = []
    multi_anno = dict()
out_dir = f'{cwd:a}/{z_score:bn}'.replace('.', '_')<jupyter_output><empty_output><jupyter_text>Utility steps  Convert variants from z-score file to bed format<jupyter_code># Auxiliary step to get variant in bed format based on variant ID in z-score file
[zscore2bed_1]
parameter: in_file = path()
input: in_file
output: f'{_input:n}.bed.unsorted'
R: expand = "${ }", container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout'
    library(readr)
    library(stringr)
    library(dplyr)
    var_file <- ${_input:r}
    out_file <- ${_output:r}

    variants <- read_tsv(var_file, col_names=FALSE, col_types='ccd')
    colnames(variants) = c('variant', 'ld', 'zscore')
    var_info <- str_split(variants$variant, ":")
    variants <- mutate(variants, chr = paste0("chr", sapply(var_info, function(x){x[1]})), 
                                 pos = sapply(var_info, function(x){x[2]})) %>%
                mutate(start = as.numeric(pos), stop=as.numeric(pos) + 1) %>%
                select(chr, start, stop, variant)
    options(scipen=1000) # So that positions are always fully written out)
    write.table(variants, file=out_file, quote=FALSE, col.names=FALSE, row.names=FALSE, sep="\t")

[zscore2bed_2]
output: f'{_input:n}'
bash: expand = True, container = 'gaow/atac-gwas', workdir = cwd
     sort-bed {_input} > {_output}

[get_variants: provides = '{data}.bed']
output: f'{data}.bed'
sos_run('zscore2bed', in_file = f'{_output:n}.gz')<jupyter_output><empty_output><jupyter_text>Merge annotations for joint torus analysis<jupyter_code># Auxiliary step to merge multiple annotations
[merge_annotations]
parameter: out_file=path()
parameter: data_files=paths()
input: data_files
output: out_file
R: expand = '${ }', container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout'

    library(readr)
    library(dplyr)
    library(stringr)

    out_name <- ${_output:r}
    annots <- c(${_input:r,})
    var_annot <- read_tsv(annots[1], col_names=T)
    for(i in 2:length(annots)){
        var_annot2 <- read_tsv(annots[i], col_names=T)
        stopifnot(all(var_annot$SNP == var_annot2$SNP))
        var_annot <- cbind(var_annot, var_annot2[,2])
    }
    write.table(var_annot, file=gzfile(out_name),
                row.names=FALSE, quote=FALSE, sep="\t")<jupyter_output><empty_output><jupyter_text>Prepare `torus` format binary annotations<jupyter_code># Get variants in data that falls in target region
[range2var_annotation_1]
depends: f'{z_score:n}.bed'
input: set(paths(single_anno + list(multi_anno.values()))), group_by = 1
output: f'{out_dir}/{_input:bn}.{z_score:bn}.bed'
bash: expand = True, container = 'gaow/atac-gwas', workdir = cwd, volumes = [f'{_input:d}:{_input:d}', f'{_output:d}:{_output:d}']
    bedops -e {z_score:n}.bed {_input} > {_output}
# Make binary annotation file
[range2var_annotation_2]
parameter: discrete = 1
depends: z_score
output: f'{_input:n}.gz'
R: expand = "${ }", container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout', volumes = [f'{_input:d}:{_input:d}', f'{_output:d}:{_output:d}']
    library(readr)
    library(dplyr)
    library(stringr)

    variant_tsv <- ${z_score:r}
    annotation_var_bed <- ${_input:r}
    annot_name <- ${_input:bnr} %>% str_replace(paste0(".",${z_score:bnr}), "")
    out_name <- ${_output:r}

    vars <- read_tsv(variant_tsv, col_names=FALSE, col_types='ccd')[,1]
    annot_vars = read_tsv(annotation_var_bed, col_names=FALSE, col_types='cddc')
    names(vars) <- "SNP"
    vars <- vars %>%
            mutate(annot_d = case_when(SNP %in% annot_vars$X4 ~ 1,
                                                        TRUE ~ 0))
    names(vars)[2] <- paste0(annot_name, '${"_d" if discrete else "_c"}')
    write.table(vars, file=gzfile(out_name),
                col.names=TRUE, row.names=FALSE, sep="\t", quote=FALSE)
# Obtain multiple annotations per SNP for enrichment analysis
[range2var_annotation_3]
for k, value in multi_anno.items():
    sos_run('merge_annotations', out_file = f'{out_dir}/{k}.{z_score:bn}.gz', data_files = [f'{out_dir}/{v:bn}.{z_score:bn}.gz' for v in paths(value)])<jupyter_output><empty_output><jupyter_text>Enrichment analysis via `torus`<jupyter_code># Run torus with annotation + z-score file
[enrichment_1 (run torus)]
depends: z_score
parameter: blocks = path()
parameter: snps = path()
input_files = [f'{out_dir}/{value:bn}.{z_score:bn}.gz' for value in paths(single_anno)] + [f'{out_dir}/{k}.{z_score:bn}.gz' for k in multi_anno]
fail_if(len(input_files) == 0, msg = "No annotations to use! Please use ``--single-annot`` and ``--multi-annot`` to pass annotation lists")
input: input_files, group_by = 1
output: f'{_input:n}.torus'

bash: container = 'gaow/atac-gwas', workdir = cwd, expand = True, stdout = f'{_output:n}.stdout', stderr = f'{_output:n}.stderr', volumes = [f'{_input:d}:{_input:d}', f'{_output:d}:{_output:d}']
    rm -rf {_output:n}_prior 
    torus -d {z_score} -smap {snps} -gmap {blocks} -annot {_input} -est --load_zval -dump_prior {_output:n}_prior > {_output}
# Consolidate all torus result into one table
[enrichment_2 (make output table)]
input: group_by = 'all'
output: f'{out_dir}/{z_score:bn}.torus.merged.csv'
python: expand = '${ }', workdir = cwd
    # get log2 folds of enrichment
    import pandas as pd
    def get_torus_res(fn):
        import math, os
        res = pd.read_table(fn, sep = "\s+", header = None, names = ["annotation", "log2_fold_enrichment", "CI_low", "CI_high"])
        res = res.iloc[1:]
        res["source"] = [os.path.split(fn)[-1] for x in range(res.shape[0])]
        res["annotation"] = res["annotation"].apply(lambda x: x.rsplit('.',1)[0])
        res["log2_fold_enrichment"] = [round(math.log2(math.exp(x)),4) for x in res["log2_fold_enrichment"]]
        res["CI_low"] = [round(math.log2(math.exp(x)),4) for x in res["CI_low"]]
        res["CI_high"] = [round(math.log2(math.exp(x)),4) for x in res["CI_high"]]
        return res
    res = pd.concat([get_torus_res(x) for x in [${_input:r,}]])
    res.to_csv(${_output:r}, index=False)<jupyter_output><empty_output><jupyter_text>Apply pre-trained `deepsea` model to variantsCredits to Yanyu Liang. Required inputs are:- Path to hg19 reference genome- Path to model HDF5 file- Path to list of variants, in the format of:```chr1	68090	68091	G	Tchr1	68090	68091	G	Achr1	68090	68091	G	Cchr1	68091	68092	C	Gchr1	68091	68092	C	Tchr1	68091	68092	C	Achr1	68092	68093	T	Achr1	68092	68093	T	Cchr1	68092	68093	T	Gchr1	68093	68094	A	C``` To run this workflow:    ```sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants /path/to/variant/list.file``````sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/scz.swe.var.txt```It requires much memory. So we test for first 10K line of variant and it worked.```sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/test.var.txtsos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/test100K.var.txt```Break 9898079 snps into two halves and run deepsea separately.```less scz.swe.var.txt | head -4949039 | gzip > half1.var.txt``````sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/half1.var.txtsos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/half2.var.txt```<jupyter_code>[deepsea_apply_1 (prepare config file)]
parameter: variants = path()
input: variants
output: f'{_input:n}.config'
report: expand = True, output = _output
    var_list: {_input:r}
    pred_model:
      path: '{deepsea_model}'
      label:
        CN: 1
        DN: 2
        GA: 3
        ips: 4
        NSC: 5
    reference_genome: '{hg19}'
    script_dir: '/opt/deepann'
    out_dir: {_output:dr}
    out_prefix: '{_output:bn}_deepsea'

    ##### some default setup #####
    #### usually don't change ####
    check_allele: False
    design: '499-1-500'
[deepsea_apply_2 (apply deepsea weights)]
output: f'{_input:n}.deepsea.list'
bash: expand = True, container = 'gaow/deepann', volumes = [f'{hg19:d}:{hg19:d}', f'{deepsea_model:d}:{deepsea_model:d}', f'{os.path.expanduser("~")}:/home/$USER'], extra_args = "-e HOME=/home/$USER -e USER=$USER"
    snakemake --snakefile /opt/deepann/Snakefile --configfile {_input} && ls {_input:n}_deepsea/output__*.gz > {_output}<jupyter_output><empty_output><jupyter_text>Enrichment analysis via GREGOR To properly perform enrichment analysis we want to match the control SNPs with the SNPs of interest -- that is, SNPs inside CS -- in terms of LD, distance to nearest gene and MAF. The [GREGOR](http://csg.sph.umich.edu/GREGOR/index.php/site/download) software can generate list of matched SNPs. I will use SNPs inside CS as input and expect a list of output SNPs matching these inputs.GREGOR is release under University of Michigan license so I'll not make it into a docker image. So path to GREGOR directory is required. Also we need reference files, prepared by:```cat \    GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz.part.00 \    GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz.part.01 \    > GREGOR.AFR.ref.r2.greater.than.0.7.tar.gztar zxvf GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz```MD5SUM check:```AFR.part.0	( MD5: 9926904128dd58d6bf1ad4f1e90638af )AFR.part.1	( MD5: c1d30aff89a584bfa8c1fa1bdc197f21 )```Same for EUR data-set.  To run GREGOR`PGC3````sos run workflow/gwas_enrichment.ipynb gregor --gregor-input ~/Documents/GWAS_ATAC/scz2_zscore/gregor.scz2.txt.gz --single-annot data/gregor.txt````PGC2````sos run workflow/gwas_enrichment.ipynb gregor --gregor-input ~/Documents/GWAS_ATAC/pgc2_zscore/gregor.pgc2.txt.gz --single-annot data/gregor.txt```<jupyter_code>[gregor_1 (make SNP index)]
depends: executable('zcat')
parameter: gregor_input = path()
input: gregor_input
output: f'{_input:nn}.rsid.txt', f'{_input:nn}.annotations.list'
bash: expand = '${ }'
    zcat ${_input} | cut -f 2,3 -d "_" | sed 's/_/:/g' > ${_output[0]}

with open(_output[1], 'w') as f:
    f.write('\n'.join(single_anno))

[gregor_2 (make configuration file)]
depends: executable('sed')
parameter: gregor_db = path('~/Documents/hg19/GREGOR_DB')
parameter: pop = 'EUR'
output: f'{_input[0]:nn}.gregor.conf'
report: output = f'{_output}', expand = True
    ##############################################################################
    # CHIPSEQ ENRICHMENT CONFIGURATION FILE
    # This configuration file contains run-time configuration of
    # CHIP_SEQ ENRICHMENT
    ###############################################################################
    ## KEY ELEMENTS TO CONFIGURE : NEED TO MODIFY
    ###############################################################################
    INDEX_SNP_FILE = {_input[0]}
    BED_FILE_INDEX = {_input[1]} 
    REF_DIR = {gregor_db}
    R2THRESHOLD = 0.7 ## must be greater than 0.7
    LDWINDOWSIZE = 10000 ## must be less than 1MB; these two values define LD buddies
    OUT_DIR = {_output:nn}_gregor_output
    MIN_NEIGHBOR_NUM = 10 ## define the size of neighborhood
    BEDFILE_IS_SORTED = true  ## false, if the bed files are not sorted
    POPULATION = {pop}  ## define the population, you can specify EUR, AFR, AMR or ASN
    TOPNBEDFILES = 2 
    JOBNUMBER = 10
    ###############################################################################
    #BATCHTYPE = mosix ##  submit jobs on MOSIX
    #BATCHOPTS = -E/tmp -i -m2000 -j10,11,12,13,14,15,16,17,18,19,120,122,123,124,125 sh -c
    ###############################################################################
    #BATCHTYPE = slurm   ##  submit jobs on SLURM
    #BATCHOPTS = --partition=broadwl --account=pi-mstephens --time=0:30:0
    ###############################################################################
    BATCHTYPE = local ##  run jobs on local machine

bash: expand = True
    sed -i '/^$/d' {_output}<jupyter_output><empty_output><jupyter_text>GREGOR is written in `perl`. Some libraries are required before one can run GREGOR:```sudo apt-get install libdbi-perl libswitch-perl libdbd-sqlite3-perl```<jupyter_code>[gregor_3 (run gregor)]
depends: executable('perl')
parameter: gregor_path = path('~/Documents/GREGOR')
output: f'{_input:nn}_gregor_output/StatisticSummaryFile.txt'
bash: expand = True
    perl {gregor_path}/script/GREGOR.pl --conf {_input} && touch {_output}

[gregor_4 (format output)]
depends: executable('sed')
output: f'{_input:n}.csv'
bash: expand = True
    sed 's/\t/,/g' {_input} > {_output}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1539391594525318, 'fraction_numerical': 0.03259660492334478, 'mean_word_length': 4.253556910569106, 'pattern_counts': {'":': 1, '<': 58, '<?xml version=': 0, '>': 52, 'https://': 11, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '40784', 'license': 'MIT', 'n_tokens_mistral': 7753, 'n_tokens_neox': 7060, 'n_words': 1826, 'path': 'code/integrative_analysis/enrichment/gwas_enrichment.ipynb', 'repo_name': 'seriousbamboo/xqtl-pipeline', 'size': '20677'}