gwas_3_scz_original_vs_ours.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    index   Rank                snp   A12  frq_case  frq_control chr  \\\n",
      "0       1   57.0     chr1_8424984_D   12D    0.3190       0.3010   1   \n",
      "1      13   78.0   chr1_243881945_I   12D    0.6380       0.6190   1   \n",
      "2      17   62.0   chr2_146436222_I   12D    0.1760       0.1630   2   \n",
      "3      18   95.0   chr2_149429178_D   12D    0.9550       0.9610   2   \n",
      "4      24   10.0   chr2_200825237_I   12D    0.7410       0.7540   2   \n",
      "5      33   28.0   chr3_180594593_I   I2D    0.1960       0.2080   3   \n",
      "6      45  128.0   chr5_140143664_I  I12D    0.4860       0.4750   5   \n",
      "7      53   53.0    chr6_84280274_D   12D    0.5240       0.5050   6   \n",
      "8      55    7.0     chr7_2025096_I   DI3    0.4050       0.4230   7   \n",
      "9      56  109.0    chr7_24747494_D   DI3    0.1040       0.0959   7   \n",
      "10     74   13.0  chr10_104957618_I   120    0.0654       0.0756  10   \n",
      "11     76   26.0   chr11_46350213_D   12D    0.8350       0.8500  11   \n",
      "12    111   43.0   chr18_52749216_D   I2D    0.5890       0.5670  18   \n",
      "13    121   41.0   chr22_39987017_D   I2D    0.5910       0.5940  22   \n",
      "\n",
      "                        pos left_border_prev right_border_prev  left_border  \\\n",
      "0       8,411,184-8,638,984        8,411,184         8,638,984      8411184   \n",
      "1   243,690,945-244,002,945      243,690,945       244,002,945    243690945   \n",
      "2   146,416,922-146,441,832      146,416,922       146,441,832    146416922   \n",
      "3   149,390,778-149,520,178      149,390,778       149,520,178    149390778   \n",
      "4   200,715,237-200,848,037      200,715,237       200,848,037    200715237   \n",
      "5   180,588,843-180,954 593      180,588,843       180,954 593    180588843   \n",
      "6   140,023,664-140,222,664      140,023,664       140,222,664    140023664   \n",
      "7     84,279,922-84,407,274       84,279,922        84,407,274     84279922   \n",
      "8       1,896,096-2,190,096        1,896,096         2,190,096      1896096   \n",
      "9     24,619,494-24,832,094       24,619,494        24,832,094     24619494   \n",
      "10  104,957,618-104,957,618      104,957,618       104,957,618    104957618   \n",
      "11    46,342,943-46,751,213       46,342,943        46,751,213     46342943   \n",
      "12    52,747,686-52,752,696       52,747,686        52,752,696     52747686   \n",
      "13    39,975,317-40,016,817       39,975,317        40,016,817     39975317   \n",
      "\n",
      "    right_border pos_actual  \n",
      "0      8638984.0    8424984  \n",
      "1    244002945.0  243881945  \n",
      "2    146441832.0  146436222  \n",
      "3    149520178.0  149429178  \n",
      "4    200848037.0  200825237  \n",
      "5    180954593.0  180594593  \n",
      "6    140222664.0  140143664  \n",
      "7     84407274.0   84280274  \n",
      "8      2190096.0    2025096  \n",
      "9     24832094.0   24747494  \n",
      "10   104957618.0  104957618  \n",
      "11    46751213.0   46350213  \n",
      "12    52752696.0   52749216  \n",
      "13    40016817.0   39987017  \n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "original_df = pd.read_csv(\"gwas_3_scz_original_files/gwas_3_scz_supp_table_2.csv\")\n",
    "\n",
    "original_df['left_border'] = original_df['left_border'] .str.replace(' ', '', regex=False)\n",
    "original_df['right_border']  = original_df['right_border'] .str.replace(' ', '', regex=False)\n",
    "original_df['left_border'] = pd.to_numeric(original_df['left_border'])\n",
    "original_df['right_border'] = pd.to_numeric(original_df['right_border'])\n",
    "\n",
    "starts_with_chr = original_df[original_df[\"snp\"].astype(str).str.startswith('chr')].reset_index()\n",
    "starts_with_chr['pos_actual'] = starts_with_chr['snp'].str.extract(r'_(\\d+)_')\n",
    "print(starts_with_chr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          snp_x   chr        pos_x  index   Rank              snp_y   A12  \\\n",
      "0    rs34269918   1.0    8424984.0      1   57.0     chr1_8424984_D   12D   \n",
      "1     rs5782266   1.0  243881945.0     13   78.0   chr1_243881945_I   12D   \n",
      "2    rs56807175   2.0  146436222.0     17   62.0   chr2_146436222_I   12D   \n",
      "3   rs200327371   2.0  149429178.0     18   95.0   chr2_149429178_D   12D   \n",
      "4           NaN   2.0          NaN     24   10.0   chr2_200825237_I   12D   \n",
      "5    rs11411529   3.0  180594593.0     33   28.0   chr3_180594593_I   I2D   \n",
      "6   rs111896713   5.0  140143664.0     45  128.0   chr5_140143664_I  I12D   \n",
      "7   rs540827191   6.0   84280274.0     53   53.0    chr6_84280274_D   12D   \n",
      "8    rs10650434   7.0    2025096.0     55    7.0     chr7_2025096_I   DI3   \n",
      "9   rs149009306   7.0   24747494.0     56  109.0    chr7_24747494_D   DI3   \n",
      "10          NaN  10.0          NaN     74   13.0  chr10_104957618_I   120   \n",
      "11   rs61126341  11.0   46350213.0     76   26.0   chr11_46350213_D   12D   \n",
      "12    rs5825114  18.0   52749216.0    111   43.0   chr18_52749216_D   I2D   \n",
      "13          NaN  22.0          NaN    121   41.0   chr22_39987017_D   I2D   \n",
      "\n",
      "    frq_case  frq_control                    pos_y left_border_prev  \\\n",
      "0     0.3190       0.3010      8,411,184-8,638,984        8,411,184   \n",
      "1     0.6380       0.6190  243,690,945-244,002,945      243,690,945   \n",
      "2     0.1760       0.1630  146,416,922-146,441,832      146,416,922   \n",
      "3     0.9550       0.9610  149,390,778-149,520,178      149,390,778   \n",
      "4     0.7410       0.7540  200,715,237-200,848,037      200,715,237   \n",
      "5     0.1960       0.2080  180,588,843-180,954 593      180,588,843   \n",
      "6     0.4860       0.4750  140,023,664-140,222,664      140,023,664   \n",
      "7     0.5240       0.5050    84,279,922-84,407,274       84,279,922   \n",
      "8     0.4050       0.4230      1,896,096-2,190,096        1,896,096   \n",
      "9     0.1040       0.0959    24,619,494-24,832,094       24,619,494   \n",
      "10    0.0654       0.0756  104,957,618-104,957,618      104,957,618   \n",
      "11    0.8350       0.8500    46,342,943-46,751,213       46,342,943   \n",
      "12    0.5890       0.5670    52,747,686-52,752,696       52,747,686   \n",
      "13    0.5910       0.5940    39,975,317-40,016,817       39,975,317   \n",
      "\n",
      "   right_border_prev  left_border  right_border  pos_actual  \n",
      "0          8,638,984      8411184     8638984.0     8424984  \n",
      "1        244,002,945    243690945   244002945.0   243881945  \n",
      "2        146,441,832    146416922   146441832.0   146436222  \n",
      "3        149,520,178    149390778   149520178.0   149429178  \n",
      "4        200,848,037    200715237   200848037.0   200825237  \n",
      "5        180,954 593    180588843   180954593.0   180594593  \n",
      "6        140,222,664    140023664   140222664.0   140143664  \n",
      "7         84,407,274     84279922    84407274.0    84280274  \n",
      "8          2,190,096      1896096     2190096.0     2025096  \n",
      "9         24,832,094     24619494    24832094.0    24747494  \n",
      "10       104,957,618    104957618   104957618.0   104957618  \n",
      "11        46,751,213     46342943    46751213.0    46350213  \n",
      "12        52,752,696     52747686    52752696.0    52749216  \n",
      "13        40,016,817     39975317    40016817.0    39987017  \n"
     ]
    }
   ],
   "source": [
    "first_file_df = pd.read_csv('./gwas_3_scz_matching/result_SAD0.csv', low_memory = False)\n",
    "starts_with_chr['chr'] = starts_with_chr['chr'].astype(int)\n",
    "starts_with_chr['pos_actual'] = starts_with_chr['pos_actual'].astype(int)\n",
    "merged_ids = pd.merge( first_file_df[['snp', 'chr', 'pos']], starts_with_chr, left_on=['chr', 'pos'], right_on = ['chr', 'pos_actual'] , how='right')\n",
    "print(merged_ids)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Rank             snp  A12  frq_case  frq_control chr  \\\n",
      "0  54.0       rs4648845   TC     0.533        0.527   1   \n",
      "1  57.0  chr1_8424984_D  12D     0.319        0.301   1   \n",
      "2  65.0       rs1498232   TC     0.311        0.296   1   \n",
      "3  50.0      rs11210892   AG     0.659        0.677   1   \n",
      "4  22.0      rs12129573   AC     0.377        0.358   1   \n",
      "\n",
      "                      pos left_border_prev right_border_prev  left_border  \\\n",
      "0         2372401-2402501          2372401           2402501      2372401   \n",
      "1     8,411,184-8,638,984        8,411,184         8,638,984      8411184   \n",
      "2   30,412,551-30,437,271       30,412,551        30,437,271     30412551   \n",
      "3  44,029,384-44, 128,084       44,029,384       44, 128,084     44029384   \n",
      "4   73,766,426-73,991 366       73,766,426        73,991 366     73766426   \n",
      "\n",
      "   right_border          P  \n",
      "0     2402501.0  4.033e-09  \n",
      "1     8638984.0        NaN  \n",
      "2    30437271.0  1.284e-09  \n",
      "3    44128084.0   4.97e-10  \n",
      "4    73991366.0  2.346e-10  \n",
      "   Rank         snp  A12  frq_case  frq_control chr                     pos  \\\n",
      "0  54.0   rs4648845   TC     0.533        0.527   1         2372401-2402501   \n",
      "1  57.0  rs34269918  12D     0.319        0.301   1     8,411,184-8,638,984   \n",
      "2  65.0   rs1498232   TC     0.311        0.296   1   30,412,551-30,437,271   \n",
      "3  50.0  rs11210892   AG     0.659        0.677   1  44,029,384-44, 128,084   \n",
      "4  22.0  rs12129573   AC     0.377        0.358   1   73,766,426-73,991 366   \n",
      "\n",
      "  left_border_prev right_border_prev  left_border  right_border             P  \n",
      "0          2372401           2402501      2372401     2402501.0  4.033000e-09  \n",
      "1        8,411,184         8,638,984      8411184     8638984.0           NaN  \n",
      "2       30,412,551        30,437,271     30412551    30437271.0  1.284000e-09  \n",
      "3       44,029,384       44, 128,084     44029384    44128084.0  4.970000e-10  \n",
      "4       73,766,426        73,991 366     73766426    73991366.0  2.346000e-10  \n"
     ]
    }
   ],
   "source": [
    "def replace_chr_with_snp_x(row, merged_ids):\n",
    "    if row.startswith('chr'):\n",
    "        # Lookup the corresponding 'snp_x' value in merged_ids DataFrame\n",
    "        replacement = merged_ids.loc[merged_ids['snp_y'] == row, 'snp_x']\n",
    "        if not replacement.empty:\n",
    "            return replacement.iloc[0]\n",
    "    return row\n",
    "\n",
    "# Create a dictionary from first_file_df\n",
    "p_value_dict = first_file_df.set_index('snp')['P'].to_dict()\n",
    "\n",
    "# Map the P values to the original_df using the snp column\n",
    "original_df['P'] = 0\n",
    "\n",
    "original_df['P'] = original_df['snp'].map(p_value_dict)\n",
    "print(original_df.head())\n",
    "\n",
    "original_df['snp'] = original_df['snp'].astype(str).apply(replace_chr_with_snp_x, merged_ids=merged_ids)\n",
    "original_df[\"P\"] = pd.to_numeric(original_df[\"P\"])\n",
    "print(original_df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_df.head()\n",
    "\n",
    "def merge_loci(df):\n",
    "    # Sort the DataFrame by chromosome and left_border\n",
    "    df = df.sort_values(by=['chr', 'left_border']).reset_index(drop=True)\n",
    "    merged_loci = []\n",
    "    current_locus = df.iloc[0].copy()  # Make a copy to avoid modifying the original DataFrame\n",
    "    \n",
    "    for i in range(1, len(df)):\n",
    "        next_locus = df.iloc[i]\n",
    "        \n",
    "        if (current_locus['chr'] == next_locus['chr']) and \\\n",
    "        (next_locus['left_border'] <= current_locus['right_border'] + 250000):\n",
    "            \n",
    "            # Update borders\n",
    "            current_locus['left_border'] = min(current_locus['left_border'], next_locus['left_border'])\n",
    "            current_locus['right_border'] = max(current_locus['right_border'], next_locus['right_border'])\n",
    "            \n",
    "            # Update SNP and P value if next_locus has a lower P value\n",
    "            if next_locus['P'] < current_locus['P']:\n",
    "                current_locus['snp'] = next_locus['snp']\n",
    "                current_locus['P'] = next_locus['P']\n",
    "        \n",
    "        else:\n",
    "            # Append the merged locus to the list\n",
    "            merged_loci.append(current_locus)\n",
    "            current_locus = next_locus.copy()  # Start a new current locus\n",
    "    \n",
    "    # Append the last locus\n",
    "    merged_loci.append(current_locus)\n",
    "    \n",
    "    # Convert the list of merged loci back to a DataFrame\n",
    "    merged_df = pd.DataFrame(merged_loci)\n",
    "    \n",
    "    return merged_df\n",
    "\n",
    "original_df_loci = merge_loci(original_df)\n",
    "original_df_loci.to_csv('gwas_3_scz_intermediate_files/original_df_merging')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95\n"
     ]
    }
   ],
   "source": [
    "our_df = pd.read_csv('./gwas_3_scz_result_files/filtered_snps_sd=0.0.csv')\n",
    "our_df.head()\n",
    "print(len(our_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "109\n",
      "95\n",
      "Number of potential matches: 104\n"
     ]
    }
   ],
   "source": [
    "def count_overlaps(df1, df2):\n",
    "        \n",
    "        overlap_count = 0\n",
    "    \n",
    "        # Iterate over all pairs of intervals\n",
    "        for _, row1 in df1.iterrows():\n",
    "            chr1 = int(row1['chr'])\n",
    "            for _, row2 in df2.iterrows():\n",
    "                chr2 = int(row2['chr'])\n",
    "                # Check if the intervals overlap on the same chromosome\n",
    "                if chr1==chr2 and row1['right_border'] >= row2['left_border'] and row1['left_border'] <= row2['right_border']:\n",
    "                    overlap_count += 1\n",
    "                    \n",
    "        return overlap_count\n",
    "\n",
    "original_df_loci[\"chr\"] = original_df_loci[\"chr\"].replace(\"X\", 23)\n",
    "print(len(original_df_loci))\n",
    "print(len(our_df))\n",
    "\n",
    "num_matches = original_df_loci['snp'].isin(first_file_df['snp']).sum()\n",
    "print(f\"Number of potential matches: {num_matches}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "88"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_overlaps(original_df_loci, our_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'rs117074560', 'rs2068012', 'rs11210892', 'rs1702294', 'rs9607782', 'rs7432375', 'rs6984242', 'rs6065094', 'rs950169', 'rs7819570', 'rs12826178', 'rs1106568', 'rs10791097', 'rs12129573', 'rs2514218', 'rs10503253', 'rs75968099', 'rs6704641', 'rs4648845', 'rs2332700', 'rs11411529', 'rs13240464', 'rs1498232', 'rs4129585', 'rs3735025', 'rs77149735', 'rs36068923', 'rs73229090', 'rs77502336', 'rs11027857', 'rs3849046', 'rs2973155', 'rs12325245', 'rs6434928', 'rs16867576', 'rs2851447', 'rs6704768', 'rs7801375', 'rs6002655', 'rs2693698', 'rs7893279', 'rs4766428', 'rs59979824', 'rs4523957', 'rs2053079', 'rs11682175', 'rs12903146', 'rs11693094', 'rs8044995', 'rs2007044', 'rs12887734', 'rs7405404', 'rs2535627', 'rs10520163', 'rs12704290', 'rs34269918', 'rs9922678', 'rs12691307', 'rs75059851', 'rs8042374', 'rs6466055', 'rs1501357', 'rs56807175', 'rs4702', 'rs17194490', 'rs11139497', 'rs215411', 'rs4391122', 'rs55661361', 'rs8082590', 'rs56205728', 'rs10650434', 'rs200327371', 'rs11685299', 'rs11740474', 'rs140505938', 'rs2905426', 'rs832187', 'rs2909457', 'rs61126341'}\n",
      "80\n"
     ]
    }
   ],
   "source": [
    "print(set(our_df['snp']).intersection(set(original_df_loci['snp'])))\n",
    "print(len(set(our_df['snp']).intersection(set(original_df_loci['snp']))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "newenv",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}