notebooks/multispecies_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-species Hi-C data downloading\n",
    "\n",
    "This notebook illustrates how to download the multi-species Hi-C file from DNAZoo data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Get the url list for downloading the hic files from dnazoo\n",
    "import os\n",
    "import sys\n",
    "import json \n",
    "import requests\n",
    " \n",
    "def download_file(url, filename):\n",
    "    response = requests.get(url)\n",
    "    if response.status_code == 200:\n",
    "        with open(filename, 'wb') as f:\n",
    "            f.write(response.content)\n",
    "        print(f\"Finish download {filename}\")\n",
    "    else:\n",
    "        print(f\"Request fail: {response.status_code}\")\n",
    "\n",
    "\n",
    "input_txt = \"../data/species_list.txt\"\n",
    "output_dir = \"../data/dnazoo_json/\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "species_list = []\n",
    "with open(input_txt) as f:\n",
    "    for line in f:\n",
    "        species_list.append(line.strip(\"\\n\"))\n",
    "#https://dnazoo.s3.wasabisys.com/Acyrthosiphon_kondoi/README.json\n",
    "final_url_list = []\n",
    "root_url = \"https://dnazoo.s3.wasabisys.com/\"\n",
    "miss_list=[]\n",
    "for species in species_list:\n",
    "    if \"/\" not in species:\n",
    "        species += \"/\"\n",
    "    readme_url = root_url + species + \"README.json\"\n",
    "\n",
    "    readme_local = os.path.join(output_dir, species.replace(\"/\",\"\") + \"_readme.json\")\n",
    "    if not os.path.exists(readme_local) or os.path.getsize(readme_local) <= 10:\n",
    "        #os.system(f\"wget {readme_url} -O {readme_local}\")\n",
    "        download_file(readme_url, readme_local)\n",
    "    if not os.path.exists(readme_local) or os.path.getsize(readme_local) <= 10:\n",
    "        print(f\"Error: {readme_local} does not exist or is empty\")\n",
    "        miss_list.append(species)\n",
    "        continue\n",
    "    json_info = json.load(open(readme_local,\"r\"))\n",
    "    hic_name = json_info['links']['cMap']['name']\n",
    "    hic_url = root_url + species + hic_name\n",
    "    if \"hicUrl\" in json_info['links']['cMap']:\n",
    "        \n",
    "        dropbox_url = json_info['links']['cMap'][\"hicUrl\"]\n",
    "        dropbox_url = dropbox_url.split(\"?\")[0]\n",
    "    else:\n",
    "        dropbox_url = hic_url\n",
    "    #final_url_list.append(hic_url)\n",
    "    new_dict={species:[hic_url, dropbox_url]}\n",
    "    final_url_list.append(new_dict)\n",
    "import pickle\n",
    "pickle.dump(final_url_list, open(\"../data/dnazoo_url_list.pkl\",\"wb\"))\n",
    "print(\"Done\")\n",
    "print(\"Missed species:\", miss_list)\n",
    "#clean the json file directories\n",
    "os.rmdir(output_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Download the hic files from the url list\n",
    "import os\n",
    "import sys\n",
    "import pickle\n",
    "import requests\n",
    " \n",
    "def download_file(url, filename):\n",
    "    response = requests.get(url)\n",
    "    if response.status_code == 200:\n",
    "        with open(filename, 'wb') as f:\n",
    "            f.write(response.content)\n",
    "        print(f\"Finish download {filename}\")\n",
    "    else:\n",
    "        print(f\"Request fail: {response.status_code}\")\n",
    "input_pkl = \"../data/dnazoo_url_list.pkl\"\n",
    "output_dir =\"../data/multi_species_hic/\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "data = pickle.load(open(input_pkl, 'rb'))\n",
    "def download_url(url,file_path):\n",
    "    os.system(f\"wget {url} -O {file_path}\")\n",
    "\n",
    "from multiprocessing import Pool\n",
    "#you can consider use the above choice to download multiple Hi-C simultaneously to accelerate downloading speed.\n",
    "for species in data:\n",
    "    url_list = data[species]\n",
    "    file_path = os.path.join(output_dir, f\"{species}.hic\")\n",
    "    if isinstance(url_list, list):\n",
    "        for url in url_list:\n",
    "            download_url(url, file_path)\n",
    "            if os.path.exists(file_path) and os.path.getsize(file_path) > 0:\n",
    "                break\n",
    "    elif isinstance(url_list, str):\n",
    "        download_url(url_list, file_path)\n",
    "    else:\n",
    "        print(f\"Error: {species} url list is not valid\")\n",
    "print(\"Done\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Request fail: 404\n",
      "Error: ../data/dnazoo_json/Chinchilla_x_readme.json does not exist or is empty\n",
      "Request fail: 404\n",
      "Error: ../data/dnazoo_json/Preliminary_readme.json does not exist or is empty\n",
      "Request fail: 404\n",
      "Error: ../data/dnazoo_json/annotations_readme.json does not exist or is empty\n",
      "Done\n",
      "Missed species: ['Chinchilla_x/', 'Preliminary/', 'annotations/']\n"
     ]
    }
   ],
   "source": [
    "# 3. download all associated fasta files \n",
    "import os\n",
    "import sys\n",
    "import json \n",
    "import requests\n",
    " \n",
    "def download_file(url, filename):\n",
    "    response = requests.get(url)\n",
    "    if response.status_code == 200:\n",
    "        with open(filename, 'wb') as f:\n",
    "            f.write(response.content)\n",
    "        print(f\"Finish download {filename}\")\n",
    "    else:\n",
    "        print(f\"Request fail: {response.status_code}\")\n",
    " \n",
    "\n",
    "input_txt = \"../data/species_list.txt\"\n",
    "output_dir = \"../data/dnazoo_json/\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "species_list = []\n",
    "with open(input_txt) as f:\n",
    "    for line in f:\n",
    "        species_list.append(line.strip(\"\\n\"))\n",
    "#https://dnazoo.s3.wasabisys.com/Acyrthosiphon_kondoi/README.json\n",
    "final_url_list = []\n",
    "root_url = \"https://dnazoo.s3.wasabisys.com/\"\n",
    "miss_list=[]\n",
    "for species in species_list:\n",
    "    if \"/\" not in species:\n",
    "        species += \"/\"\n",
    "    readme_url = root_url + species + \"README.json\"\n",
    "\n",
    "    readme_local = os.path.join(output_dir, species.replace(\"/\",\"\") + \"_readme.json\")\n",
    "    if not os.path.exists(readme_local) or os.path.getsize(readme_local) <= 10:\n",
    "    #os.system(f\"wget {readme_url} -O {readme_local}\")\n",
    "        download_file(readme_url, readme_local)\n",
    "    if not os.path.exists(readme_local) or os.path.getsize(readme_local) <= 10:\n",
    "        print(f\"Error: {readme_local} does not exist or is empty\")\n",
    "        miss_list.append(species)\n",
    "        continue\n",
    "    json_info = json.load(open(readme_local,\"r\"))\n",
    "    hic_name = json_info['links']['cMap']['name']\n",
    "    species_name = hic_name.replace(\".rawchrom.hic\",\"\")\n",
    "    hic_url = root_url + species + species_name+\"_HiC.fasta.gz\"\n",
    "    if \"cFasta\" in json_info['links']:\n",
    "        \n",
    "        dropbox_url = json_info['links']['cFasta']\n",
    "        dropbox_url = dropbox_url.split(\"?\")[0]\n",
    "    else:\n",
    "        dropbox_url = hic_url\n",
    "    #final_url_list.append(hic_url)\n",
    "    new_dict={species:[hic_url, dropbox_url]}\n",
    "    final_url_list.append(new_dict)\n",
    "import pickle\n",
    "pickle.dump(final_url_list, open(\"../data/dnazoo_fasta_url_list.pkl\",\"wb\"))\n",
    "print(\"Done\")\n",
    "print(\"Missed species:\", miss_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4 download all fasta files\n",
    "\n",
    "import sys\n",
    "import pickle\n",
    "import requests\n",
    " \n",
    "def download_file(url, filename):\n",
    "    response = requests.get(url)\n",
    "    if response.status_code == 200:\n",
    "        with open(filename, 'wb') as f:\n",
    "            f.write(response.content)\n",
    "        print(f\"Finish download {filename}\")\n",
    "    else:\n",
    "        print(f\"Request fail: {response.status_code}\")\n",
    "input_pkl = \"../data/dnazoo_fasta_url_list.pkl\"\n",
    "output_dir =\"../data/multi_species_hic/\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "data = pickle.load(open(input_pkl, 'rb'))\n",
    "def download_url(url,file_path):\n",
    "    os.system(f\"wget {url} -O {file_path}\")\n",
    "\n",
    "from multiprocessing import Pool\n",
    "#you can consider use the above choice to download multiple Hi-C simultaneously to accelerate downloading speed.\n",
    "for species in data:\n",
    "    url_list = data[species]\n",
    "    file_path = os.path.join(output_dir, f\"{species}.hic\")\n",
    "    if isinstance(url_list, list):\n",
    "        for url in url_list:\n",
    "            download_url(url, file_path)\n",
    "            if os.path.exists(file_path) and os.path.getsize(file_path) > 0:\n",
    "                break\n",
    "    elif isinstance(url_list, str):\n",
    "        download_url(url_list, file_path)\n",
    "    else:\n",
    "        print(f\"Error: {species} url list is not valid\")\n",
    "print(\"Done\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5 Generate the chromosme size file from the fasta file\n",
    "\n",
    "# please first install samtools in your system: https://github.com/samtools/samtools‌\n",
    "import os\n",
    "import sys\n",
    "input_dir = \"../data/multi_species_hic/\"\n",
    "listfiles= [f for f in os.listdir(input_dir) if f.endswith(\".fasta\") or f.endswith(\".fasta.gz\")]\n",
    "from multiprocessing import Pool\n",
    "p=Pool(10)\n",
    "def process_file(cur_path):\n",
    "    item = os.path.basename(cur_path)\n",
    "    if item.endswith(\".gz\"):\n",
    "        fasta_path = cur_path.replace(\".gz\",\"\")\n",
    "        if os.path.exists(fasta_path):\n",
    "            os.remove(fasta_path)\n",
    "        os.system(f\"gzip -d {cur_path}\")\n",
    "        cur_path = cur_path.replace(\".gz\",\"\")\n",
    "    command_line = \"samtools faidx \" + cur_path\n",
    "    expect_output = cur_path + \".fai\"\n",
    "    if os.path.exists(expect_output) and os.path.getsize(expect_output) > 100:\n",
    "        print(f\"{cur_path} already have fai file, skip.\")\n",
    "        return\n",
    "    os.system(command_line)\n",
    "\n",
    "for item in listfiles:\n",
    "    cur_path = os.path.join(input_dir, item)\n",
    "    p.apply_async(process_file,args=(cur_path,))\n",
    "p.close()\n",
    "p.join()\n",
    "\n",
    "# This is to reorganize Hi-C data into a more organized format with data stored in each chromosome\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6. Convert the hic files to .pkl format\n",
    "# This is based on the utils/hic2array.py script\n",
    "import os\n",
    "import sys\n",
    "input_dir = \"../data/multi_species_hic/\"\n",
    "listfiles= [f for f in os.listdir(input_dir) if f.endswith(\".hic\")]\n",
    "resolution = 10000\n",
    "script_path= \"../utils/hic2array.py\"\n",
    "#connvet script\n",
    "for item in listfiles:\n",
    "    cur_hic = os.path.join(input_dir, item)\n",
    "    output_path = cur_hic.replace(\".hic\",\".pkl\")\n",
    "    if os.path.exists(output_path) and os.path.getsize(output_path) > 100:\n",
    "        print(f\"{output_path} already exists, skip.\")\n",
    "        continue\n",
    "    command_line= f\"python3 hic2array.py {cur_hic} {output_path} {resolution} 0 2\"\n",
    "    #no normalization and only use the cis contact\n",
    "    os.system(command_line)\n",
    "print(\"Done\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 7. Reformat the pkl to chromosome-based pkl based on the chromosisze file\n",
    "\n",
    "import os\n",
    "import pickle\n",
    "import numpy as np\n",
    "\n",
    "def filter_region(current_data, current_start, current_end):\n",
    "    from scipy.sparse import coo_matrix\n",
    "    row = current_data.row\n",
    "    col = current_data.col\n",
    "    data = current_data.data\n",
    "    select_index1 = (row>=current_start) & (row<current_end)\n",
    "    select_index2 = (col>=current_start) & (col<current_end)\n",
    "    select_index = select_index1 & select_index2\n",
    "    new_row = row[select_index] - current_start\n",
    "    new_col = col[select_index] - current_start\n",
    "    new_data = data[select_index]\n",
    "    new_array = coo_matrix((new_data, (new_row, new_col)), shape=(current_end-current_start, current_end-current_start))\n",
    "    return new_array\n",
    "def reformat_pkl(input_pkl, input_chrom_file, resolution):\n",
    "    data=pickle.load(open(input_pkl, 'rb'))\n",
    "    #parse the chrom size file to get start and end position\n",
    "    bounding_index_info={}\n",
    "    with open(input_chrom_file, 'r') as f:\n",
    "        for line in f:\n",
    "            info= line.strip().split(\"\\t\")\n",
    "            chrom_name = info[0]\n",
    "            current_size = int(np.floor(int(info[1])/resolution/2)) #the recorded fasta file is double than the actual size\n",
    "            current_start = int(int(info[2])//resolution/2)\n",
    "            bounding_index_info[chrom_name] = [current_start,current_size]\n",
    "    final_data = {}\n",
    "    for chrom in data:\n",
    "        if \"assembly\" not in chrom:\n",
    "            final_data[chrom] = data[chrom]\n",
    "            continue\n",
    "        print(f\"Start processing {chrom}\")\n",
    "        print(f\"Chrom size: {data[chrom].shape}\")\n",
    "        #read the chrom info to get the corresponding sub array\n",
    "        current_data = data[chrom]\n",
    "        current_max = current_data.shape[0]\n",
    "        for chrom_tmp in bounding_index_info:\n",
    "            current_start, current_size = bounding_index_info[chrom_tmp]\n",
    "            if current_start > current_max:\n",
    "                print(f\"Error: {chrom_tmp} size is larger than the data size.\")\n",
    "                continue\n",
    "            current_end = min(current_max, current_start+current_size)\n",
    "            actual_size = current_end - current_start\n",
    "            if actual_size<=50:#skip the small chrom\n",
    "                print(f\"Warning: {chrom_tmp} size is smaller than 50, skip.\")\n",
    "                continue\n",
    "            final_data[chrom_tmp] = filter_region(current_data, current_start, current_end)\n",
    "            print(f\"Finish processing {chrom_tmp} in {chrom}\")\n",
    "            print(f\"filter Chrom size: {final_data[chrom_tmp].shape}\")\n",
    "    return final_data\n",
    "if __name__ == '__main__':\n",
    "    import sys\n",
    "    input_pkl_dir = \"../data/multi_species_hic/\"\n",
    "    input_chrom_size_dir = \"../data/multi_species_hic/\"\n",
    "    output_pkl_dir = \"../data/multi_species_final_hic/\"\n",
    "    os.makedirs(output_pkl_dir, exist_ok=True)\n",
    "    resolution = 10000\n",
    "    listfiles= [f for f in os.listdir(input_pkl_dir) if f.endswith(\".pkl\")]\n",
    "    listfiles.sort()\n",
    "    for item in listfiles:\n",
    "        input_pkl = os.path.join(input_pkl_dir, item)\n",
    "\n",
    "        input_chrom_file = os.path.join(input_chrom_size_dir, item.replace(\".pkl\",\".fasta.fai\"))\n",
    "        output_pkl = os.path.join(output_pkl_dir, item)\n",
    "        if not os.path.exists(input_chrom_file):\n",
    "            print(f\"Error: {input_chrom_file} does not exist.\")\n",
    "            continue\n",
    "        if os.path.exists(output_pkl) or os.path.exists(output_pkl.replace(\".pkl\",\"\")):\n",
    "            print(f\"Skip: {output_pkl} already exists.\")\n",
    "            continue\n",
    "        \n",
    "        reformat_data = reformat_pkl(input_pkl, input_chrom_file, resolution)\n",
    "        pickle.dump(reformat_data, open(output_pkl, 'wb'))\n",
    "        print(f\"Finish processing {item}\")\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}