{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the following utilities\n",
    "- http://hgdownload.soe.ucsc.edu/admin/exe/macOSX.x86_64/bigWigToBedGraph\n",
    "- http://hgdownload.soe.ucsc.edu/admin/exe/macOSX.x86_64/hgWiggle\n",
    "- bedtools\n",
    "- pandas (Python package)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the following annotations\n",
    "- https://public.hoffman2.idre.ucla.edu/ernst/R0RG6/LECIF/mm10.LECIFv1.1.bw\n",
    "- http://hgdownload.cse.ucsc.edu/goldenpath/mm10/phastCons60way/mm10.60way.phastCons.bw\n",
    "- http://hgdownload.cse.ucsc.edu/goldenpath/mm10/phastCons60way/mm10.60way.phastCons60wayPlacental.bw\n",
    "- http://hgdownload.cse.ucsc.edu/goldenpath/mm10/phyloP60way/mm10.60way.phyloP60way.bw\n",
    "- http://hgdownload.cse.ucsc.edu/goldenpath/mm10/phyloP60way/mm10.60way.phyloP60wayPlacental.bw\n",
    "- https://drive.google.com/file/d/1TbqDipE2zcw2-FEs5aVWKkQ8zYwjOsQq/view?usp=sharing\n",
    "- https://github.com/shorvath/MammalianMethylationConsortium/blob/main/Annotations%2C%20Amin%20Haghani/Mammals/Mus_musculus.grcm38.100.HorvathMammalMethylChip40.v1.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write a bed file with all mm10 coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "input_filename = 'Mus_musculus.grcm38.100.HorvathMammalMethylChip40.v1.txt'\n",
    "df = pd.read_csv(input_filename,usecols=['CGid','seqnames','CGstart','CGend']).dropna()\n",
    "df = df.sort_values(by=['seqnames','CGstart'])\n",
    "df['chrom'] = \"chr\"+df['seqnames'] \n",
    "df['CGstart'] = df['CGstart'].astype(int)\n",
    "df['CGend'] = df['CGend'].astype(int)\n",
    "df[['chrom','CGstart','CGend','CGid']].to_csv('mm10_coord.bed',sep='\\t',header=False,index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert bigWig files to gzipped bed files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "./bigWigToBedGraph mm10.LECIFv1.1.bw mm10.LECIF.bed\n",
    "gzip mm10.LECIF.bed\n",
    "\n",
    "./bigWigToBedGraph mm10.60way.phastCons.bw mm10.phastCons60way.bed\n",
    "gzip mm10.phastCons60way.bed\n",
    "\n",
    "./bigWigToBedGraph mm10.60way.phastCons60wayPlacental.bw mm10.phastCons60wayPlacental.bed\n",
    "gzip mm10.phastCons60wayPlacental.bed\n",
    "\n",
    "./bigWigToBedGraph mm10.60way.phyloP60way.bw mm10.phyloP60way.bed\n",
    "gzip mm10.phyloP60way.bed\n",
    "\n",
    "./bigWigToBedGraph mm10.60way.phyloP60wayPlacental.bw mm10.phyloP60wayPlacental.bed\n",
    "gzip mm10.phyloP60wayPlacental.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map LECIF score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! bedtools map -a mm10_coord.bed -b mm10.LECIF.bed.gz -c 4 -o mean > mm10_coord.LECIF.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map sequence constraint scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "for f in phastCons60way phyloP60way phastCons60wayPlacental phyloP60wayPlacental; do\n",
    "    bedtools map -a mm10_coord.bed -b mm10.$f.bed.gz -c 4 -o mean > mm10_coord.$f.bed\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map ConsHMM states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! bedtools map -a mm10_coord.bed -b 60_states_segmentation.bed.gz -c 4 -o first > mm10_coord.ConsHMM.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map ENCODE ChromHMM states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "path_prefix=http://hgdownload.soe.ucsc.edu/gbdb/mm10/encode3/chromHmm/\n",
    "output_filename=mm10_coord.ChromHMM.bed\n",
    "\n",
    "cp mm10_coord.bed $output_filename\n",
    "for i in `cat mm10_ChromHMM_filenames.txt`; do\n",
    "    ./bigBedToBed $path_prefix/$i stdout | bedtools map -a $output_filename -b - -c 4 -o first > tmp.bed\n",
    "    mv tmp.bed $output_filename\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Combine them into a formatted csv file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "coord_filename = 'mm10_coord.bed'\n",
    "bed_filenames = ['mm10_coord.LECIF.bed',\n",
    "                 'mm10_coord.phastCons60way.bed',\n",
    "                 'mm10_coord.phastCons60wayPlacental.bed',\n",
    "                 'mm10_coord.phyloP60way.bed',\n",
    "                 'mm10_coord.phyloP60wayPlacental.bed',\n",
    "                 'mm10_coord.ConsHMM.bed',\n",
    "                 'mm10_coord.ChromHMM.bed']\n",
    "\n",
    "position_col_names = ['chrom','start','end','probeID']\n",
    "\n",
    "chromhmm_col_name_prefix = 'ChromHMMState_ENCODE_GorkinEtAl2017bioRxiv_'\n",
    "chromhmm_col_names = [chromhmm_col_name_prefix+i.strip().replace('.bb','').replace('encode3RenChromHmm','') for i in open('mm10_ChromHMM_filenames.txt')]\n",
    "\n",
    "names = [position_col_names + n for n in [['LECIFScore_KwonErnst2012NatureCommunications'],\n",
    "                                          ['PhastConsScore_SiepelEtAl2005GenomeResearch'],\n",
    "                                          ['PhyloPScore_PollardEtAl2009GenomeResearch'],\n",
    "                                          ['PhastConsPlacentalScore_SiepelEtAl2005GenomeResearch'],\n",
    "                                          ['PhyloPPlacentalScore_PollardEtAl2009GenomeResearch'],\n",
    "                                          ['ConsHMMState_Multiz60way_ArnesonEtAl2020NARGenomBioinform'],\n",
    "                                          chromhmm_col_names]]\n",
    "\n",
    "df = pd.read_table(coord_filename,header=None,names=position_col_names)\n",
    "for i in range(len(bed_filenames)):\n",
    "    input_df = pd.read_table(bed_filenames[i],header=None,names=names[i])\n",
    "    df = df.merge(input_df,on=position_col_names,how='left')\n",
    "    \n",
    "df.to_csv('mm10_annotErnstLab.csv',index=False)"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}