<jupyter_start><jupyter_text>Fine-mapping with PolyFun  Aim The purpose of this notebook ipmlements commands for [a functionally-informed fine-mapping workflow using the PolyFun method](https://github.com/omerwe/polyfun/wiki).  Methods Overview `PolyFun` offers the following features:1. Using and/or creating Functional Annotations2. Estimating Functional Enrichment Using `S-LDSC`3. Using and/or computating Prior Causal Probabilities from 14. Functionally Informed Fine Mapping with Finemapper5. Polygenic Localization with `PolyLoc`**Notice: this workflow does not implements 5 the `PolyLoc`.**  Input  1) GWAS summary statistics including the following variables:     - variant_id - variant ID     - P - p-value     - CHR - chromosome number     - BP - base pair position    - A1 - The effect allele (i.e., the sign of the effect size is with respect to A1)    - A2 - the second allele     - MAF - minor allele frequency     - BETA - effect size     - SE - effect size standard error2) Functional annotation files including the following columns:     - CHR - chromosome number    - BP base pair position (in hg19 coordinates)    - SNP - dbSNP reference number     - A1 - The effect allele     - A2 - the second allele    - Arbitrary additional columns representing annotations3) A `.l2.M` white-space delimited file containing a single line with the sums of the columns of each annotation4) LD-score files     - Strongly recommended that LD-score files include A1,A2 columns5) LD information, taken from one of three possible data sources:    - plink files with genotypes from a reference panel    - bgen file with genotypes from a reference panel    - pre-computed LD matrix    Optional if (4) is obtained and no plans to compute prior causal probabilities non-parametrically 6) Ld-score weights files.    - Strongly recommended that weight files include A1,A2 columns  OutputA `.gz` file containing input summary statistics columns and additionally the following columns:- PIP - posterior causal probability- BETA_MEAN - posterior mean of causal effect size (in standardized genotype scale)- BETA_SD - posterior standard deviation of causal effect size (in standardized genotype scale)- CREDIBLE_SET - the index of the first (typically smallest) credible set that the SNP belongs to (0 means none).  Workflow Step 1 and 2 are optional if using pre-computed prior causal probabilities  Step 1: Obtain functional annotations  For each chromosome, the following files need to be obtained: 1) A `.gz` or `.parquet` annotations file containing the following columns:- CHR - chromosome number- BP base pair position- SNP - dbSNP reference number - A1 - The effect allele - A2 - the second allele- Arbitrary additional columns representing annotations 2) A `.l2.M` white-space delimited file containing a single line with the sums of the columns of each annotation3) (Optional) A `l2.M_5_50` file that is the `.l2.M` file but only containing common SNPS (MAF between 5% and 50%)  The above files can be obtained either by using existing function annotation files, or by creating your own through other software such as `TORUS`.Existing function annotation files example: functional annotations for ~19 million UK Biobank imputed SNPs with MAF>0.1%, based on the baseline-LF 2.2.UKB annotations.Download (30G): https://data.broadinstitute.org/alkesgroup/LDSCORE/baselineLF_v2.2.UKB.polyfun.tar.gz  Step 2: Compute LD-scores for annotations  Precomputed LD-score files can be used. LD-score files can also be generated through the methods below:  Method 1: Compute with reference panel of sequenced individuals  Reference panel should have at least 3000 sequenced individuals from target population.<jupyter_code>[global]
parameter: container = "/mnt/mfs/statgen/containers/xqtl_pipeline_sif/polyfun.sif"
parameter: wd = path("./")
parameter: exe_dir = "/usr/local/bin/"
parameter: name = "demo"
parameter: genoFile = path("./")
parameter: annot_file = path("./")
parameter: sumstats = path("./")
[ld_score]
input: annot_file, genoFile
output: f'{wd:a}/{name}.ref.ldscore.parquet'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/compute_ldscores.py \
        --bfile $[_input[1]:n] \
        --annot $[_input[0]] \
        --out $[_output]<jupyter_output><empty_output><jupyter_text>Method 2: Compute with pre-computed UK Biobank LD matrices  Matrices download: https://data.broadinstitute.org/alkesgroup/UKBB_LD<jupyter_code>[ld_score_ukb]
input: annot_file
output: f'{wd:a}/{name}.ukb.ldscore.parquet'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/compute_ldscores_from_ld.py \
        --annot $[_input[0]] \
        --ukb \
        --out $[_output]<jupyter_output><empty_output><jupyter_text>Method 3: Compute with own pre-computed LD matrices Own pre-computed LD matrices should be in `.bcor` format.<jupyter_code>[ld_score_own]
parameter: sample_size = int
parameter: bcor_files = paths
input: annot_file,bcor_files
output: f'{wd:a}/{name}.original.ldscore.parquet'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/compute_ldscores_from_ld.py $[_input[1]] \
        --annot $[_input[0]] \
        --out $[_output] \
        --n $[sample_size] \<jupyter_output><empty_output><jupyter_text>Step 3: Compute Prior Causal Probabilities  Method 1: Use precomputed prior causal probabilities Use precomputed prior causal probabilities of 19 million imputed UK Biobank SNPs with MAF>0.1%, based on a meta-analysis of 15 UK Biobank traits.<jupyter_code>[prior_causal_prob]
input: sumstats
output: f'{wd:a}/{name}.pcp.gz'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container
    python $[exe_dir]/extract_snpvar.py \
        --sumstats $[_input] \
        --out $[_output] \
        --allow-missing<jupyter_output><empty_output><jupyter_text>Method 2: Compute via L2-regularized extension of S-LDSC (preferred) Compute via an L2-regularized extension of stratified LD-score regression (S-LDSC). Use the annotation and LD-score files produced in Step1.  1) Create a munged summary statistics file in a PolyFun-friendly parquet format.<jupyter_code>[munged_sumstats]
parameter: sample_size = 472868 
parameter: min_info = 0.6
parameter: min_maf = 0.01
input: sumstats
output: f'{wd:a}/{name}.sumstats_munged.parquet'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '127G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/munge_polyfun_sumstats.py \
      --sumstats $[_input] \
      --n $[sample_size] \
      --out $[_output] \
      --min-info $[min_info] \
      --min-maf $[min_maf]<jupyter_output><empty_output><jupyter_text>2) Run PolyFun with L2-regularized S-LDSC- Require at least 45 GB of mem<jupyter_code>[L2_SLDSC]
# a ld score file with surfix l2.ldscore.parquet
parameter: ref_ld = path
# another ld score file with surfix l2.ldscore.parquet, different from ref_ld
parameter: ref_wgt = ref_ld
parameter: partitions = ""
input: ref_ld, ref_wgt,output_from("munged_sumstats")
# parameter: sumstat = _input[2]
output: f'{wd:a}/{name}.ldsldsc.parquent'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '127G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/polyfun.py \
        --compute-h2-L2 \
        --output-prefix $[_output] \
        --sumstats $[_input[2]] \
        --ref-ld-chr $[_input[0]:nnnn].\
        --w-ld-chr $[_input[1]:nnnn]. \
        --allow-missing $["" if partitions else "--no-partitions"]<jupyter_output><empty_output><jupyter_text>Method 3: Compute Non-parametrically 1) Create a munged summary statistics file in a PolyFun-friendly parquet format.Duplicated cells are commented out, the input of [ld_snpbin] is the output from [L2_regu_SLDSC]<jupyter_code>#[munged_sumstats2]
#parameter: sumstats = AD_sumstats_Jansenetal_2019sept.txt.gz
#parameter: sample_size = int
#parameter: container = none
#bash: container = container 
#    mkdir -p SLDSC_output
#    python munge_polyfun_sumstats.py \
#      --sumstats sumstats \
#      --n sample_size \
#      --out /SLDSC_output/sumstats_munged.parquet \
#      --min-info 0 \
#      --min-maf 0<jupyter_output><empty_output><jupyter_text>2) Run PolyFun with L2-regularized S-LDSC<jupyter_code># [L2_regu_SLDSC2]
# 
# parameter: container = none
# paramter: ref_ld = example_data/annotations.
# parameter: ref_wgt = example_data/weights.
# bash: container=container
#     python polyfun.py \
#     --compute-h2-L2 \
#     --output-prefix output/testrun \
#     --sumstats example_data/sumstats.parquet \
#     --ref-ld-chr ref_ld \
#     --w-ld-chr ref_wgt<jupyter_output><empty_output><jupyter_text>3) Compute LD-scores for each SNP bin<jupyter_code>[ld_snpbin]
depends: sos_step("L2_regu_SLDSC")
parameter: chrom = int
input: annot_file, genoFile
output: f'{wd:a}/{name}.snpbin.ldscore.parquet'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
     python $[exe_dir]/polyfun.py \
        --compute-ldscores \
        --bfile-chr $[_input[1]:n] \
        --output-prefix $[_output] \
        --chr $[chrom]<jupyter_output><empty_output><jupyter_text>4) Re-estimate per-SNP heritabilities via S-LDSC<jupyter_code>#[re_SLDSC]
#bash:
#    python polyfun.py \
#    --compute-h2-bins \
#    --output-prefix output/testrun \
#    --sumstats example_data/sumstats.parquet \
#    --w-ld-chr example_data/weights.

[L2_SLDSC_bins]
paramter: ref_ld = path
parameter: ref_wgt = ref_ld
parameter: partitions = ""
input: ref_ld, ref_wgt,output_from("munged_sumstats")
output: f'{wd:a}/{name}.txt.gz'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/polyfun.py \
        --compute-h2-bins \
        --output-prefix $[_output] \
        --sumstats $[_input[2]] \
        --ref-ld-chr $[_input[0]]\
        --w-ld-chr $[_input[1]]<jupyter_output><empty_output><jupyter_text>Step 4: Functionally informed fine mapping with finemapper Input summary statistics file must have `SNPVAR` column (per-SNP heritability) to perform functionally-informed fine-mapping. To fine-map without annotations, use additional parameter `--non-funct`. The summary statistical file then will not require the `SNPVAR` column.<jupyter_code>[fine_mapping]
parameter: sample_size = 383290
parameter: chrom = 1
parameter: start = 46000001
parameter: end = 49000001
#parameter: output_path = "output/finemap.1.46000001.49000001.gz"
parameter: max_num_causal = 5
input: genoFile,sumstats
output: f'{wd:a}/output/finemap.{chrom}.{start}.{end}.gz'
task: trunk_workers = 1, trunk_size = 1, walltime = '24h',  mem = '30G', tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout' , container = container
    python $[exe_dir]/finemapper.py \
        --geno $[_input[0]:n] \
        --sumstats $[_input[1]]  \
        --n $[sample_size] \
        --chr $[chrom] \
        --start $[start] \
        --end $[end] \
        --method susie \
        --max-num-causal $[max_num_causal] \
        --cache-dir $[_output:d]/cache \
        --out $[_output]<jupyter_output><empty_output><jupyter_text>Minimal Working Example  Example 1: Functionally-informed fine-mapping using summary statistics file with precomputed prior causal probabilities<jupyter_code>nohup sos run ~/GIT/xqtl-pipeline/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb prior_causal_prob \
    --sumstats /home/at3535/polyfun/AD_sumstats_Jansenetal_2019sept.txt.gz  \
    -J 200 -q csg \
    -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &
nohup sos run ~/GIT/xqtl-pipeline/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb fine_mapping \
    --sumstats /home/at3535/polyfun/AD_sumstats_Jansenetal_2019sept.txt.gz  \
    --genoFile /mnt/mfs/statgen/ROSMAP_xqtl/dataset/snvCombinedPlink/chr1.bed \
    -J 200 -q csg \
    -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &<jupyter_output><empty_output><jupyter_text>Example 2: Functionally-informed fine-mapping using summary statistics file generated from pre-obtained annotation and LD-score files<jupyter_code>nohup sos run ~/GIT/xqtl-pipeline/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb munged_sumstats  \
--sumstats /home/at3535/polyfun/GCST90012877_buildGRCh37_colrenamed.txt.gz \
-J 200 -q csg -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &
nohup sos run ~/GIT/xqtl-pipeline/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb L2_SLDSC  \
--sumstats demo.sumstats_munged.parquet \
--ref_ld /mnt/mfs/statgen/tl3030/baselineLF2.2.UKB/baselineLF2.2.UKB.1.l2.ldscore.parquet \
-J 200 -q csg -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &
nohup sos run ~/GIT/xqtl-pipeline/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb fine_mapping \
    --sumstats demo.ldsldsc.parquent  \
    --genoFile /mnt/mfs/statgen/ROSMAP_xqtl/dataset/snvCombinedPlink/chr1.bed \
    -J 200 -q csg -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &<jupyter_output><empty_output><jupyter_text>Summary  Example 1:<jupyter_code>bash:
    gzcat output/finemap.1.46000001.49000001.gz | head
import numpy as np
import pandas as pd

data = pd.read_csv('output/finemap.1.46000001.49000001.gz', sep="\t")

data.head(5)
    
num_var_cs = np.count_nonzero(data['CREDIBLE_SET'])
total_cs = len(data.CREDIBLE_SET.unique())- 1
avg_var_cs = float(num_var_cs) / total_cs
pip50 = sum(1 for i in data['PIP'] if i >0.5)
pip95 = sum(1 for i in data['PIP'] if i >0.95)

result = "Number of variants with PIP > 0.5: " + str(pip50) + "\n" + "Number of variants with PIP > 0.95: " + str(pip95) + "\n" \
    + "Number of variants that have credible sets: " + str(num_var_cs) + "\n" \
    + "Number of unique credible sets: " + str(total_cs) + "\n" \
    + "Average number of variants per credible set: " + str(avg_var_cs) 


with open('results.txt', 'a') as the_file:
    the_file.write(result)

with open('results.txt') as f:
    contents = f.readlines()
    print(contents)<jupyter_output>['Number of variants with PIP > 0.5: 3\n', 'Number of variants with PIP > 0.95: 3\n', 'Number of variants that have credible sets: 3\n', 'Number of unique credible sets: 3\n', 'Average number of variants per credible set: 1.0']<jupyter_text>Example 2:<jupyter_code>bash:
    gzcat output/finemap.1.460000010.49000000.gz | head
import numpy as np
import pandas as pd

data = pd.read_csv('output/finemap.1.46000000.49000000.gz', sep="\t")

data.head(5)
    
num_var_cs = np.count_nonzero(data['CREDIBLE_SET'])
total_cs = len(data.CREDIBLE_SET.unique())- 1
avg_var_cs = float(num_var_cs) / total_cs
pip50 = sum(1 for i in data['PIP'] if i >0.5)
pip95 = sum(1 for i in data['PIP'] if i >0.95)

result = "Number of variants with PIP > 0.5: " + str(pip50) + "\n" + "Number of variants with PIP > 0.95: " + str(pip95) + "\n" \
    + "Number of variants that have credible sets: " + str(num_var_cs) + "\n" \
    + "Number of unique credible sets: " + str(total_cs) + "\n" \
    + "Average number of variants per credible set: " + str(avg_var_cs) 


with open('results.txt', 'a') as the_file:
    the_file.write(result)

with open('results.txt') as f:
    contents = f.readlines()
    print(contents)<jupyter_output>['Number of variants with PIP > 0.5: 3\n', 'Number of variants with PIP > 0.95: 3\n', 'Number of variants that have credible sets: 3\n', 'Number of unique credible sets: 3\n', 'Average number of variants per credible set: 1.0']
{'fraction_non_alphanumeric': 0.12934960830752415, 'fraction_numerical': 0.027691747130624885, 'mean_word_length': 4.005471124620061, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 73, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '31410', 'license': 'MIT', 'n_tokens_mistral': 6137, 'n_tokens_neox': 5605, 'n_words': 1512, 'path': 'pipeline/polyfun.ipynb', 'repo_name': 'seriousbamboo/xqtl-pipeline', 'size': '16467'}