5_4_manhattan.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\n\\njupyter nbconvert 5_association.ipynb --to script\\npython 5_association.py albumin 2 0\\n\\nfor i in {00..101};do python 5_association.py $i 2 1;done\\n\\n\\nfor i in {00..10};do python 5_association.py $i 2 0;done\\nfor i in {11..20};do python 5_association.py $i 2 0;done\\nfor i in {21..30};do python 5_association.py $i 2 0;done\\nfor i in {31..40};do python 5_association.py $i 2 0;done\\nfor i in {41..50};do python 5_association.py $i 2 0;done\\nfor i in {51..60};do python 5_association.py $i 2 0;done\\nfor i in {61..70};do python 5_association.py $i 2 0;done\\nfor i in {71..80};do python 5_association.py $i 2 0;done\\nfor i in {91..90};do python 5_association.py $i 2 0;done\\nfor i in {91..101};do python 5_association.py $i 2 0;done\\nfor i in {51..60};do python 5_association.py $i 2 0;done\\nfor i in {51..60};do python 5_association.py $i 2 0;done\\nfor i in {51..60};do python 5_association.py $i 2 0;done\\n\\nfor i in {00..101};do python 5_association.py $i 1 0,2;done\\nfor i in {00..101};do python 5_association.py $i 1 1;done\\n\\nfor i in {00..20};do python 5_association.py $i 1 1;done\\nfor i in {21..40};do python 5_association.py $i 1 1;done\\nfor i in {41..60};do python 5_association.py $i 1 1;done\\nfor i in {61..80};do python 5_association.py $i 1 1;done\\nfor i in {81..101};do python 5_association.py $i 1 1;done\\n\\n'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "import logging\n",
    "import re\n",
    "import pathlib\n",
    "\n",
    "import datetime\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from pyplink import PyPlink\n",
    "import seaborn as sns\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "import matplotlib.patches as patches\n",
    "import matplotlib\n",
    "\n",
    "from basic_tools import *\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "jupyter nbconvert 5_4_manhattan.ipynb --to script\n",
    "\n",
    "for i in {00..20};do python 5_4_manhattan.py $i;done\n",
    "for i in {21..40};do python 5_4_manhattan.py $i;done\n",
    "for i in {41..60};do python 5_4_manhattan.py $i;done\n",
    "for i in {61..80};do python 5_4_manhattan.py $i;done\n",
    "for i in {81..98};do python 5_4_manhattan.py $i;done\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plink_KCHIP_HLA_AA_SNP_1000G=PyPlink(plink_KCHIP_HLA_AA_SNP_1000G_path)\n",
    "plink_KCHIP_HLA_AA_SNP_1000G_fam=plink_KCHIP_HLA_AA_SNP_1000G.get_fam().astype({'fid':str,'iid':str}).rename(columns={'fid':'FID','iid':'IID'})\n",
    "plink_KCHIP_HLA_AA_SNP_1000G_bim=plink_KCHIP_HLA_AA_SNP_1000G.get_bim()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "97\n"
     ]
    }
   ],
   "source": [
    "phenotypes=pd.read_csv(pheno_all_file_path,sep='\\t')\n",
    "phenotypes=phenotypes.set_index('ID').loc[plink_KCHIP_HLA_AA_SNP_1000G_fam['IID']]\n",
    "\n",
    "binary_continuous_traits=phenotypes.columns.difference(['age','sex','cohort','diabetes'])\n",
    "print(len(binary_continuous_traits))\n",
    "\n",
    "if 'ipykernel' in sys.argv[0]:\n",
    "    ipykernel=True\n",
    "    phenotype_name='allergic_disease'\n",
    "    #phenotype_name='height'\n",
    "else:\n",
    "    ipykernel=False\n",
    "    phenotype_name=sys.argv[1]\n",
    "    \n",
    "if phenotype_name.isdigit():\n",
    "    phenotype_name=int(phenotype_name)\n",
    "    phenotype_name=binary_continuous_traits[phenotype_name]      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_out_assoc_phenotype_path=data_out_assoc_path+phenotype_name+'/'\n",
    "pathlib.Path(data_out_assoc_phenotype_path).mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'binary'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(209, 30)\n",
      "(207, 30)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ch6845/tools/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:14: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "gene_bed_path='data/mart_export_grch37.txt'\n",
    "gene_bed=pd.read_csv(gene_bed_path,sep='\\t')\n",
    "gene_bed=gene_bed.drop(columns='Exon stable ID')\n",
    "gene_bed=gene_bed[(gene_bed['Gene start (bp)']>=plink_KCHIP_HLA_AA_SNP_1000G_bim.pos.min())&(gene_bed['Gene end (bp)']<=plink_KCHIP_HLA_AA_SNP_1000G_bim.pos.max())]\n",
    "gene_bed=gene_bed[(gene_bed['Transcript type']=='protein_coding')]\n",
    "gene_bed=gene_bed[~gene_bed.duplicated(['Gene name','Gene start (bp)','Gene end (bp)'])]\n",
    "print(gene_bed.shape)\n",
    "gene_bed=gene_bed[~gene_bed.duplicated(['Gene name'])]\n",
    "print(gene_bed.shape)\n",
    "\n",
    "gene_assign=plink_KCHIP_HLA_AA_SNP_1000G_bim[['pos']]\n",
    "\n",
    "for idx,row in gene_bed.iterrows():\n",
    "    gene_assign[row['Gene name']]=0\n",
    "    \n",
    "for idx,row in gene_bed.iterrows():    \n",
    "    gene_assign[row['Gene name']][(gene_assign['pos']>=row['Gene start (bp)'])&(gene_assign['pos']<=row['Gene end (bp)'])]=1\n",
    "\n",
    "#gene_assign.columns=gene_assign.columns.str.replace('HLA-','HLA_')        \n",
    "    \n",
    "HLA_names=np.unique([i[0].split('_')[1] for i in plink_KCHIP_HLA_AA_SNP_1000G_bim[plink_KCHIP_HLA_AA_SNP_1000G_bim.index.str.contains('HLA_')].index.str.split('*')])\n",
    "\n",
    "for HLA_name in HLA_names:\n",
    "    gene_select=gene_assign[gene_assign.index.str.contains('HLA_'+HLA_name)|gene_assign.index.str.contains('SNPS_'+HLA_name)|gene_assign.index.str.contains('AA_'+HLA_name)]#print(gene_select.sort_values('pos').iloc[0],gene_select.sort_values('pos').iloc[-1])\n",
    "    HLA_name='HLA-{}'.format(HLA_name)\n",
    "    gene_assign[HLA_name][(gene_assign['pos']>=gene_select['pos'].min())&(gene_assign['pos']<=gene_select['pos'].max())]=1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "log = logging.getLogger('logger')\n",
    "log.setLevel(logging.DEBUG)\n",
    "\n",
    "log_file_name=datetime.datetime.now().strftime('%Y%m%d_%H%M%S')+'.log'\n",
    "log_file_path=data_out_assoc_phenotype_path+log_file_name\n",
    "fileHandler = logging.FileHandler(log_file_path)\n",
    "streamHandler = logging.StreamHandler()\n",
    "\n",
    "formatter = logging.Formatter('%(message)s')\n",
    "fileHandler.setFormatter(formatter)\n",
    "streamHandler.setFormatter(formatter)\n",
    "\n",
    "log.addHandler(fileHandler)\n",
    "log.addHandler(streamHandler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "log.info_head=lambda x: log.info('-'*int((100-len(x))/2)+x+'-'*int((100-len(x))/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "----------------phenotype_name: allergic_disease, phenotype_type:binary , Step : 1 ----------------\n",
      "/home/ch6845/tools/miniconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3049: DtypeWarning: Columns (12) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n",
      "----------------phenotype_name: allergic_disease, phenotype_type:binary , Step : 2 ----------------\n",
      "/home/ch6845/tools/miniconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3049: DtypeWarning: Columns (12) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n",
      "----------------phenotype_name: allergic_disease, phenotype_type:binary , Step : 3 ----------------\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "meets end 3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1440x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (20,5)\n",
    "plt.rcParams[\"font.size\"] = 15\n",
    "plt.rcParams['font.family']='Arial'\n",
    "\n",
    "\n",
    "for step_idx_sub in range(1,100):\n",
    "    log.info_head(\"phenotype_name: {}, phenotype_type:{} , Step : {} \".format(phenotype_name,phenotype_type,step_idx_sub))\n",
    "    plt.clf()\n",
    "    #print(step_idx_sub)\n",
    "    if os.path.exists(data_out_assoc_phenotype_path+'step_{:02d}.cond.stop'.format(step_idx_sub)):\n",
    "        print('meets end',step_idx_sub)\n",
    "        break\n",
    "    result_merge=pd.read_csv(data_out_assoc_phenotype_path+'step_{:02d}.merge.result.tsv'.format(step_idx_sub),sep='\\t')\n",
    "\n",
    "\n",
    "\n",
    "    ##########################################\n",
    "    #############Plotting#####################\n",
    "    ##########################################\n",
    "\n",
    "    data=result_merge[['marker_name','POS','P']].copy()\n",
    "\n",
    "\n",
    "\n",
    "    data['-log10_P']=-np.log10(data['P'])\n",
    "    data['size']=-np.log10(data['P'])*5\n",
    "    data['size']=350*(data['size']-data['size'].min())/(data['size'].max()-data['size'].min())+0.001\n",
    "\n",
    "    POS_MB_min=27.9\n",
    "    POS_MB_max=35.1\n",
    "\n",
    "\n",
    "    tick_pad=20\n",
    "    gene_annot_pad=data['-log10_P'].max()/30\n",
    "    gene_annot_height=data['-log10_P'].max()/30\n",
    "\n",
    "    data['POS_MB']=data['POS']/1000000\n",
    "    data['check']=False\n",
    "\n",
    "    data=data.sort_values('P',ascending=False).reset_index()\n",
    "\n",
    "    ax = plt.subplot(111)\n",
    "\n",
    "\n",
    "    data_select_list=[]\n",
    "\n",
    "    for gene_name in ['HLA-A','HLA-B','HLA-C']:\n",
    "        pos_min=gene_assign[gene_assign[gene_name]==1]['pos'].min()/1000000\n",
    "        pos_max=gene_assign[gene_assign[gene_name]==1]['pos'].max()/1000000\n",
    "        color='#fa3c14'\n",
    "\n",
    "        rect = patches.Rectangle((pos_min,-gene_annot_height-gene_annot_pad),pos_max-pos_min,gene_annot_height,linewidth=1,edgecolor=None,facecolor=color)\n",
    "        ax.add_patch(rect)\n",
    "        #plt.annotate(gene_name.replace('_','-'),xy=((pos_min+pos_max)/2,-5),xytext=((pos_min+pos_max)/2,-10),fontsize=3,ha='center',arrowprops=dict(arrowstyle=\"- >\",connectionstyle=\"arc3,rad=0\",ls='dashed'))#,bbox= dict(boxstyle=\"round,pad=0.3\", fc=(1,1,1,0.5), ec=\"black\", lw=0.3),size=15)\n",
    "\n",
    "        select=(data['POS_MB']>pos_min)&(data['POS_MB']<pos_max)&(data['check']==False)\n",
    "        data.loc[select,'check']=True\n",
    "        data_select=data.loc[select].copy()   \n",
    "        data_select['color']=color\n",
    "\n",
    "        data_select_list.append(data_select.sort_values('-log10_P'))\n",
    "\n",
    "\n",
    "\n",
    "    for gene_name in ['HLA-DPA1', 'HLA-DPB1','HLA-DQA1', 'HLA-DQB1','HLA-DRB1','HLA-DRA']:\n",
    "        pos_min=gene_assign[gene_assign[gene_name]==1]['pos'].min()/1000000\n",
    "        pos_max=gene_assign[gene_assign[gene_name]==1]['pos'].max()/1000000\n",
    "        color='#2850c8'\n",
    "\n",
    "        rect = patches.Rectangle((pos_min,-gene_annot_height-gene_annot_pad),pos_max-pos_min,gene_annot_height,linewidth=1,edgecolor=None,facecolor=color)\n",
    "        ax.add_patch(rect)\n",
    "        #plt.annotate(gene_name.replace('_','-'),xy=((pos_min+pos_max)/2,-5),xytext=((pos_min+pos_max)/2,-10),fontsize=3,ha='center',arrowprops=dict(arrowstyle=\"- >\",connectionstyle=\"arc3,rad=0\",ls='dashed'))#,bbox= dict(boxstyle=\"round,pad=0.3\", fc=(1,1,1,0.5), ec=\"black\", lw=0.3),size=15)\n",
    "\n",
    "        select=(data['POS_MB']>pos_min)&(data['POS_MB']<pos_max)&(data['check']==False)\n",
    "        data.loc[select,'check']=True\n",
    "        data_select=data.loc[select].copy()       \n",
    "        data_select['color']=color\n",
    "        data_select_list.append(data_select.sort_values('-log10_P'))\n",
    "\n",
    "\n",
    "    #for gene_name in gene_assign.columns[gene_assign.columns.str.contains('HLA')].difference(['HLA_A','HLA_B','HLA_C','HLA_DPA1','HLA_DPB1','HLA_DQA1', 'HLA_DQB1', 'HLA_DRB1','HLA-DRA']):#['HLA_A29.1']):\n",
    "    for gene_name in ['HLA-F', 'HLA-G', 'HLA-E', 'HLA-DRB5', 'HLA-DQA2', 'HLA-DQB2', 'HLA-DOB', 'HLA-DMB', 'HLA-DMA', 'HLA-DOA', 'MICA','MICB','TAP2','TAP1']:\n",
    "    #for gene_name in ['HLA-F', 'HLA-G',  'HLA-E', 'HLA-DRB5','HLA-DMB', 'HLA-DMA', 'HLA-DOA', 'HLA-DOB','HLA-DQA2', 'HLA-DQB2']:\n",
    "        pos_min=gene_assign[gene_assign[gene_name]==1]['pos'].min()/1000000\n",
    "        pos_max=gene_assign[gene_assign[gene_name]==1]['pos'].max()/1000000\n",
    "        color='#28c828'\n",
    "\n",
    "        rect = patches.Rectangle((pos_min,-gene_annot_height-gene_annot_pad),pos_max-pos_min,gene_annot_height,linewidth=1,edgecolor=None,facecolor=color)\n",
    "        ax.add_patch(rect)\n",
    "        #plt.annotate(gene_name.replace('_','-'),xy=((pos_min+pos_max)/2,-5),xytext=((pos_min+pos_max)/2,-10),fontsize=3,ha='center',arrowprops=dict(arrowstyle=\"- >\",connectionstyle=\"arc3,rad=0\",ls='dashed'))#,bbox= dict(boxstyle=\"round,pad=0.3\", fc=(1,1,1,0.5), ec=\"black\", lw=0.3),size=15)\n",
    "\n",
    "\n",
    "        select=(data['POS_MB']>pos_min)&(data['POS_MB']<pos_max)&(data['check']==False)\n",
    "        data.loc[select,'check']=True\n",
    "        data_select=data.loc[select].copy()       \n",
    "        data_select['color']=color\n",
    "        data_select_list.append(data_select.sort_values('-log10_P'))  \n",
    "\n",
    "\n",
    "    for gene_name in gene_assign.columns[~gene_assign.columns.str.contains('HLA')]:\n",
    "        #print(gene_name)\n",
    "        pos_min=gene_assign[gene_assign[gene_name]==1]['pos'].min()/1000000\n",
    "        pos_max=gene_assign[gene_assign[gene_name]==1]['pos'].max()/1000000\n",
    "        color='#fafac8'\n",
    "        #color=matplotlib.colors.rgb2hex(0.5*np.array(plt.cm.rainbow(np.random.rand())[:3]))\n",
    "        rect = patches.Rectangle((pos_min,-gene_annot_height-gene_annot_pad),pos_max-pos_min,gene_annot_height,linewidth=1,edgecolor='grey',facecolor=color,alpha=0.5)\n",
    "        #rect = patches.Rectangle((pos_min,-5.5),pos_max-pos_min,3,linewidth=1,edgecolor='grey',facecolor=color,alpha=0.5)\n",
    "        ax.add_patch(rect)\n",
    "        #plt.annotate(gene_name.replace('_','-'),xy=((pos_min+pos_max)/2,-5),xytext=((pos_min+pos_max)/2,-10),fontsize=3,ha='center',arrowprops=dict(arrowstyle=\"- >\",connectionstyle=\"arc3,rad=0\",ls='dashed'))#,bbox= dict(boxstyle=\"round,pad=0.3\", fc=(1,1,1,0.5), ec=\"black\", lw=0.3),size=15)\n",
    "\n",
    "\n",
    "        select=(data['POS_MB']>pos_min)&(data['POS_MB']<pos_max)&(data['check']==False)\n",
    "        data.loc[select,'check']=True\n",
    "        data_select=data.loc[select].copy()       \n",
    "        data_select['color']=color\n",
    "        data_select_list.append(data_select.sort_values('-log10_P'))     \n",
    "\n",
    "\n",
    "    color='#fafac8'\n",
    "    select=(data['check']==False)\n",
    "    data.loc[select,'check']=True\n",
    "    data_select=data.loc[select].copy()    \n",
    "    data_select['color']=color\n",
    "    data_select_list.append(data_select.sort_values('-log10_P'))  \n",
    "    #data_select=data.loc[select].sort_values('-log10_P')\n",
    "    #plt.scatter(x=data_select['POS_MB'],y=data_select['-log10_P'],sizes=data_select['size'],marker='d',linewidth=0.3,edgecolor='black',color=color)   \n",
    "\n",
    "    for data_select in data_select_list[::-1]:\n",
    "        plt.scatter(x=data_select['POS_MB'],y=data_select['-log10_P'],sizes=data_select['size'],marker='d',linewidth=0.3,edgecolor='black',color=data_select['color'])   \n",
    "\n",
    "\n",
    "\n",
    "    plt.plot(np.linspace(POS_MB_min,POS_MB_max,100),-np.log10(5e-8)*np.ones_like(np.linspace(POS_MB_min,POS_MB_max,100)),'--',color='black',alpha=0.7,linewidth=0.7)    \n",
    "    #ax.spines['left'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)  \n",
    "    ax.spines['bottom'].set_position('zero')#ax.spines['bottom'].set_visible(False) \n",
    "\n",
    "    ax.tick_params(axis='x', which='major', pad=tick_pad)\n",
    "\n",
    "    plt.xlim(POS_MB_min,POS_MB_max)\n",
    "    #ax.set_ylim(bottom=0.)\n",
    "    #ax.yaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "    #plt.savefig('savefig_polygenicity.png',bbox_inches=\"tight\")        \n",
    "\n",
    "\n",
    "    ##########################################\n",
    "    #############Plotting#####################\n",
    "    ##########################################\n",
    "    #plt.savefig('savefig_polygenicity.png',bbox_inches=\"tight\") \n",
    "    #plt.show()\n",
    "    plt.savefig(data_out_assoc_phenotype_path+'step_{:02d}.merge.manhattan.png'.format(step_idx_sub),bbox_inches=\"tight\")\n",
    "    #plt.show()\n",
    "    #break\n",
    "#break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "tar -czvf s_gwas_all.tar data/out_assoc/*/*merge*.png\n",
    "tar -czvf s_figure_manhattan.tar data/out_assoc/*/*merge*"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "#for HLA_name in gene_assign.columns[gene_assign.columns.str.contains('HLA')].difference(['HLA_A29.1']):\n",
    "#    gene_assign[gene_assign[HLA_name]==1]['pos'].min(),gene_assign[gene_assign[HLA_name]==1]['pos'].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ax.tick_params?"
   ]
  }
 ],
 "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": 4
}