7_1_gcta_bivar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "jupyter nbconvert 7_1_gcta_bivar.ipynb --to script\n",
    "\n",
    "\n",
    "for i in {00..20};do python 7_1_gcta_bivar.py $i;done\n",
    "for i in {21..40};do python 7_1_gcta_bivar.py $i;done\n",
    "for i in {41..60};do python 7_1_gcta_bivar.py $i;done\n",
    "for i in {61..80};do python 7_1_gcta_bivar.py $i;done\n",
    "for i in {81..101};do python 7_1_gcta_bivar.py $i;done\n",
    "\n",
    "\n",
    "\n",
    "python 7_1_gcta_bivar.py $WINDOW\n",
    "\n",
    "for i in {00..101};do python 6_gcta_uni.py $i;done\n",
    "\n",
    "for i in {00..10};do python 4_association.py $i;done\n",
    "\n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6292650/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "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",
    "\n",
    "from basic_tools import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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_gcta_path+log_file_name\n",
    "fileHandler = logging.FileHandler(log_file_path)\n",
    "streamHandler = logging.StreamHandler()\n",
    "\n",
    "formatter = logging.Formatter(' %(asctime)s [%(levelname)s] %(lineno)d > %(message)s')\n",
    "fileHandler.setFormatter(formatter)\n",
    "streamHandler.setFormatter(formatter)\n",
    "\n",
    "log.addHandler(fileHandler)\n",
    "log.addHandler(streamHandler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# load plink, aa and check integrity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "grm_path='data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#final_plink_aa_grm_path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# load phenotype and check integrity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ID  age  sex  asthma  chronic_bronchitis  allergic_disease  stomach_cancer  liver_cancer  colorectal_cancer  breast_cancer  cervical_cancer  lung_cancer  thyroid_cancer  prostate_cancer  bladder_cancer  colon_polyps  angina_pectoris_myocardial_infarction  hypertension  stroke  transient_ischemic_attacks  acute_liver_disease  diabetes  t2_diabetes  hyperlipidemia  thyroid_disease  arthritis  rheumatoid_arthritis  osteoporosis  parkinson  pwi_total  chronic_gastritis_stomach_ulcer  duodenal_ulcer  urinary_tract_infection  height  hip_circumference  waist_circumference  weight  systolic_blood_pressure  diastolic_blood_pressure  pulse  hysterectomy  menarche_onset_age  menopause_age  menopause_cycle  ovarian_extraction  child_delivery_number  big_baby  small_baby  child_present  child_dead  hematocrit  hemoglobin  platelet  rbc  wbc  mcv  mch  mchc  total_bilirubin  bun  creatinine  blood_in_urine  glucose_in_urine  alt  ast  gtp  total_cholesterol  hdl_cholesterol  hba1c  freet4  tsh  triglyceride  albumin  hscrp  protein_in_urine  heart_failure  peripheral_vascular_diseases  cataract  glaucoma  renal_failure  renin  fatty_liver  copd  tuberculosis  benign_prostatic_hyperplasia  gout  fracture  cholelithiasis_cholecystitis  periodontal_disease  insomnia  sleep_duration  neurosis  grip_strength  FVC_predicted  FEV_predicted  FEV_over_FVC_predicted  MVV_predicted  bone_density_distal_sos  bone_density_midshaft_sos  protein_in_blood  cohort  bmi'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phenotypes=pd.read_csv(pheno_all_file_path,sep='\\t')\n",
    "'  '.join(phenotypes.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(125673, 101)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phenotypes=phenotypes.set_index('ID').loc[plink_KCHIP_HLA_AA_SNP_1000G_fam['IID']]\n",
    "phenotypes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#np.all(phenotypes['ALP'].isnull())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert (phenotypes.index!=plink_KCHIP_HLA_AA_SNP_1000G_fam['IID']).sum()==0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Index(['FEV_over_FVC_predicted', 'FEV_predicted', 'FVC_predicted',\n",
       "        'MVV_predicted', 'acute_liver_disease', 'albumin', 'allergic_disease',\n",
       "        'alt', 'angina_pectoris_myocardial_infarction', 'arthritis', 'ast',\n",
       "        'asthma', 'benign_prostatic_hyperplasia', 'big_baby', 'bladder_cancer',\n",
       "        'blood_in_urine', 'bmi', 'bone_density_distal_sos',\n",
       "        'bone_density_midshaft_sos', 'breast_cancer', 'bun', 'cataract',\n",
       "        'cervical_cancer', 'child_dead', 'child_delivery_number',\n",
       "        'child_present', 'cholelithiasis_cholecystitis', 'chronic_bronchitis',\n",
       "        'chronic_gastritis_stomach_ulcer', 'colon_polyps', 'colorectal_cancer',\n",
       "        'copd', 'creatinine', 'diabetes', 'diastolic_blood_pressure',\n",
       "        'duodenal_ulcer', 'fatty_liver', 'fracture', 'freet4', 'glaucoma',\n",
       "        'glucose_in_urine', 'gout', 'grip_strength', 'gtp', 'hba1c',\n",
       "        'hdl_cholesterol', 'heart_failure', 'height', 'hematocrit',\n",
       "        'hemoglobin', 'hip_circumference', 'hscrp', 'hyperlipidemia',\n",
       "        'hypertension', 'hysterectomy', 'insomnia', 'liver_cancer',\n",
       "        'lung_cancer', 'mch', 'mchc', 'mcv', 'menarche_onset_age',\n",
       "        'menopause_age', 'menopause_cycle', 'neurosis', 'osteoporosis',\n",
       "        'ovarian_extraction', 'parkinson', 'periodontal_disease',\n",
       "        'peripheral_vascular_diseases', 'platelet', 'prostate_cancer',\n",
       "        'protein_in_blood', 'protein_in_urine', 'pulse', 'pwi_total', 'rbc',\n",
       "        'renal_failure', 'renin', 'rheumatoid_arthritis', 'sleep_duration',\n",
       "        'small_baby', 'stomach_cancer', 'stroke', 'systolic_blood_pressure',\n",
       "        't2_diabetes', 'thyroid_cancer', 'thyroid_disease', 'total_bilirubin',\n",
       "        'total_cholesterol', 'transient_ischemic_attacks', 'triglyceride',\n",
       "        'tsh', 'tuberculosis', 'urinary_tract_infection', 'waist_circumference',\n",
       "        'wbc', 'weight'],\n",
       "       dtype='object'), 98)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_continuous_traits=phenotypes.columns.difference(['age','sex','cohort'])\n",
    "binary_continuous_traits,len(binary_continuous_traits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# parse parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 2020-07-09 04:14:04,122 [INFO] 9 > 0\n"
     ]
    }
   ],
   "source": [
    "if 'ipykernel' in sys.argv[0]:\n",
    "    ipykernel=True\n",
    "    i=0\n",
    "    #phenotype_name='height'\n",
    "else:\n",
    "    ipykernel=False\n",
    "    i=int(sys.argv[1])  \n",
    "    \n",
    "log.info(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#i=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#(pheno['pheno_x']!=-9).sum(),(pheno['pheno_y']!=-9).sum(),((pheno['pheno_x']!=-9) & (pheno['pheno_y']!=-9)).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pd.read_csv(data_out_gcta_path+phenotype_name1+'-'+phenotype_name2+'.phe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for j in range(i+1,len(binary_continuous_traits)):\n",
    "    phenotype_name1=binary_continuous_traits[i]\n",
    "    phenotype_name2=binary_continuous_traits[j]\n",
    "    \n",
    "    if os.path.exists(data_out_gcta_path+phenotype_name1+'-'+phenotype_name2+'.HEreg') or os.path.exists(data_out_gcta_path+phenotype_name2+'-'+phenotype_name1+'.HEreg'):\n",
    "        continue\n",
    "\n",
    "\n",
    "\n",
    "    pheno1=pd.read_csv(data_out_pheno_path+phenotype_name1+'.phe',sep='\\t',header=None,names=['FID','IID','pheno'])\n",
    "    pheno2=pd.read_csv(data_out_pheno_path+phenotype_name2+'.phe',sep='\\t',header=None,names=['FID','IID','pheno'])\n",
    "    pheno=pheno1.merge(right=pheno2,left_on=['FID','IID'],right_on=['FID','IID'])\n",
    "\n",
    "    pheno_filter=pheno[(pheno['pheno_x']!=-9) & (pheno['pheno_y']!=-9)]\n",
    "    pheno_filter.to_csv(data_out_gcta_path+phenotype_name1+'-'+phenotype_name2+'.phe',sep='\\t',index=None,header=None)\n",
    "\n",
    "    log.info(\"phenotype_name1: {}, phenotype_name2:{}\".format(phenotype_name1,phenotype_name2))\n",
    "    log.info('pheno1 not missing {}'.format((pheno['pheno_x']!=-9).sum()))\n",
    "    log.info('pheno2 not missing {}'.format((pheno['pheno_y']!=-9).sum()))\n",
    "\n",
    "    if os.path.exists(data_out_gcta_path+phenotype_name1+'-'+phenotype_name2+'.HEreg'):\n",
    "        log.info(data_out_gcta_path+phenotype_name1+'-'+phenotype_name2+'.HEreg exists')\n",
    "        break\n",
    "    \n",
    "    log.info(\"#########################################  Run GCTA  #########################################\")\n",
    "    #Run omnibus association test\n",
    "    command='gcta64 --HEreg-bivar 1 2 --grm {} --pheno {} --out {} --thread-num 10'.format(grm_path,\n",
    "                                                                                 data_out_gcta_path+phenotype_name1+'-'+phenotype_name2+'.phe',\n",
    "                                                                                 data_out_gcta_path+phenotype_name1+'-'+phenotype_name2,                                         \n",
    "                                                                                )\n",
    "    log.info(command)\n",
    "    stdout,stderr=run_subprocess(command,dry=False)\n",
    "    log.info(stdout)\n",
    "    log.error(stderr)          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}