6_1_gcta_uni.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "jupyter nbconvert 6_1_gcta_uni.ipynb --to script\n",
    "\n",
    "for i in {00..101};do python 6_1_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": "markdown",
   "metadata": {},
   "source": [
    "# load plink, aa and check integrity"
   ]
  },
  {
   "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": [],
   "source": [
    "grm_path='data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#final_plink_aa_grm_path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# load phenotype and check integrity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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": 5,
     "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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(125673, 101)"
      ]
     },
     "execution_count": 6,
     "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#np.all(phenotypes['ALP'].isnull())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert (phenotypes.index!=plink_KCHIP_HLA_AA_SNP_1000G_fam['IID']).sum()==0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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', 'diastolic_blood_pressure', 'duodenal_ulcer',\n",
       "        'fatty_liver', 'fracture', 'freet4', 'glaucoma', 'glucose_in_urine',\n",
       "        'gout', 'grip_strength', 'gtp', 'hba1c', 'hdl_cholesterol',\n",
       "        'heart_failure', 'height', 'hematocrit', 'hemoglobin',\n",
       "        'hip_circumference', 'hscrp', 'hyperlipidemia', 'hypertension',\n",
       "        'hysterectomy', 'insomnia', 'liver_cancer', 'lung_cancer', 'mch',\n",
       "        'mchc', 'mcv', 'menarche_onset_age', 'menopause_age', 'menopause_cycle',\n",
       "        'neurosis', 'osteoporosis', 'ovarian_extraction', 'parkinson',\n",
       "        'periodontal_disease', 'peripheral_vascular_diseases', 'platelet',\n",
       "        'prostate_cancer', 'protein_in_blood', 'protein_in_urine', 'pulse',\n",
       "        'pwi_total', 'rbc', 'renal_failure', 'renin', 'rheumatoid_arthritis',\n",
       "        'sleep_duration', 'small_baby', 'stomach_cancer', 'stroke',\n",
       "        'systolic_blood_pressure', 't2_diabetes', 'thyroid_cancer',\n",
       "        'thyroid_disease', 'total_bilirubin', 'total_cholesterol',\n",
       "        'transient_ischemic_attacks', 'triglyceride', 'tsh', 'tuberculosis',\n",
       "        'urinary_tract_infection', 'waist_circumference', 'wbc', 'weight'],\n",
       "       dtype='object'), 97)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_continuous_traits=phenotypes.columns.difference(['age','sex','cohort','diabetes'])\n",
    "binary_continuous_traits,len(binary_continuous_traits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# parse parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'ipykernel' in sys.argv[0]:\n",
    "    ipykernel=True\n",
    "    phenotype_name='diabetes'\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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'binary'"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pheno=pd.read_csv(data_out_pheno_path+phenotype_name+'.phe',sep='\\t',header=None,names=['FID','IID','pheno'])\n",
    "phenotype_type='binary' if len(pheno['pheno'][pheno['pheno']!=-9].value_counts())<3 else 'continuous'\n",
    "phenotype_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37\n",
      "freet4 False\n",
      "68\n",
      "peripheral_vascular_diseases False\n",
      "84\n",
      "t2_diabetes False\n"
     ]
    }
   ],
   "source": [
    "for a,i in enumerate(binary_continuous_traits):\n",
    "    if not os.path.exists(data_out_gcta_path+i+'.HEreg'):\n",
    "        print(a)\n",
    "        #print(a)\n",
    "        print(i,os.path.exists(data_out_gcta_path+i+'.HEreg'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 2020-07-09 03:57:01,697 [INFO] 1 > phenotype_name: diabetes, phenotype_type:binary\n"
     ]
    }
   ],
   "source": [
    "log.info(\"phenotype_name: {}, phenotype_type:{}\".format(phenotype_name,phenotype_type))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno[pheno['pheno']!=-9].to_csv(data_out_gcta_path+phenotype_name+'.phe',sep='\\t',index=None,header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prev 0.07085054068893079\n"
     ]
    }
   ],
   "source": [
    "if phenotype_type=='binary':\n",
    "    with open(data_out_pheno_path+phenotype_name+'.phe'+'.prev','r') as f:\n",
    "        prev=float(f.read())\n",
    "    print('prev',prev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 2020-07-09 03:57:04,018 [INFO] 1 > #########################################  Run GCTA  #########################################\n",
      " 2020-07-09 03:57:04,021 [INFO] 7 > gcta64 --HEreg --grm data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm --pheno data/out_gcta/diabetes.phe --out data/out_gcta/diabetes --thread-num 40\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------RUN-----\n",
      "gcta64 --HEreg --grm data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm --pheno data/out_gcta/diabetes.phe --out data/out_gcta/diabetes --thread-num 40\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 2020-07-09 04:03:10,455 [INFO] 9 > *******************************************************************\n",
      "* Genome-wide Complex Trait Analysis (GCTA)\n",
      "* version 1.92.2 beta Linux\n",
      "* (C) 2010-2019, The University of Queensland\n",
      "* Please report bugs to: Jian Yang <jian.yang@uq.edu.au>\n",
      "*******************************************************************\n",
      "Analysis started at 18:57:04 UTC on Wed Jul 08 2020.\n",
      "Hostname: h3\n",
      "\n",
      "Accepted options:\n",
      "--HEreg\n",
      "--grm data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm\n",
      "--pheno data/out_gcta/diabetes.phe\n",
      "--out data/out_gcta/diabetes\n",
      "--thread-num 40\n",
      "\n",
      "Note: the program will be running on 40 threads.\n",
      "\n",
      "Reading IDs of the GRM from [data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm.grm.id].\n",
      "125673 IDs read from [data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm.grm.id].\n",
      "Reading phenotypes from [data/out_gcta/diabetes.phe].\n",
      "Non-missing phenotypes of 89355 individuals are included from [data/out_gcta/diabetes.phe].\n",
      "89355 individuals are in common in these files.\n",
      "\n",
      "Performing Haseman-Elston regression ...\n",
      "\n",
      "Standardising the phenotype ...\n",
      "Constructing ordinary least squares equations ...\n",
      "\n",
      "Left-hand side of OLS equations (X'X)\n",
      "3.99211e+09     8642.88\n",
      "    8642.88 9.00677e+07\n",
      "\n",
      "HE-CP\n",
      "Coefficient     Estimate        SE_OLS          SE_Jackknife    P_OLS           P_Jackknife     \n",
      "Intercept       -1.11926e-05    1.58262e-05     2.00222e-07     0.479427        0               \n",
      "V(G)/Vp         0.000602419     0.000105364     0.00027065      1.08107e-08     0.0260262       \n",
      "\n",
      "HE-SD\n",
      "Coefficient     Estimate        SE_OLS          SE_Jackknife    P_OLS           P_Jackknife     \n",
      "Intercept       -1              3.38452e-05     0.00894295      0               0               \n",
      "V(G)/Vp         0.000558495     0.000225328     0.000520975     0.0131903       0.28371         \n",
      "\n",
      "Results from Haseman-Elston regression have been saved in [data/out_gcta/diabetes.HEreg].\n",
      "\n",
      "Analysis finished at 19:03:10 UTC on Wed Jul 08 2020\n",
      "Computational time: 6 minutes 6 second(s).\n",
      "\n",
      " 2020-07-09 04:03:10,457 [ERROR] 10 > \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------ERROR-----\n",
      "\n",
      "------OUTPUT-----\n",
      "*******************************************************************\n",
      "* Genome-wide Complex Trait Analysis (GCTA)\n",
      "* version 1.92.2 beta Linux\n",
      "* (C) 2010-2019, The University of Queensland\n",
      "* Please report bugs to: Jian Yang <jian.yang@uq.edu.au>\n",
      "*******************************************************************\n",
      "Analysis started at 18:57:04 UTC on Wed Jul 08 2020.\n",
      "Hostname: h3\n",
      "\n",
      "Accepted options:\n",
      "--HEreg\n",
      "--grm data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm\n",
      "--pheno data/out_gcta/diabetes.phe\n",
      "--out data/out_gcta/diabetes\n",
      "--thread-num 40\n",
      "\n",
      "Note: the program will be running on 40 threads.\n",
      "\n",
      "Reading IDs of the GRM from [data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm.grm.id].\n",
      "125673 IDs read from [data/genotype/4_merge/KCHIP_HLA_AA_SNP_1000G.grm.grm.id].\n",
      "Reading phenotypes from [data/out_gcta/diabetes.phe].\n",
      "Non-missing phenotypes of 89355 individuals are included from [data/out_gcta/diabetes.phe].\n",
      "89355 individuals are in common in these files.\n",
      "\n",
      "Performing Haseman-Elston regression ...\n",
      "\n",
      "Standardising the phenotype ...\n",
      "Constructing ordinary least squares equations ...\n",
      "\n",
      "Left-hand side of OLS equations (X'X)\n",
      "3.99211e+09     8642.88\n",
      "    8642.88 9.00677e+07\n",
      "\n",
      "HE-CP\n",
      "Coefficient     Estimate        SE_OLS          SE_Jackknife    P_OLS           P_Jackknife     \n",
      "Intercept       -1.11926e-05    1.58262e-05     2.00222e-07     0.479427        0               \n",
      "V(G)/Vp         0.000602419     0.000105364     0.00027065      1.08107e-08     0.0260262       \n",
      "\n",
      "HE-SD\n",
      "Coefficient     Estimate        SE_OLS          SE_Jackknife    P_OLS           P_Jackknife     \n",
      "Intercept       -1              3.38452e-05     0.00894295      0               0               \n",
      "V(G)/Vp         0.000558495     0.000225328     0.000520975     0.0131903       0.28371         \n",
      "\n",
      "Results from Haseman-Elston regression have been saved in [data/out_gcta/diabetes.HEreg].\n",
      "\n",
      "Analysis finished at 19:03:10 UTC on Wed Jul 08 2020\n",
      "Computational time: 6 minutes 6 second(s).\n",
      "\n"
     ]
    }
   ],
   "source": [
    "log.info(\"#########################################  Run GCTA  #########################################\")\n",
    "#Run omnibus association test\n",
    "command='gcta64 --HEreg --grm {} --pheno {} --out {} --thread-num 40'.format(grm_path,\n",
    "                                                                             data_out_gcta_path+phenotype_name+'.phe',\n",
    "                                                                             data_out_gcta_path+phenotype_name,                                         \n",
    "                                                                            )\n",
    "log.info(command)\n",
    "stdout,stderr=run_subprocess(command,dry=False)\n",
    "log.info(stdout)\n",
    "log.error(stderr)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09964747356051704"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pheno0=pd.read_csv(data_out_assoc_path+phenotype_list[0]+'/'+'phenotype.phe',header=None,sep='\\t',names=['FID','IID','pheno0'])\n",
    "#pheno1=pd.read_csv(data_out_assoc_path+phenotype_list[1]+'/'+'phenotype.phe',header=None,sep='\\t',names=['FID','IID','pheno1']);pheno0['pheno1']=pheno1['pheno1']\n",
    "#pheno2=pd.read_csv(data_out_assoc_path+phenotype_list[1]+'/'+'phenotype.phe',header=None,sep='\\t',names=['FID','IID','pheno2']);pheno0['pheno2']=pheno2['pheno2']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "command='gcta64 --HEreg-bivar 1 2 --grm {} --pheno {} --out {} --thread-num 40'.format(final_plink_aa_grm_path,\n",
    "                                                                            'temp.phe',\n",
    "                                                                             'testout2'  \n",
    "                                                                            )\n",
    "command"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "gcta64 --HEreg-bivar 1 2 --grm data/genotype/4_merge/grm --pheno temp2.phe --out data/out_assoc/height/HE2 --thread-num 40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(pheno0['pheno0']!=-9).sum(),(pheno0['pheno1']!=-9).sum(),(pheno0['pheno2']!=-9).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!cat data/out_gcta/ALP.phe"
   ]
  },
  {
   "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
}