4_make_annot_ldscore.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os\n",
    "import glob\n",
    "import time\n",
    "import subprocess as sp\n",
    "import multiprocessing as mp\n",
    "import psutil\n",
    "\n",
    "from path_configure import *\n",
    "\n",
    "debug=False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_annot_chr_cm(pheno_code,chrN,scale=1,mode='cm'):\n",
    "    dest=annot_path.format(pheno_code,chrN)\n",
    "    if os.path.exists(dest):\n",
    "        print(\"Congratulations!. annot of\",pheno_code,chrN,\"exists. passed.\")\n",
    "        return\n",
    "    tosave=pd.read_csv(bim_path.format(chrN),sep='\\t',names=['CHR','SNP','CM','BP','A1','A2'])\n",
    "\n",
    "    if mode=='cm':\n",
    "        cm_series=tosave['CM']\n",
    "    elif mode=='bp':\n",
    "        cm_series=tosave['BP']/1000000\n",
    "    cm_max=int(np.ceil(np.max(cm_series)/scale)*scale)\n",
    "    cm_min=int(np.floor(np.min(cm_series)/scale)*scale)\n",
    "    print(pheno_code,\"chr:\",chrN,\"scale:\",scale,\"{} min,max\".format(mode),cm_min,cm_max)\n",
    "    category=np.zeros(len(cm_series),dtype=np.int32)\n",
    "    category_name=[]\n",
    "    category_count=0\n",
    "    for i in range(cm_min,cm_max,scale):\n",
    "        cm_index=np.logical_and(cm_series>=(i),cm_series<(i+scale))\n",
    "        print(pheno_code,\"chr:\",chrN,\"SNPs in {} range\".format(mode),i,i+scale,\"are\",sum(cm_index),end='')  \n",
    "        if sum(cm_index)>100:\n",
    "            category_count=category_count+1\n",
    "            category[cm_index]=category_count\n",
    "            category_name.append(str(chrN)+'.'+str(i))\n",
    "            print('-> included')\n",
    "        else:\n",
    "            print('-> ignored')\n",
    "            #print(pheno_code,\"chr:\",chrN,\"SNPs in {} range\".format(mode),i,i+scale,\"are ignored\")  \n",
    "\n",
    "    category_onehot=pd.DataFrame(np.vstack([np.zeros(category_count),np.eye(category_count)])[category],columns=category_name,dtype='int32')\n",
    "\n",
    "    tosave.reset_index(drop=False,inplace=True)\n",
    "    tosave=pd.concat([tosave,category_onehot],axis=1, join='inner')\n",
    "    print(tosave[tosave['22.32']==1])\n",
    "    #tosave.to_csv(dest,sep='\\t',columns=['CHR','BP','SNP','CM']+category_name,index=False)\n",
    "    print(\"Saved\",dest)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Process level parallelism for shell commands\n",
    "def estimate_ldscore_chrN(annot_name,chrN,mode='annot'):\n",
    "    \"\"\"Defines the work unit on an input file\"\"\"\n",
    "    #ldsc.py --l2 --bfile /scratch/ch6845/ukb_ld/ukb_imp_plink_gwas/ukb_imp_chr1_v3 --ld-wind-cm 1 --annot /scratch/ch6845/ukb_ld/annot/23115.chr1.annot --out --print-snps /scratch/ch6845/ukb_ld/hapmap3_snps/hm.1.snp\n",
    "    #ldsc.py --l2 --bfile /data01/ch6845/ukb_imp_plink_gwas/ukb_imp_chr9_v3 --ld-wind-cm 1 --annot /scratch/ch6845/ukb_ld/annot/23115.9.annot --out test2 --print-snps /scratch/ch6845/ukb_ld/hapmap3_snps/hm.9.snp\n",
    "    if mode=='annot':\n",
    "        if os.path.exists(ld_path.format(annot_name,chrN)+'.l2.ldscore.gz'):\n",
    "            print(\"Congratulations!. ld score file of\",annot_name,chrN,\"exists. passed.\")\n",
    "            return 0\n",
    "        scripts=['ldsc.py','--l2','--bfile',bfile_path.format(chrN), '--ld-wind-cm', '1','--annot',annot_path.format(annot_name,chrN),'--out',ld_path.format(annot_name,chrN),'--print-snps',print_snps_path.format(chrN)]\n",
    "    elif mode=='univariate':\n",
    "        annot_name='uni'\n",
    "        if os.path.exists(ld_path.format(annot_name,chrN)+'.l2.ldscore.gz'):\n",
    "            print(\"Congratulations!. ld score file of\",chrN,\"exists. passed.\")\n",
    "            return 0\n",
    "        scripts=['ldsc.py','--l2','--bfile',bfile_path.format(chrN), '--ld-wind-cm', '1','--out',ld_path.format(annot_name,chrN),'--print-snps',print_snps_path.format(chrN)]\n",
    "    print(' '.join(scripts))\n",
    "    sp.call(scripts)    \n",
    "    print(\"estimating ld score of \",annot_name,chrN,\"finished\")\n",
    "    return 0\n",
    "\n",
    "\n",
    "def estimate_ldscore(annot_name):\n",
    "    for chrN in range(1,22+1):\n",
    "        estimate_ldscore_chrN(chrN)\n",
    "    print(\"estimating LD score of\",annot_name,\"for all chr finished\")\n",
    "\n",
    "\n",
    "def estimate_ldscore_parallel(annot_name,mode='annot'):\n",
    "    pool = mp.Pool(processes=22)\n",
    "    if mode=='annot':\n",
    "        pool.starmap(estimate_ldscore_chrN,[[annot_name,chrN] for chrN in range(1,22+1)])\n",
    "    elif mode=='univariate':\n",
    "        pool.starmap(estimate_ldscore_chrN,[[annot_name,chrN,'univariate'] for chrN in range(1,22+1)])\n",
    "    print(\"estimating LD score of\",annot_name,\"finished\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "invalid literal for int() with base 10: '/run/user/1119/jupyter/kernel-617452a9-e175-4e42-8078-2060393890b1.json'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-e802b2571e46>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margv\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mscale\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margv\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: invalid literal for int() with base 10: '/run/user/1119/jupyter/kernel-617452a9-e175-4e42-8078-2060393890b1.json'"
     ]
    }
   ],
   "source": [
    "mode=sys.argv[1]\n",
    "scale=int(sys.argv[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "jupyter nbconvert 4_make_annot_ldscore.ipynb --to script && python 4_make_annot_ldscore.py 64 bp\n",
    "python 4_make_annot_ldscore.py bp 300\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "if scale!=0:\n",
    "    for chrN in range(1,22+1):\n",
    "        save_annot_chr_cm(mode+str(scale),chrN,scale=scale,mode=mode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr1 --ld-wind-cm 1 --out data/out_annot_/uni.1 --print-snps data/hapmap3_snps/hm.1.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr22 --ld-wind-cm 1 --out data/out_annot_/uni.22 --print-snps data/hapmap3_snps/hm.22.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr2 --ld-wind-cm 1 --out data/out_annot_/uni.2 --print-snps data/hapmap3_snps/hm.2.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr3 --ld-wind-cm 1 --out data/out_annot_/uni.3 --print-snps data/hapmap3_snps/hm.3.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr4 --ld-wind-cm 1 --out data/out_annot_/uni.4 --print-snps data/hapmap3_snps/hm.4.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr5 --ld-wind-cm 1 --out data/out_annot_/uni.5 --print-snps data/hapmap3_snps/hm.5.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr6 --ld-wind-cm 1 --out data/out_annot_/uni.6 --print-snps data/hapmap3_snps/hm.6.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr7 --ld-wind-cm 1 --out data/out_annot_/uni.7 --print-snps data/hapmap3_snps/hm.7.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr8 --ld-wind-cm 1 --out data/out_annot_/uni.8 --print-snps data/hapmap3_snps/hm.8.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr9 --ld-wind-cm 1 --out data/out_annot_/uni.9 --print-snps data/hapmap3_snps/hm.9.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr10 --ld-wind-cm 1 --out data/out_annot_/uni.10 --print-snps data/hapmap3_snps/hm.10.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr11 --ld-wind-cm 1 --out data/out_annot_/uni.11 --print-snps data/hapmap3_snps/hm.11.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr12 --ld-wind-cm 1 --out data/out_annot_/uni.12 --print-snps data/hapmap3_snps/hm.12.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr15 --ld-wind-cm 1 --out data/out_annot_/uni.15 --print-snps data/hapmap3_snps/hm.15.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr13 --ld-wind-cm 1 --out data/out_annot_/uni.13 --print-snps data/hapmap3_snps/hm.13.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr16 --ld-wind-cm 1 --out data/out_annot_/uni.16 --print-snps data/hapmap3_snps/hm.16.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr17 --ld-wind-cm 1 --out data/out_annot_/uni.17 --print-snps data/hapmap3_snps/hm.17.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr14 --ld-wind-cm 1 --out data/out_annot_/uni.14 --print-snps data/hapmap3_snps/hm.14.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr19 --ld-wind-cm 1 --out data/out_annot_/uni.19 --print-snps data/hapmap3_snps/hm.19.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr18 --ld-wind-cm 1 --out data/out_annot_/uni.18 --print-snps data/hapmap3_snps/hm.18.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr20 --ld-wind-cm 1 --out data/out_annot_/uni.20 --print-snps data/hapmap3_snps/hm.20.snp\n",
      "ldsc.py --l2 --bfile data/1000G_plink_EUR/EUR_phase3_chr21 --ld-wind-cm 1 --out data/out_annot_/uni.21 --print-snps data/hapmap3_snps/hm.21.snp\n",
      "estimating ld score of  uni 21 finished\n",
      "estimating ld score of  uni 22 finished\n",
      "estimating ld score of  uni 20 finished\n",
      "estimating ld score of  uni 19 finished\n",
      "estimating ld score of  uni 18 finished\n",
      "estimating ld score of  uni 17 finished\n",
      "estimating ld score of  uni 15 finished\n",
      "estimating ld score of  uni 16 finished\n",
      "estimating ld score of  uni 14 finished\n",
      "estimating ld score of  uni 13 finished\n",
      "estimating ld score of  uni 9 finished\n",
      "estimating ld score of  uni 12 finished\n",
      "estimating ld score of  uni 10 finished\n",
      "estimating ld score of  uni 11 finished\n",
      "estimating ld score of  uni 7 finished\n",
      "estimating ld score of  uni 8 finished\n",
      "estimating ld score of  uni 5 finished\n",
      "estimating ld score of  uni 1 finished\n",
      "estimating ld score of  uni 6 finished\n",
      "estimating ld score of  uni 4 finished\n",
      "estimating ld score of  uni 3 finished\n",
      "estimating ld score of  uni 2 finished\n",
      "estimating LD score of uni finished\n"
     ]
    }
   ],
   "source": [
    "if scale==0:\n",
    "    estimate_ldscore_parallel('uni','univariate')\n",
    "else:\n",
    "    estimate_ldscore_parallel(mode+str(scale),'annot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "for i in {1..22}; do gzip -d \"./bp300.${i}.l2.ldscore.gz\"; sed -i -e \"s/CHR\\tSNP\\tBP\\tL2/CHR\\tSNP\\tBP\\t${i}.0L2/g\" \"./bp300.${i}.l2.ldscore\"; gzip \"./bp300.${i}.l2.ldscore\" ; done\n",
    "\n",
    "for i in {1..22}; do gzip -d \"./bp128.${i}.l2.ldscore.gz\"; sed -i -e \"s/CHR\\tSNP\\tBP\\tL2/CHR\\tSNP\\tBP\\t${i}.0L2/g\" \"./bp128.${i}.l2.ldscore\"; gzip \"./bp128.${i}.l2.ldscore\" ; done\n",
    "\n",
    "for i in {1..22}; do gzip -d \"./bp64.${i}.l2.ldscore.gz\"; sed -i -e \"s/CHR\\tSNP\\tBP\\tL2/CHR\\tSNP\\tBP\\t${i}.0L2/g\" \"./bp64.${i}.l2.ldscore\"; gzip \"./bp64.${i}.l2.ldscore\" ; done\n",
    "\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}