Raylab_analysis_v2g_codecombinedforpaper_09142024.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8cfac8e-00b9-442b-9c6d-f782b0061f26",
   "metadata": {},
   "outputs": [],
   "source": [
    "##### rsID to var ID conversion using genopyc tool "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90bae75b-04c0-4097-9eca-c95ba87719fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import genopyc as gp \n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dd6cbbf-7051-4442-91f6-835c326ce1b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"/nfs/amishra/Arpit/rsID18k_onecol.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a88ec086-4010-46ff-b5a3-de512ad6d2d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[[\"rsid\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d48558e8-ca11-4d13-a5f8-1bbee4e1213f",
   "metadata": {},
   "outputs": [],
   "source": [
    "variants_to_study_ot = gp.variantId_mapping(df.rsid, source='rsid',target= 'variantid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c54d0c4c-a0fb-4d22-87ec-11731f26716d",
   "metadata": {},
   "outputs": [],
   "source": [
    "variantdf = pd.DataFrame(variants_to_study_ot)\n",
    "variantdf.columns = [\"variant\"]\n",
    "variantdf.to_csv(\"/nfs/amishra/Arpit/18k_rsid_converted_to15kvariant.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea12a477-5a18-44b0-a71a-c878d87ba566",
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Var id to target gene using V2G otargen pipeline this is R code snippet "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bbeae10-932f-4a51-a8a9-285d9d075273",
   "metadata": {},
   "source": [
    "# Load necessary libraries\n",
    "library(otargen)\n",
    "library(purrr)\n",
    "library(dplyr)\n",
    "library(readr)\n",
    "\n",
    "\n",
    "\n",
    "# Extract the variant IDs into a vector\n",
    "variant_ids <- df$variant\n",
    "\n",
    "# Define a function to retrieve gene information for a variant ID\n",
    "get_gene_info <- function(variant_id) {\n",
    "  genesForVariant(variant_id = variant_id)\n",
    "}\n",
    "\n",
    "# Wrap the function with purrr::safely to handle errors gracefully\n",
    "safe_get_gene_info <- safely(get_gene_info)\n",
    "\n",
    "# Apply the safe function to the list of variant IDs\n",
    "results <- map(variant_ids, safe_get_gene_info)\n",
    "\n",
    "# Extract the results (successful calls)\n",
    "successful_results <- map(results, \"result\")\n",
    "\n",
    "# Optionally, filter out the NULL results (where the function failed)\n",
    "successful_results_filtered <- compact(successful_results)\n",
    "\n",
    "\n",
    "# Combine results into a data frame if needed\n",
    "final_results <- bind_rows(successful_results_filtered)\n",
    "\n",
    "# Print or save the results\n",
    "print(final_results)\n",
    "saveRDS(successful_results_filtered, file = \"/nfs/amishra/Arpit/15kv2g.rds\")\n",
    "flat_results <- unlist(successful_results_filtered, recursive = FALSE)\n",
    "results_df <- bind_rows(flat_results)\n",
    "\n",
    "# Write results to a CSV file\n",
    "write.csv(results_df, file = \"/nfs/amishra/Arpit/MS_atac_cutandrun_diffpeak/15kmprabackground_multivariantcall4.csv\", row.names = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46d7f9e0-12e6-4dba-930d-4ccea3935d07",
   "metadata": {},
   "outputs": [],
   "source": [
    "####list of genes expressin in any t cell types from DICE to filter t cell expressing genes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3300bbba-bb3c-4bea-992f-2fdf2ce04108",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "from urllib.parse import urlparse\n",
    "\n",
    "# List of file URLs\n",
    "file_urls = [\n",
    "    \"https://dice-database.org/download/CD4_NAIVE_TPM.csv\",\n",
    "    \"https://dice-database.org/download/CD8_NAIVE_TPM.csv\",\n",
    "    \"https://dice-database.org/download/CD8_STIM_TPM.csv\",\n",
    "    \"https://dice-database.org/download/CD4_STIM_TPM.csv\",\n",
    "    \"https://dice-database.org/download/TH2_TPM.csv\",\n",
    "    \"https://dice-database.org/download/TH17_TPM.csv\",\n",
    "    \"https://dice-database.org/download/THSTAR_TPM.csv\",\n",
    "    \"https://dice-database.org/download/TREG_MEM_TPM.csv\",\n",
    "    \"https://dice-database.org/download/TREG_NAIVE_TPM.csv\",\n",
    "    \"https://dice-database.org/download/TFH_TPM.csv\"]\n",
    "# Initialize an empty list to store DataFrame for each file\n",
    "data_frames = []\n",
    "\n",
    "for url in file_urls:\n",
    "    # Extract file name from URL\n",
    "    parsed_url = urlparse(url)\n",
    "    file_name = os.path.basename(parsed_url.path)\n",
    "    \n",
    "    # Download or read the CSV file into a DataFrame\n",
    "    # For illustration, assuming already downloaded files are being read\n",
    "    df = pd.read_csv(file_name, index_col=None)\n",
    "    \n",
    "    # Process the DataFrame as per the example\n",
    "    df[[\"Genename\", \"Type\"]] = df.Additional_annotations.str.split(\";\", expand=True)\n",
    "    dfTPM = df.iloc[:, 3:-1]  # Adjust columns selection based on your data\n",
    "    dfTPM[\"MedianTPM\"] = dfTPM.median(numeric_only=True, axis=1)\n",
    "    df_medianTPM = dfTPM[[\"Genename\", \"MedianTPM\"]]\n",
    "    df_medianTPM_grtoreq1 = df_medianTPM[df_medianTPM[\"MedianTPM\"] >= 1]\n",
    "    df_medianTPM_grtoreq1.columns = ['geneSymbol', 'MedianTPM']\n",
    "    \n",
    "    # Add a column indicating the source DataFrame\n",
    "    df_medianTPM_grtoreq1['Source'] = file_name  # Or any identifier you prefer\n",
    "    \n",
    "    # Append to the list of data frames\n",
    "    data_frames.append(df_medianTPM_grtoreq1)\n",
    "\n",
    "# Merge all the data frames into a single DataFrame\n",
    "merged_df = pd.concat(data_frames, ignore_index=True)\n",
    "\n",
    "# Optionally, sort by geneSymbol or any other column if needed\n",
    "merged_df = merged_df.sort_values(by='geneSymbol')\n",
    "\n",
    "# Print the merged DataFrame\n",
    "print(merged_df)\n",
    "\n",
    "# Write the merged DataFrame to a CSV file\n",
    "merged_df.to_csv(\"/nfs/amishra/Arpit/Tcell_DICEpopulation_merged_medianTPM_data.csv\", index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0013fd14-565e-4314-9bc6-2a858ce1d187",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### polars to filter v2g gene list based on relevant cell type expression \n",
    "\n",
    "import polars as pl\n",
    "\n",
    "\n",
    "# Load the CSV file with Polars\n",
    "df = pl.read_csv(\"/nfs/amishra/Arpit/MS_atac_cutandrun_diffpeak/15kmprabackground_multivariantcall4.csv\")\n",
    "\n",
    "# Filter out rows where 'typeId' is not null\n",
    "dfx = df.filter(pl.col('typeId').is_not_null())\n",
    "\n",
    "# Select relevant columns\n",
    "dfy = dfx.select(['gene.symbol', 'variant', 'tissues_distance', 'sourceId', 'tissues_name', 'typeId'])\n",
    "\n",
    "# Create the 'Combined' column\n",
    "dfy = dfy.with_column((pl.col('variant').cast(pl.Utf8) + \"_\" + pl.col('gene.symbol').cast(pl.Utf8)).alias('Combined'))\n",
    "\n",
    "\n",
    "# List of tissues names to match\n",
    "tissue_names = [\n",
    "    \"Blood (eQTLGen)\", \"Blood (GTEx v8)\", \"Blood (Lepik 2017)\", \"Blood (TwinsUK)\", \"Blood plasma (Sun2018)\",\n",
    "    \"CD4+ T cell (Blueprint)\", \"CD4+ T cell (CEDAR)\", \"CD4+ T cell (Kasela 2017)\", \"CD4+ T cell (Schmiedel 2018)\",\n",
    "    \"CD4+ T cell anti-CD3-CD28 4h (Schmiedel 2018)\", \"CD8+ T cell (CEDAR)\", \"CD8+ T cell (Kasela 2017)\",\n",
    "    \"CD8+ T cell (Schmiedel 2018)\", \"CD8+ T cell anti-CD3-CD28 4h (Schmiedel 2018)\", \"Foetal thymus\",\n",
    "    \"Gtex-sqtl-whole blood\", \"Naive CD4\", \"Naive CD8\", \"T cell (GENCORD)\", \"Tfh cell (Schmiedel 2018)\",\n",
    "    \"Th1 cell (Schmiedel 2018)\", \"Th17 cell (Schmiedel 2018)\", \"Th2 cell (Schmiedel 2018)\", \"Thymus\",\n",
    "    \"Total CD4 Activated\", \"Total CD4 MF\", \"Total CD4 NonActivated\", \"Total CD8\", \"Treg (Schmiedel 2018)\",\n",
    "    \"Treg memory (Schmiedel 2018)\"\n",
    "]\n",
    "\n",
    "# Type IDs to filter\n",
    "type_ids = [\"eqtl\", \"pqtl\", \"sqtl\", \"pchic\"]\n",
    "\n",
    "# Filter rows\n",
    "filtered_df = dfy.filter(\n",
    "    ((pl.col(\"typeId\").is_in(type_ids)) & (pl.col(\"tissues_name\").is_in(tissue_names))) |\n",
    "    (~pl.col(\"typeId\").is_in(type_ids))\n",
    ")### initially thought of filtering here but it will throw away the distance v2g because they dont have tissue name thus took both filtered and non filtered and filtered later in pandas to remove NAs \n",
    "\n",
    "# Display the filtered DataFrame \n",
    "print(filtered_df)\n",
    "filtered_df.write_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/15kvariant_v2goutput.csv\")\n",
    "unique_genes = filtered_df.select(pl.col(\"gene.symbol\").unique())\n",
    "unique_genes.write_csv(\"/nfs/amishra/Arpit/raylab15kvariant_uniquegenesbasedontcellbloodfilter.csv\")\n",
    "\n",
    "### by overlapping v2g output with T cell expressed genes we generate master v2g list. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dbb1737-32b5-4ba3-b282-278041a9a713",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### using master v2g output difference background and foreground creation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dea2fe2-6388-412b-a815-0b06adfa4f54",
   "metadata": {},
   "outputs": [],
   "source": [
    "##### run this chunk with genopyc kernel\n",
    "import genopyc as gp \n",
    "import pandas as pd\n",
    "### prepare v2g df from v2g pipeline output, and prepare var_emvar_dhs etc input from Ray lab file \"240626_TGWAS_paper_rsids_per_screen_results.xlsx\"\"######\n",
    "v2g15k = pd.read_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/15kvariant_v2goutput.csv\")\n",
    "v2g15k= v2g15k[v2g15k['typeId'].notna()].copy()#### remove NAs \n",
    "dhs_emvars = pd.read_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/sanity_checkrun07282024/240626_TGWAS_paper_rsids_per_screen_results_JR (1).csv\")\n",
    "##### prepare expression filter using DICE database \n",
    "dice_exp = pd.read_csv(\"/nfs/amishra/Arpit/Tcell_DICEpopulation_merged_medianTPM_data.csv\")\n",
    "dice_exp.columns = ['gene.symbol', 'MedianTPM', 'Source']\n",
    "\n",
    "##### prepare proliferation network background (proliferationbg) using emvars_dhs_in_sc using all 56 scCRISPR v2g hits\n",
    "dhs_emvars_56 = dhs_emvars[[\"emvars_dhs_in_sc_screens\"]]\n",
    "dhs_emvars_56 = dhs_emvars_56.dropna()\n",
    "dhs_emvars_56_varid = gp.variantId_mapping(dhs_emvars_56.emvars_dhs_in_sc_screens, source=\"rsid\", target=\"variantid\")\n",
    "dhs_emvars_56_varid=list(filter(None,dhs_emvars_56_varid))\n",
    "dhs_emvars_56_varid_df = pd.DataFrame(dhs_emvars_56_varid)\n",
    "dhs_emvars_56_varid_df.columns = [\"variant\"]\n",
    "dhs_emvars_56_varid_v2g = v2g15k.merge(dhs_emvars_56_varid_df, how=\"inner\", on=\"variant\")\n",
    "dhs_emvars_56_varid_v2g_unique_genes = dhs_emvars_56_varid_v2g[\"gene.symbol\"].unique()\n",
    "dhs_emvars_56_varid_v2g_unique_genes_df = pd.DataFrame(dhs_emvars_56_varid_v2g_unique_genes, columns=[\"gene.symbol\"])\n",
    "proliferation_bg_dice = dice_exp.merge(dhs_emvars_56_varid_v2g_unique_genes_df, how=\"inner\", on=\"gene.symbol\")\n",
    "proliferation_bg_dice_uniquegenes = proliferation_bg_dice[\"gene.symbol\"].unique()\n",
    "proliferation_bg_dice_uniquegenes_df = pd.DataFrame(proliferation_bg_dice_uniquegenes, columns=[\"gene.symbol\"])\n",
    "proliferation_bg_dice_uniquegenes_df.to_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/sanity_checkrun07282024/sanitycheck_newlistJR_1/proliferation_hit_dice_uniquegenesbg.csv\")\n",
    "\n",
    "##### prepare proliferation network foreground (proliferationfg) using \"emVars+DHS+Prolif (FDR<0.1)\" for only prol hits \n",
    "dhs_emvars_prol17 = dhs_emvars[[\"emVars+DHS+Prolif (FDR<0.1)\"]]\n",
    "dhs_emvars_prol17 = dhs_emvars_prol17.dropna()\n",
    "dhs_emvars_prol17_varid = gp.variantId_mapping(dhs_emvars_prol17[\"emVars+DHS+Prolif (FDR<0.1)\"], source='rsid',target = 'variantid')\n",
    "dhs_emvars_prol17_varid_df = pd.DataFrame(dhs_emvars_prol17_varid)\n",
    "dhs_emvars_prol17_varid_df.columns = [\"variant\"]\n",
    "dhs_emvars_prol17_varid_v2g = v2g15k.merge(dhs_emvars_prol17_varid_df, how=\"inner\", on=\"variant\")\n",
    "dhs_emvars_prol17_uniquegenes = dhs_emvars_prol17_varid_v2g[\"gene.symbol\"].unique()\n",
    "dhs_emvars_prol17_unique_genes_df = pd.DataFrame(dhs_emvars_prol17_uniquegenes, columns=[\"gene.symbol\"])\n",
    "#dhs_emvars_prol17_unique_genes_df.to_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/dhs_emvars_prol17_filtered_v2g.csv\")\n",
    "proliferation_fg_dice = dice_exp.merge(dhs_emvars_prol17_unique_genes_df, how=\"inner\", on=\"gene.symbol\")\n",
    "proliferation_fg_dice_uniquegenes = proliferation_fg_dice[\"gene.symbol\"].unique()\n",
    "proliferation_fg_dice_uniquegenes_df = pd.DataFrame(proliferation_fg_dice_uniquegenes, columns=[\"gene.symbol\"])\n",
    "proliferation_fg_dice_uniquegenes_df.to_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/sanity_checkrun07282024/sanitycheck_newlistJR_1/proliferation_hit_dice_uniquegenesfg.csv\")\n",
    "\n",
    "##### prepare big network background(bignetworkbg) using all the variants in DHS v2g output \n",
    "dhs_var = dhs_emvars[[\"all_t_cell_DHS_var\"]]\n",
    "dhs_var_varid = gp.variantId_mapping(dhs_var[\"all_t_cell_DHS_var\"], source='rsid',target = 'variantid')\n",
    "dhs_var_varid_df = pd.DataFrame(dhs_var_varid)\n",
    "dhs_var_varid_df.columns = [\"variant\"]\n",
    "dhs_var_varid_df_v2g = v2g15k.merge(dhs_var_varid_df, how=\"inner\", on=\"variant\")\n",
    "dhs_var_varid_df_v2g_uniquegenes = dhs_var_varid_df_v2g[\"gene.symbol\"].unique()\n",
    "dhs_var_varid_df_v2g_uniquegenes_df = pd.DataFrame(dhs_var_varid_df_v2g_uniquegenes, columns=[\"gene.symbol\"])\n",
    "#dhs_var_varid_df_v2g_uniquegenes_df.to_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/tcellsdhs_vars_filtered_v2g.csv\")\n",
    "var_inDHS_andDICE = dice_exp.merge(dhs_var_varid_df_v2g_uniquegenes_df, how=\"inner\", on=\"gene.symbol\")\n",
    "var_inDHS_andDICEv2g_uniquegenes = var_inDHS_andDICE[\"gene.symbol\"].unique()\n",
    "var_inDHS_andDICEv2g_uniquegenes_df = pd.DataFrame(var_inDHS_andDICEv2g_uniquegenes, columns=[\"gene.symbol\"])\n",
    "var_inDHS_andDICEv2g_uniquegenes_df.to_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/sanity_checkrun07282024/sanitycheck_newlistJR_1/var_inDHSandDICE_uniquegenes_bignetworkbg.csv\")\n",
    "\n",
    "###### prepare big network foreground (bignetworkfg) using emvars in DHS v2g output \n",
    "emvar_in_DHS = dhs_emvars[[\"emVars+DHS\"]]\n",
    "emvar_in_DHS = emvar_in_DHS.dropna()\n",
    "emvar_in_DHSvarid = gp.variantId_mapping(emvar_in_DHS[\"emVars+DHS\"], source='rsid',target = 'variantid')\n",
    "emvar_in_DHSvarid_df = pd.DataFrame(emvar_in_DHSvarid)\n",
    "emvar_in_DHSvarid_df.columns = [\"variant\"]\n",
    "emvar_in_DHSvarid_df_v2g = v2g15k.merge(emvar_in_DHSvarid_df, how=\"inner\", on=\"variant\")\n",
    "emvar_in_DHSvarid_df_v2g_uniquegenes = emvar_in_DHSvarid_df_v2g[\"gene.symbol\"].unique()\n",
    "emvar_in_DHSvarid_df_v2g_uniquegenes_df = pd.DataFrame(emvar_in_DHSvarid_df_v2g_uniquegenes, columns=[\"gene.symbol\"])\n",
    "emvar_inDHS_andDICE = dice_exp.merge(emvar_in_DHSvarid_df_v2g_uniquegenes_df, how=\"inner\",on=\"gene.symbol\")\n",
    "emvar_inDHS_andDICE_uniquegenes = emvar_inDHS_andDICE[\"gene.symbol\"].unique()\n",
    "emvar_inDHS_andDICE_uniquegenes_df = pd.DataFrame(emvar_inDHS_andDICE_uniquegenes, columns=[\"gene.symbol\"])\n",
    "emvar_inDHS_andDICE_uniquegenes_df.to_csv(\"/nfs/amishra/Arpit/Raylab_Network_analysis/sanity_checkrun07282024/sanitycheck_newlistJR_1/emvar_inDHSandDICE_uniquegenes_bignetworkfg.csv\")\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base Environment",
   "language": "python",
   "name": "base"
  },
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}