genome2pathway.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "from Bio import SeqIO, Entrez\n",
    "from Bio.Blast import NCBIWWW, NCBIXML\n",
    "import requests\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing CP000975.1 in GCA_000019665.fna\n"
     ]
    }
   ],
   "source": [
    "Entrez.email = \"u7457260@anu.edu.au\"\n",
    "\n",
    "def read_fasta_file(file_path):\n",
    "    sequences = list(SeqIO.parse(file_path, \"fasta\"))\n",
    "    return sequences\n",
    "\n",
    "def blast_sequence(sequence, program=\"blastx\", database=\"nr\", e_value=0.001, max_length=10000):\n",
    "    result_handles = []\n",
    "    for i in range(0, len(sequence), max_length):\n",
    "        subseq = sequence[i:i+max_length]\n",
    "        result_handle = NCBIWWW.qblast(program, database, subseq, hitlist_size=10, expect=e_value)\n",
    "        result_handles.append(result_handle)\n",
    "    return result_handles\n",
    "\n",
    "def parse_blast_result(blast_result, e_value_threshold=1e-3, max_results=10):\n",
    "    blast_records = NCBIXML.read(blast_result)\n",
    "    best_hits = []\n",
    "\n",
    "    for alignment in blast_records.alignments[:max_results]:\n",
    "        for hsp in alignment.hsps:\n",
    "            if hsp.expect <= e_value_threshold:\n",
    "                best_hits.append((alignment.accession, hsp.expect))\n",
    "    return best_hits\n",
    "\n",
    "def get_kegg_pathways(gene_id):\n",
    "    url = f\"https://www.kegg.jp/dbget-bin/www_bget?{gene_id}\"\n",
    "    response = requests.get(url)\n",
    "\n",
    "    pathways = []\n",
    "    if \"Pathway\" in response.text:\n",
    "        for line in response.text.split(\"\\n\"):\n",
    "            if \"<td><a href=\\\"/dbget-bin/www_bget?\" in line:\n",
    "                pathway_id = line.split(\">\")[-1].split(\"<\")[0]\n",
    "                pathway_name = line.split(\"</a>\")[-1].strip()\n",
    "                pathways.append((pathway_id, pathway_name))\n",
    "    return pathways\n",
    "\n",
    "def get_protein_name_by_gene_id(gene_id):\n",
    "    search_result = Entrez.esearch(db=\"protein\", term=f\"{gene_id}[Accession]\")\n",
    "    result_ids = Entrez.read(search_result)[\"IdList\"]\n",
    "\n",
    "    protein_names = []\n",
    "    for result_id in result_ids:\n",
    "        protein_summary = Entrez.esummary(db=\"protein\", id=result_id)\n",
    "        protein_info = Entrez.read(protein_summary)[0]\n",
    "        protein_names.append(protein_info[\"Title\"])\n",
    "\n",
    "    return protein_names\n",
    "\n",
    "def main(input_folder, file_name):\n",
    "    summary_data = []\n",
    "\n",
    "    file_path = os.path.join(input_folder, file_name)\n",
    "    sequences = read_fasta_file(file_path)\n",
    "\n",
    "    for sequence in sequences:\n",
    "        print(f\"Processing {sequence.id} in {file_name}\")\n",
    "        blast_results = blast_sequence(sequence.seq)\n",
    "        for blast_result in blast_results:\n",
    "            best_hits = parse_blast_result(blast_result)\n",
    "\n",
    "            for gene_id, e_value in best_hits:\n",
    "                pathways = get_kegg_pathways(gene_id)\n",
    "                protein_names = get_protein_name_by_gene_id(gene_id)\n",
    "\n",
    "                for protein_name in protein_names:\n",
    "                    for pathway_id, pathway_name in pathways:\n",
    "                        summary_data.append({\n",
    "                            \"Bacteria\": file_name,\n",
    "                            \"GeneID\": gene_id,\n",
    "                            \"Protein\": protein_name,\n",
    "                            \"PathwayID\": pathway_id,\n",
    "                            \"PathwayName\": pathway_name\n",
    "                        })\n",
    "\n",
    "    summary_df = pd.DataFrame(summary_data)\n",
    "    output_file = f\"genome/{file_name}_summary.csv\"\n",
    "    summary_df.to_csv(output_file, index=False)\n",
    "\n",
    "\n",
    "input_folder = \"genome\"\n",
    "filenames_file = \"genome/name_list.TXT\"\n",
    "\n",
    "with open(filenames_file, 'r', encoding='ISO-8859-1') as file:\n",
    "    filenames = [line.strip() for line in file.readlines()]\n",
    "\n",
    "for file_name in filenames:\n",
    "    main(input_folder, file_name)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}