01_blast-genomes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Genome IDs\n",
    "\n",
    "For every genome, download the GenBank assembly.  Make a blast database. Run tblastn for query genomes against \n",
    "\n",
    "+ bambooshark ASM401019v1\n",
    "+ homo GRCh38.p13\n",
    "+ mus GRCm38.p6\n",
    "+ gallus GRCg6a\n",
    "+ xenopus Xenopus_tropicalis_v9.1\n",
    "+ lepisosteus LepOcu1\n",
    "+ bonytongue fSclFor1.1\n",
    "+ ictalurus IpCoco_1.2\n",
    "+ danio GRCz11\n",
    "+ esox Eluc_v4\n",
    "+ gadus gadMor3.0\n",
    "+ takifugu fTakRub1.2\n",
    "\n",
    "\n",
    "#### Approach\n",
    "\n",
    "Construct complete dataset that has:\n",
    " \n",
    " + query_species (homo/mus)\n",
    " + query_name (homo or mus gene that gave the hit)\n",
    " + query_start (position that query match started at)\n",
    " + query_end (position that query match ended at)\n",
    " + subject_species (genome that was searched)\n",
    " + block (chromosome/assembly/contig it is on, genbank numbering)\n",
    " + subject_start (position that subject started at)\n",
    " + subject_end (position that subject ended at)\n",
    " + e_value (hit e-value)\n",
    " + sequence (hit sequence)\n",
    " \n",
    " Save this to `all-hits.csv`\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set genomes to query\n",
    "subject_species = [\"bambooshark\",\n",
    "                   \"homo\",\n",
    "                   \"gallus\",\n",
    "                   \"xenopus\",\n",
    "                   \"lepisosteus\",\n",
    "                   \"bonytongue\",\n",
    "                   \"ictalurus\",\n",
    "                   \"danio\",\n",
    "                   \"esox\",\n",
    "                   \"gadus\",\n",
    "                   \"takifugu\"]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required modules\n",
    "\n",
    "import phylopandas as phy\n",
    "from phylogenetics import tools\n",
    "\n",
    "%matplotlib inline \n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from shapely.geometry.polygon import LinearRing, Polygon\n",
    "from descartes import PolygonPatch\n",
    "\n",
    "import re, pickle, glob, os, random, string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------------------------------------------------------------\n",
    "# Functions for BLASTING human and zebrafish gene lists against\n",
    "# genomes of various critters.\n",
    "# ---------------------------------------------------------------\n",
    "\n",
    "def location_from_description(description):\n",
    "    \"\"\"\n",
    "    Get the genome location by parsing the description line from an ncbi\n",
    "    genome fasta file sequence header.\n",
    "    \"\"\"\n",
    "        \n",
    "    try:\n",
    "        location = description.split(\"chromosome:\")[1].split(\" gene:\")[0].split(\":\")\n",
    "    except IndexError:\n",
    "        location = description.split(\"scaffold:\")[1].split(\" gene:\")[0].split(\":\")\n",
    "        \n",
    "    for i in range(2,5):\n",
    "        location[i] = int(location[i])\n",
    "    \n",
    "    try:\n",
    "        gene_name = description.split(\"gene_symbol:\")[1].split(\"description:\")[0]\n",
    "    except IndexError:\n",
    "        gene_name = description.split(\"gene:\")[1].split(\" \")[0]\n",
    "\n",
    "    gene_name = gene_name.strip()\n",
    "    \n",
    "    location.append(gene_name)\n",
    "    \n",
    "    return location\n",
    "\n",
    "def load_genome(fasta_file):\n",
    "    \"\"\"\n",
    "    Load a genome from an ncbi fasta file, extracting useful information\n",
    "    into columns of the data frame.\n",
    "    \"\"\"\n",
    "\n",
    "    df = phy.read_fasta(fasta_file)\n",
    "    \n",
    "    locations = [location_from_description(d) for d in df.description]\n",
    "\n",
    "    df[\"chromosome\"] = [l[1] for l in locations]\n",
    "    df[\"start\"] = [l[2] for l in locations]\n",
    "    df[\"end\"] = [l[3] for l in locations]\n",
    "    df[\"direction\"] = [l[4] for l in locations]\n",
    "    df[\"gene_name\"] = [l[5] for l in locations]\n",
    "\n",
    "    return df\n",
    "    \n",
    "def range_mask(row,min_search,max_search,chromosome):\n",
    "    \"\"\"\n",
    "    See if some part of a gene falls between min_search and max_search\n",
    "    given its start and stop positions.\n",
    "    \"\"\"\n",
    "    \n",
    "    if row['chromosome'] != str(chromosome):\n",
    "        return False\n",
    "    \n",
    "    span = [row['start'],row['end']]\n",
    "    abs_start = np.min(span)\n",
    "    abs_end = np.max(span)\n",
    "        \n",
    "    if (abs_start < min_search and abs_end >= min_search) or \\\n",
    "       (abs_start < max_search and abs_end >= max_search) or \\\n",
    "       (abs_start >= min_search and abs_end <= max_search):\n",
    "        \n",
    "        return True\n",
    "\n",
    "    return False\n",
    "\n",
    "        \n",
    "def create_region_list(df,pattern,chromosome,extend_out=3000000,take_only_longest=True):\n",
    "    \"\"\"\n",
    "    Given a genome `df`, find all genes within `extend_out` of the gene matched by\n",
    "    `pattern` regular expression and on the specified `chromosome`.  \n",
    "    For genes in the region with duplicate names, optionally take only the longest\n",
    "    gene.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Get genes that match the pattern on the appropriate chromosome\n",
    "    df_pattern = df[np.logical_and(df.description.str.contains(pattern),\n",
    "                                   df.chromosome==str(chromosome))]\n",
    "        \n",
    "    # Find the minimum and maximum extent of the gene\n",
    "    all_locations = list(df_pattern.start)\n",
    "    all_locations.extend(df_pattern.end)\n",
    "    \n",
    "    min_value = np.min(all_locations)\n",
    "    max_value = np.max(all_locations)\n",
    "\n",
    "    # Define region to search around the matched region\n",
    "    min_search = min_value - extend_out\n",
    "    max_search = max_value + extend_out\n",
    "    \n",
    "    # Find genes that overlap the range defined by min_search and max_search\n",
    "    df_near = df[df.apply(range_mask,axis=1,args=(min_search,max_search,chromosome))]\n",
    "    \n",
    "    # Make a dictionary of all genes with identical names\n",
    "    unique_set = {}\n",
    "    for i in range(len(df_near.gene_name)):\n",
    "        try:\n",
    "            unique_set[df_near.gene_name.iloc[i]][0].append(df_near.start.iloc[i])\n",
    "            unique_set[df_near.gene_name.iloc[i]][1].append(df_near.end.iloc[i])\n",
    "            unique_set[df_near.gene_name.iloc[i]][2].append(i)\n",
    "            unique_set[df_near.gene_name.iloc[i]][3].append(df_near.sequence.iloc[i])\n",
    "        except KeyError:\n",
    "            unique_set[df_near.gene_name.iloc[i]] = [[df_near.start.iloc[i]],\n",
    "                                                     [df_near.end.iloc[i]],\n",
    "                                                     [i],\n",
    "                                                     [df_near.sequence.iloc[i]]]\n",
    "    \n",
    "    # Collapse genes with same name, taking:\n",
    "    #     1) First entry for accession etc.\n",
    "    #     2) Total extent min -> max of set of genes\n",
    "    #     3) Longest protein sequence\n",
    "    output = np.zeros(len(df_near.direction),dtype=np.bool)\n",
    "    for g in unique_set.keys():\n",
    "        \n",
    "        index = unique_set[g][2][-1]\n",
    "        direction = df_near.direction.iloc[index]\n",
    "        \n",
    "        seq_lengths = [(len(s),s) for s in unique_set[g][3]]\n",
    "        seq_lengths.sort()\n",
    "        sequence = seq_lengths[-1][1]\n",
    "        \n",
    "        df_near.iloc[index].sequence = sequence\n",
    "        \n",
    "        if direction > 0:\n",
    "            df_near.iloc[index].start = np.min(unique_set[g][0])\n",
    "            df_near.iloc[index].end = np.max(unique_set[g][1])\n",
    "        if direction < 0:\n",
    "            df_near.iloc[index].end = np.min(unique_set[g][1])\n",
    "            df_near.iloc[index].start = np.max(unique_set[g][0])\n",
    "        \n",
    "        output[index] = True\n",
    "    \n",
    "    df_near = df_near[output]\n",
    "\n",
    "    return df_near\n",
    "    \n",
    "\n",
    "def scaffold_blast(region_df,db):\n",
    "    \"\"\"\n",
    "    Blast a region with a collection of genes against a blast db,\n",
    "    only returning hits that match scaffold name.\n",
    "    \"\"\"\n",
    "\n",
    "    all_hits = []\n",
    "    for i in range(len(region_df.sequence)):\n",
    "        print(i,len(region_df.sequence))\n",
    "\n",
    "        r = tools.blast.local_blast(region_df.iloc[i],db=db,blast_program=\"tblastn\")\n",
    "        all_hits.append(r)\n",
    "        \n",
    "    return all_hits\n",
    "\n",
    "\n",
    "def load_block_lengths(species):\n",
    "    \"\"\"\n",
    "    Load the lengths of the blocks (contigs, chromosomes, etc.) that\n",
    "    are being queried.\n",
    "    \"\"\"\n",
    "    \n",
    "    block_lengths = {}\n",
    "    \n",
    "    report_files = glob.glob(os.path.join(species,\"*assembly_report.txt\"))\n",
    "    for report in report_files:\n",
    "        with open(report) as f:\n",
    "            for line in f.readlines():\n",
    "                if not line.startswith(\"#\"):\n",
    "                    fields = line.split(\"\\t\")\n",
    "                    \n",
    "                    genbank_id = fields[4]\n",
    "                    try:\n",
    "                        seq_length = int(fields[8])\n",
    "                    except ValueError:\n",
    "                        seq_length = np.nan\n",
    "                    \n",
    "                    block_lengths[genbank_id] = seq_length\n",
    "                    \n",
    "    return block_lengths\n",
    "\n",
    "# -------------------------------------------------------------------\n",
    "# Functions for processing BLAST hits.\n",
    "# -------------------------------------------------------------------\n",
    "\n",
    "def merge_overlapping_genes(df_to_sort):\n",
    "    \"\"\"\n",
    "    Merge overlapping genes. \n",
    "    \"\"\"\n",
    "    \n",
    "    # Nothing to do\n",
    "    if len(df_to_sort) < 2:\n",
    "        return df_to_sort\n",
    "    \n",
    "    # Copy df and make a column that will encode the original order\n",
    "    # of the rows\n",
    "    some_df = df_to_sort.copy()\n",
    "    some_df[\"original_order\"] = np.arange(len(some_df))\n",
    "    \n",
    "    # Decide whether the intervals go left to right (not inverted)\n",
    "    # or right to left (inverted).  Encode the smaller values as \"a\"\n",
    "    # and the larger values as \"b\"\n",
    "    if some_df.iloc[0].subject_start < some_df.iloc[0].subject_end:\n",
    "        inverted = False\n",
    "        some_df[\"a\"] = some_df.subject_start\n",
    "        some_df[\"b\"] = some_df.subject_end\n",
    "    else:\n",
    "        inverted = True\n",
    "        some_df[\"b\"] = some_df.subject_start\n",
    "        some_df[\"a\"] = some_df.subject_end\n",
    "    \n",
    "    # Sort by \"a\" (the interval starts)\n",
    "    some_df = some_df.sort_values(by=[\"a\"])\n",
    "    \n",
    "    # Mask that encodes whether we will keep row or not\n",
    "    to_keep_mask = np.ones(len(some_df),dtype=np.bool)\n",
    "    \n",
    "    # Go through all rows...\n",
    "    list_of_ends = []\n",
    "    current_end = some_df.iloc[0].b\n",
    "    for i in range(1,len(some_df)):\n",
    "        \n",
    "        # Does new start overlap previous end? If so, merge it.\n",
    "        if some_df.iloc[i].a < some_df.iloc[i-1].b:\n",
    "            to_keep_mask[i] = False\n",
    "            \n",
    "            # Does end of the new overlap stick past current end?\n",
    "            if some_df.iloc[i].b > current_end:\n",
    "                current_end = some_df.iloc[i].b\n",
    "        else:\n",
    "            list_of_ends.append(current_end)\n",
    "            current_end = some_df.iloc[i].b\n",
    "    \n",
    "    # List will have all ends we've seen\n",
    "    list_of_ends.append(current_end)\n",
    "    \n",
    "    # Whack out duplicates\n",
    "    some_df = some_df[to_keep_mask]\n",
    "    \n",
    "    # Assign new ends to intervals\n",
    "    some_df.b = list_of_ends\n",
    "    \n",
    "    # Assign new \"a\" and \"b\" values back to start/end values\n",
    "    if inverted:\n",
    "        some_df.subject_start = some_df.b\n",
    "        some_df.subject_end = some_df.a\n",
    "    else:\n",
    "        some_df.subject_start = some_df.a\n",
    "        some_df.subject_end = some_df.b\n",
    "    \n",
    "    # Sort by original order\n",
    "    some_df = some_df.sort_values([\"original_order\"])\n",
    "    \n",
    "    # Drop temporary columns and return\n",
    "    return some_df.drop([\"a\",\"b\",\"original_order\"],axis=1)\n",
    "    \n",
    "                \n",
    "    \n",
    "def merge_blast_hits(subject_species_list,\n",
    "                     query_hit_list,\n",
    "                     query_region_list,\n",
    "                     query_species_list,\n",
    "                     e_value_cutoff=0.001,\n",
    "                     merge_pattern=\"TLR4\"):\n",
    "\n",
    "    all_hits = []\n",
    "    for i in range(len(query_region_list)):\n",
    "        \n",
    "        query_region = query_region_list[i]\n",
    "        query_species = query_species_list[i]\n",
    "        query_hits = query_hit_list[i]\n",
    "        \n",
    "        for j in range(len(subject_species_list)):\n",
    "            for k in range(len(query_region.gene_name)):\n",
    "\n",
    "                hit_df = query_hits[j][k].copy()\n",
    "\n",
    "                # Filter by e-value\n",
    "                hit_df = hit_df[hit_df.e_value < e_value_cutoff]\n",
    "                if len(hit_df) == 0:\n",
    "                    continue\n",
    "\n",
    "                hit_df[\"query_species\"] = query_species\n",
    "                hit_df[\"query_name\"] = query_region.gene_name.iloc[k]\n",
    "                hit_df[\"subject_species\"] = subject_species_list[j]\n",
    "                hit_df[\"block\"] = [h.split()[0] for h in hit_df.hit_def]\n",
    "\n",
    "                all_hits.append(hit_df)\n",
    "\n",
    "    all_hits = pd.concat(all_hits,ignore_index=True)\n",
    "    all_hits = all_hits.reset_index(drop=True)\n",
    "    \n",
    "    # Look for hits where the query gene name matches the specified pattern\n",
    "    compiled_pattern = re.compile(merge_pattern,flags=re.IGNORECASE)\n",
    "    pattern_mask = np.array([compiled_pattern.search(q)\n",
    "                             for q in all_hits.query_name],dtype=np.bool)\n",
    "    pattern_hits = all_hits[pattern_mask]\n",
    "    \n",
    "    # Go through unique subject-species/dna block pairs in the pattern hits\n",
    "    # and genes that overlap. \n",
    "    final_hits = []\n",
    "    pattern_hits[\"unique\"] = [\"{}-{}\".format(*pair)\n",
    "                              for pair in zip(pattern_hits.subject_species,pattern_hits.block)]\n",
    "    for k in set(pattern_hits.unique):\n",
    "        if len(pattern_hits[pattern_hits.unique == k]) > 1:\n",
    "            tmp = merge_overlapping_genes(pattern_hits[pattern_hits.unique == k])\n",
    "            tmp = tmp.drop([\"unique\"],axis=1)\n",
    "            tmp[\"query_name\"] = merge_pattern\n",
    "            final_hits.append(tmp)\n",
    "    \n",
    "    # Grab hits from the original data frame that do not match the pattern\n",
    "    final_hits.append(all_hits[np.logical_not(pattern_mask)])\n",
    "    \n",
    "    # Reconstruct the total hit df\n",
    "    all_hits = pd.concat(final_hits,ignore_index=True)\n",
    "    all_hits = all_hits.sort_values(by=[\"query_species\",\"subject_species\",\"block\",\"subject_start\"])\n",
    "    all_hits = all_hits.reset_index(drop=True)\n",
    "    \n",
    "    all_hits[\"uid\"] = [\"\".join([random.choice(string.ascii_letters) for j in range(10)])\n",
    "                                for i in range(len(all_hits))]\n",
    "    \n",
    "    return all_hits\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "# Create lists of genes in the region around TLR4 for human and zebrafish\n",
    "# Extend out so we have 22 genes in each dataset. \n",
    "hs = load_genome(\"reference/Homo_sapiens.GRCh38.pep.all.fa\")\n",
    "hs_region = create_region_list(hs,\"TLR4|toll.like.receptor.4|toll-like receptor 4\",chromosome=9,extend_out=3200000,take_only_longest=False)\n",
    "\n",
    "dr = load_genome(\"reference/Danio_rerio.GRCz11.pep.all.fa\")\n",
    "dr_region = create_region_list(dr,\"TLR4|toll.like.receptor.4|toll-like receptor 4\",chromosome=13,extend_out=500000)\n",
    "\n",
    "hs_region.to_csv(\"hs_region.csv\")\n",
    "dr_region.to_csv(\"dr_region.csv\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Blast genes in reference regions against genomes from different\n",
    "# critters. This is very slow. It saves out all_hs_hits.pickle and\n",
    "# all_dr_hits.pickle. These can be loaded in the next cell, meaning\n",
    "# you only need to run this cell once, even if you tweak the \n",
    "# downstream analysis. \n",
    "\n",
    "hs_hits = []\n",
    "dr_hits = []\n",
    "for s in subject_species:\n",
    "    \n",
    "    block_lengths = load_block_lengths(s)\n",
    "    pickle.dump(block_lengths,\n",
    "                open(\"{}_length.pickle\".format(s),\"wb\"))\n",
    "    \n",
    "    db = \"{}/{}\".format(s,s)\n",
    "    \n",
    "    print(\"human against {}\".format(s))\n",
    "    hs_hits.append(scaffold_blast(hs_region,db))\n",
    "    \n",
    "    print(\"zebrafish against {}\".format(s))\n",
    "    dr_hits.append(scaffold_blast(dr_region,db))\n",
    "\n",
    "pickle.dump(hs_hits,open(\"all_hs_hits.pickle\",\"wb\"))\n",
    "pickle.dump(dr_hits,open(\"all_dr_hits.pickle\",\"wb\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge all unique blast hits into a single .csv file\n",
    "hs_hits = pickle.load(open(\"all_hs_hits.pickle\",\"rb\"))\n",
    "dr_hits = pickle.load(open(\"all_dr_hits.pickle\",\"rb\"))\n",
    "\n",
    "all_hits = merge_blast_hits(subject_species,\n",
    "                 [hs_hits,dr_hits],\n",
    "                 [hs_region,dr_region],\n",
    "                 [\"homo\",\"danio\"])\n",
    "\n",
    "all_hits.to_csv(\"all-hits.csv\")"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}