phm/0.phm_data_preparation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specify input/output directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify LOCAL PATHS\n",
    "data_path = \"../test_data\"\n",
    "output_path = \"../1.mapped_CMs/phm\"\n",
    "\n",
    "# Create output path if it does not exist\n",
    "if not os.path.exists(output_path):\n",
    "    os.makedirs(output_path)\n",
    "\n",
    "# Specify names of datasets to analyze\n",
    "# PHM input/output file will be stored there\n",
    "dataset = \"test_data\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create dictionary with sample IDs per dataset\n",
    "vcf_samples_dict = {\n",
    "    dataset: pd.read_csv(\n",
    "        os.path.join(data_path, \"LCL_genotypes_chr22_samples.txt\"),\n",
    "        sep=\"\\t\",\n",
    "        header=None,\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create necessary input files for PHM\n",
    "\n",
    "1. BED file with peak coordinates\n",
    "2. Binary file of normalized counts\n",
    "\n",
    "\\+ Store common samples between ChIP-seq and genotype data to ensure proper ordering of samples across data modalities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in count matrices\n",
    "# CHANGE PATHS\n",
    "count_mtx_k4me1 = pd.read_csv(os.path.join(data_path, \"H3K4me1_chr22.bed\"), sep=\"\\t\")\n",
    "count_mtx_k4me1[\"mark\"] = \"H3K4me1\"\n",
    "\n",
    "count_mtx_k27ac = pd.read_csv(os.path.join(data_path, \"H3K27ac_chr22.bed\"), sep=\"\\t\")\n",
    "count_mtx_k27ac[\"mark\"] = \"H3K27ac\"\n",
    "\n",
    "all_peaks_mtx = pd.concat([count_mtx_k4me1, count_mtx_k27ac])\n",
    "all_peaks_mtx = all_peaks_mtx.sort_values([\"#Chr\", \"start\"])\n",
    "all_peaks_mtx = all_peaks_mtx.rename(columns={\"#Chr\": \"#chr\"})\n",
    "all_peaks_mtx.loc[:, \"#chr\"] = all_peaks_mtx.loc[:, \"#chr\"].astype(str)\n",
    "\n",
    "all_peaks_mtx[\"Peak\"] = np.arange(1, all_peaks_mtx.shape[0] + 1)\n",
    "\n",
    "all_peaks_mtx = all_peaks_mtx.astype({\"start\": int, \"end\": int})\n",
    "all_peaks_mtx = all_peaks_mtx.dropna(axis=\"columns\")\n",
    "\n",
    "\n",
    "chromosomes = list(\n",
    "    set(all_peaks_mtx.loc[:, \"#chr\"]) - set([\"chrX\", \"chrY\", \"chrM\", \"X\", \"Y\", \"M\"])\n",
    ")\n",
    "if not all([isinstance(el, str) for el in chromosomes]):\n",
    "    sys.exit(\"Chromosome IDs have different types! Standardize the format\")\n",
    "\n",
    "if \"chr\" in str(chromosomes[0]):\n",
    "    suffix = \"\"\n",
    "else:\n",
    "    suffix = \"chr\"\n",
    "\n",
    "# Get common samples between ChIP-seq and genotype data\n",
    "final_samples = sorted(\n",
    "    list(\n",
    "        set(all_peaks_mtx.columns).intersection(\n",
    "            set(vcf_samples_dict.get(dataset).iloc[:, 0])\n",
    "        )\n",
    "    )\n",
    ")\n",
    "\n",
    "# Store only common samples between ChIP-seq and genotype data\n",
    "pd.DataFrame(final_samples).to_csv(\n",
    "    os.path.join(output_path, dataset + \"_LCL_chr22_samples.txt\"),\n",
    "    sep=\"\\t\",\n",
    "    index=False,\n",
    "    header=False,\n",
    ")\n",
    "\n",
    "first_cols = [\"#chr\", \"start\", \"end\", \"pid\", \"gid\", \"strand\"]\n",
    "last_cols = [\"mark\", \"Peak\"]\n",
    "ordered_cols = first_cols + final_samples + last_cols\n",
    "columns_for_tsv = [\"Peak\"] + final_samples\n",
    "\n",
    "count_matrix = all_peaks_mtx.loc[:, ordered_cols].copy()\n",
    "\n",
    "count_matrices_by_chromosome_dict = {\n",
    "    suffix + str(chromosome): chromosome_count_matrix.sort_values(\"start\")\n",
    "    for chromosome, chromosome_count_matrix in count_matrix.groupby(\"#chr\")\n",
    "    if chromosome in chromosomes\n",
    "}\n",
    "\n",
    "for chromosome, chromosome_count_matrix in count_matrices_by_chromosome_dict.items():\n",
    "    if (chromosome != \"chrX\") and (chromosome != \"chrY\") and (chromosome != \"chrM\"):\n",
    "        path = os.path.join(output_path, str(chromosome))\n",
    "        if not os.path.exists(os.path.join(path, \"hm_output\")) and not os.path.exists(\n",
    "            os.path.join(path, \"phm_output\")\n",
    "        ):\n",
    "            os.makedirs(os.path.join(path, \"hm_output\"))\n",
    "            os.makedirs(os.path.join(path, \"phm_output\"))\n",
    "        chromosome_count_matrix[\"Peak\"] = np.arange(\n",
    "            1, chromosome_count_matrix.shape[0] + 1\n",
    "        )\n",
    "        chromosome_count_matrix_for_tsv = chromosome_count_matrix[columns_for_tsv]\n",
    "\n",
    "        # Create .bed file with peak coordinates\n",
    "        if type(chromosome_count_matrix[\"#chr\"].iloc[0]) != int:\n",
    "            chromosome_count_matrix[\"#chr\"] = chromosome_count_matrix[\n",
    "                \"#chr\"\n",
    "            ].str.replace(\"chr\", \"\")\n",
    "            chromosome_count_matrix[\"#chr\"] = chromosome_count_matrix[\"#chr\"].astype(\n",
    "                int\n",
    "            )\n",
    "        chromosome_count_matrix[[\"#chr\", \"start\", \"end\", \"mark\"]].to_csv(\n",
    "            os.path.join(path, chromosome + \"_peak_coordinates.bed\"),\n",
    "            sep=\"\\t\",\n",
    "            index=False,\n",
    "            header=False,\n",
    "        )\n",
    "\n",
    "        # Create a binary file from normalized counts\n",
    "        records_array = chromosome_count_matrix[final_samples].to_records(index=False)\n",
    "        records_array.tofile(\n",
    "            os.path.join(path, \"normalized_counts_\" + chromosome + \".bin\")\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IMPORTANT:\n",
    "\n",
    "1. Make sure the samples in a VCF file are orderded in the same way as samples in count matrices.\n",
    "2. Pay attention to the chromosome formats in a VCF file and count matrices, it should be the same, and formated as a string without \"chr\" prefix, i.e., chr22 -> 22.\n",
    "3. Make sure to change paths prior to executing the cell!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "dataset=\"test_data\";\n",
    "data_path=\"../test_data\";\n",
    "output_path=\"../1.mapped_CMs/phm\";\n",
    "\n",
    "## Subset samples from the VCF file\n",
    "bcftools view \\\n",
    "    -S ${output_path}/${dataset}_LCL_chr22_samples.txt \\\n",
    "    ${data_path}/LCL_genotypes_chr22.vcf.gz \\\n",
    "    > ${output_path}/${dataset}_sample_intersection.vcf\n",
    "\n",
    "## zip and index VCF\n",
    "bgzip ${output_path}/${dataset}_sample_intersection.vcf\n",
    "tabix -p vcf ${output_path}/${dataset}_sample_intersection.vcf.gz\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compresses and index files per chromosome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/pushkare/Chromatin_modules/1.mapped_CMs/phm\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "output_path=\"../1.mapped_CMs/phm\";\n",
    "cd ${output_path};\n",
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "dataset=\"test_data\";\n",
    "output_path=\"../1.mapped_CMs/phm\";\n",
    "cd ${output_path};\n",
    "input_vcf=${dataset}_sample_intersection.vcf.gz\n",
    "\n",
    "for chromosome_dir in */\n",
    "do  \n",
    "    ## Subset a chromosome from VCF file\n",
    "    tabix \\\n",
    "        ${input_vcf} \\\n",
    "        ${chromosome_dir%/} \\\n",
    "        > ${chromosome_dir%/}/${dataset}_${chromosome_dir%/}.vcf\n",
    "\n",
    "    ## If chromosome IDs in the VCF file contain \"chr\" character,\n",
    "    ## replace it with an empty string, s.t. \"chr22\" -> \"22\"\n",
    "    awk '{gsub(/^chr/,\"\"); print}' ${chromosome_dir%/}/${dataset}_${chromosome_dir%/}.vcf \\\n",
    "        > ${chromosome_dir%/}/${dataset}_${chromosome_dir%/}_no_chr.vcf\n",
    "\n",
    "\n",
    "    ## Compress and index BED, VCF files\n",
    "    for filename in ${chromosome_dir%/}/*.bed\n",
    "    do\n",
    "        bgzip $filename\n",
    "        tabix -p bed $filename.gz\n",
    "    done\n",
    "\n",
    "    for filename in ${chromosome_dir%/}/*.vcf\n",
    "    do\n",
    "        bgzip $filename\n",
    "        tabix -p vcf $filename.gz\n",
    "    done\n",
    "done\n",
    "\n",
    "\n",
    "## # If chromosome IDs in the VCF file *do not* contain \"chr\",\n",
    "## # replace lines 10-18 with:\n",
    "## chr_id=${chromosome_dir//\"chr\"/}\n",
    "## tabix \\\n",
    "##     ${input_vcf} \\\n",
    "##     ${chr_id%/} \\\n",
    "##     > ${chromosome_dir%/}/${dataset}_${chromosome_dir%/}_no_chr.vcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DONE\n"
     ]
    }
   ],
   "source": [
    "print(\"DONE\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}