genelist2sequence.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# genelist to gene sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:84: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:85: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:92: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:93: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:96: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:97: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "/opt/anaconda/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:99: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "blast_column_names = [\"qseqid\", \"sseqid\", \"pident\", \"length\", \"mismatch\", \"gapopen\", \"qstart\", \"qend\", \"sstart\", \"send\", \"evalue\", \"bitscore\"]                         \n",
    "gff_file = '/bioinfo/palm/ref/dura/dura_ref.gff'\n",
    "gff_column_names = ['Seqid', 'Source', 'Type', 'Start', 'End', 'Score', 'Strand', 'Phase', 'Attributes']\n",
    "bs = 0 #bitscore filtering\n",
    "n = 3 #hits per gene\n",
    "flanking = 2000 #flanking regions upstream and downstream\n",
    "fileA = 'ABA_dura_blastp.tbl'\n",
    "fileB = 'dura_ABA_blastp.tbl'\n",
    "\n",
    "# Read data from file A\n",
    "fileA_data = pd.read_csv(fileA, sep='\\t', names=blast_column_names)\n",
    "\n",
    "# Read data from file B\n",
    "fileB_data = pd.read_csv(fileB, sep='\\t', names=blast_column_names)\n",
    "\n",
    "# Initialize an empty list to store the results\n",
    "result_pairs = {}\n",
    "\n",
    "# Loop through unique values in 'sseqid' column in file A\n",
    "\n",
    "for sseqid_a in fileA_data['sseqid'].unique():\n",
    "    count = 0\n",
    "    # Filter rows in file A where 'sseqid' is equal to the current value in the loop\n",
    "    filtered_rows = fileA_data[fileA_data['sseqid'] == sseqid_a]\n",
    "    # Sort the filtered rows by 'bitscore' in descending order\n",
    "    sorted_rows = filtered_rows.sort_values(by='bitscore', ascending=False)\n",
    "    # Iterate through the sorted rows\n",
    "    count = 0\n",
    "    for _, max_bitscore_row in sorted_rows.iterrows():\n",
    "        if count < n:\n",
    "            # Get the value for 'qseqid' in the current row\n",
    "            qseqid_b = max_bitscore_row['qseqid']\n",
    "\n",
    "            # Check if the row exists in file B with 'sseqid' == 'b' and 'qseqid' == 'a'\n",
    "            matching_row = fileB_data[(fileB_data['sseqid'] == qseqid_b) & (fileB_data['qseqid'] == sseqid_a)]\n",
    "\n",
    "            # If the matching row exists and 'bitscore' is greater than the threshold, add 'a' and 'b' to the result_pairs dictionary\n",
    "            if not matching_row.empty and max_bitscore_row['bitscore'] > bs:\n",
    "                sseqid_a = sseqid_a.split('-mRNA')[0]\n",
    "                qseqid_b = qseqid_b.split('-mRNA')[0]\n",
    "                if sseqid_a not in result_pairs:\n",
    "                    result_pairs[sseqid_a] = []\n",
    "                elif qseqid_b not in result_pairs[sseqid_a]:\n",
    "                    result_pairs[sseqid_a].append(qseqid_b)                    \n",
    "                    count += 1\n",
    "    \n",
    "# get the local gene list\n",
    "def invert_dict_with_lists(original_dict):\n",
    "    inverted_dict = {}\n",
    "\n",
    "    for key, values in original_dict.items():\n",
    "        for value in values:\n",
    "            if value not in inverted_dict:\n",
    "                inverted_dict[value] = []\n",
    "            inverted_dict[value].append(key)\n",
    "\n",
    "    return inverted_dict\n",
    "\n",
    "\n",
    "inverted_dict = invert_dict_with_lists(result_pairs)          \n",
    "string_dict = {key: '/'.join(values) for key, values in inverted_dict.items()}\n",
    "# Now we need to get the position of genes from the gff file.      \n",
    "\n",
    "# Read data from the file\n",
    "df = pd.read_csv(gff_file, sep='\\t',comment='#', names=gff_column_names)\n",
    "# Function to parse the 'Attributes' column and convert it into a dictionary of key-value pairs\n",
    "def parse_attributes(attribute_str):\n",
    "    attribute_pairs = attribute_str.split(';')\n",
    "    attribute_dict = {}\n",
    "    for pair in attribute_pairs:\n",
    "        if len(pair.split('=')) == 2:\n",
    "            key, value = pair.split('=')\n",
    "        else:\n",
    "            key = pair.split('=')[0]\n",
    "            value = ''\n",
    "        attribute_dict[key] = value\n",
    "    return attribute_dict\n",
    "df1 = df['Attributes'].apply(parse_attributes).apply(pd.Series)\n",
    "file_data = pd.concat([df.drop('Attributes', axis=1), df1], axis=1)\n",
    "\n",
    "# Filter rows where value from column A is in the list\n",
    "filtered_rows = file_data[file_data['ID'].isin(list(string_dict))]\n",
    "filtered_rows['genes'] = filtered_rows['ID'].map(string_dict)\n",
    "filtered_rows['ID_genes'] = filtered_rows['ID'] + '_' + filtered_rows['genes']\n",
    "\n",
    "# Extract values from column B and C from the filtered rows\n",
    "result_values = filtered_rows[['Seqid', 'Start','End','ID_genes']]\n",
    "\n",
    "\n",
    "# Convert columns A and B to numeric\n",
    "result_values['Start'] = pd.to_numeric(result_values['Start'])\n",
    "result_values['End'] = pd.to_numeric(result_values['End'])\n",
    "\n",
    "# Add flanking regions\n",
    "result_values['Start'] = result_values['Start'] - flanking\n",
    "result_values['Start'] = result_values['Start'].apply(lambda x: max(0, x))  # Set negative values to 0\n",
    "\n",
    "result_values['End'] = result_values['End'] + flanking\n",
    "bedfile = fileA.split('_blastp.tbl')[0] + '_' + str(flanking) + '.bed'\n",
    "result_values.to_csv(bedfile, sep='\\t', index=False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genelist to cds (dura_transcripts.fasta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "blast_column_names = [\"qseqid\", \"sseqid\", \"pident\", \"length\", \"mismatch\", \"gapopen\", \"qstart\", \"qend\", \"sstart\", \"send\", \"evalue\", \"bitscore\"]                         \n",
    "gff_file = '/bioinfo/palm/ref/dura/dura_ref.gff'\n",
    "transcript_file = '/bioinfo/palm/ref/dura/dura_transcripts.fasta'\n",
    "gff_column_names = ['Seqid', 'Source', 'Type', 'Start', 'End', 'Score', 'Strand', 'Phase', 'Attributes']\n",
    "bs = 0 #bitscore filtering\n",
    "n = 3 #hits per gene\n",
    "flanking = 2000 #flanking regions upstream and downstream\n",
    "fileA = 'ABA_dura_blastp.tbl'\n",
    "fileB = 'dura_ABA_blastp.tbl'\n",
    "# Read data from file A\n",
    "fileA_data = pd.read_csv(fileA, sep='\\t', names=blast_column_names)\n",
    "\n",
    "# Read data from file B\n",
    "fileB_data = pd.read_csv(fileB, sep='\\t', names=blast_column_names)\n",
    "\n",
    "# Initialize an empty list to store the results\n",
    "result_pairs = {}\n",
    "\n",
    "# Loop through unique values in 'sseqid' column in file A\n",
    "\n",
    "for sseqid_a in fileA_data['sseqid'].unique():\n",
    "    count = 0\n",
    "    # Filter rows in file A where 'sseqid' is equal to the current value in the loop\n",
    "    filtered_rows = fileA_data[fileA_data['sseqid'] == sseqid_a]\n",
    "    # Sort the filtered rows by 'bitscore' in descending order\n",
    "    sorted_rows = filtered_rows.sort_values(by='bitscore', ascending=False)\n",
    "    # Iterate through the sorted rows\n",
    "    count = 0\n",
    "    for _, max_bitscore_row in sorted_rows.iterrows():\n",
    "        if count < n:\n",
    "            # Get the value for 'qseqid' in the current row\n",
    "            qseqid_b = max_bitscore_row['qseqid']\n",
    "\n",
    "            # Check if the row exists in file B with 'sseqid' == 'b' and 'qseqid' == 'a'\n",
    "            matching_row = fileB_data[(fileB_data['sseqid'] == qseqid_b) & (fileB_data['qseqid'] == sseqid_a)]\n",
    "\n",
    "            # If the matching row exists and 'bitscore' is greater than the threshold, add 'a' and 'b' to the result_pairs dictionary\n",
    "            if not matching_row.empty and max_bitscore_row['bitscore'] > bs:\n",
    "                if sseqid_a not in result_pairs:\n",
    "                    result_pairs[sseqid_a] = []\n",
    "                elif qseqid_b not in result_pairs[sseqid_a]:\n",
    "                    result_pairs[sseqid_a].append(qseqid_b)                    \n",
    "                    count += 1\n",
    "    \n",
    "# get the local gene list\n",
    "def invert_dict_with_lists(original_dict):\n",
    "    inverted_dict = {}\n",
    "\n",
    "    for key, values in original_dict.items():\n",
    "        for value in values:\n",
    "            if value not in inverted_dict:\n",
    "                inverted_dict[value] = []\n",
    "            inverted_dict[value].append(key)\n",
    "\n",
    "    return inverted_dict\n",
    "\n",
    "\n",
    "inverted_dict = invert_dict_with_lists(result_pairs)          \n",
    "string_dict = {key: '/'.join(values) for key, values in inverted_dict.items()}\n",
    "# Now we get the transcript sequences for those genes (key in string_dict) from dura_transcripts.fasta\n",
    "from Bio import SeqIO\n",
    "out_file = '/bioinfo/palm/ref/ABA/ABA_dura_transcripts.fasta'\n",
    "with open(out_file, 'w') as out_fh:\n",
    "    for record in SeqIO.parse(transcript_file, \"fasta\"):\n",
    "        if record.id in string_dict:\n",
    "            record.description = string_dict[record.id]\n",
    "            SeqIO.write(record, out_fh, \"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Egu017379-mRNA-1': 'LEC1',\n",
       " 'Egu002768-mRNA-1': 'LEC1',\n",
       " 'Egu002765-mRNA-1': 'LEC1',\n",
       " 'Egu023380-mRNA-1': 'BBM',\n",
       " 'Egu002214-mRNA-2': 'BBM',\n",
       " 'Egu000576-mRNA-1': 'BBM',\n",
       " 'Egu019891-mRNA-5': 'ABI5',\n",
       " 'Egu023685-mRNA-1': 'ABI5',\n",
       " 'Egu015962-mRNA-1': 'ABI5',\n",
       " 'Egu019493-mRNA-2': 'ABI3/LEC2/FUS3',\n",
       " 'Egu006921-mRNA-1': 'ABI3/LEC2',\n",
       " 'Egu020554-mRNA-1': 'ABI3/LEC2/FUS3',\n",
       " 'Egu019493-mRNA-1': 'FUS3',\n",
       " 'Egu000535-mRNA-1': 'AGL15',\n",
       " 'Egu003576-mRNA-1': 'AGL15',\n",
       " 'Egu017102-mRNA-2': 'AGL15',\n",
       " 'Egu009827-mRNA-1': 'WUS',\n",
       " 'Egu012330-mRNA-1': 'WUS',\n",
       " 'Egu028948-mRNA-1': 'WUS',\n",
       " 'Egu024837-mRNA-1': 'LEA'}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "string_dict\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "string_dict = {key: '/'.join(values) for key, values in inverted_dict.items()}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Egu017379',\n",
       " 'Egu002768',\n",
       " 'Egu002765',\n",
       " 'Egu023380',\n",
       " 'Egu002214',\n",
       " 'Egu000576',\n",
       " 'Egu019891',\n",
       " 'Egu023685',\n",
       " 'Egu015962',\n",
       " 'Egu019493',\n",
       " 'Egu006921',\n",
       " 'Egu020554',\n",
       " 'Egu005708',\n",
       " 'Egu000535',\n",
       " 'Egu003576',\n",
       " 'Egu017102',\n",
       " 'Egu009827',\n",
       " 'Egu012330',\n",
       " 'Egu028948',\n",
       " 'Egu024837']"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(string_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'LEC1': ['Egu015173', 'Egu017379', 'Egu002768'],\n",
       " 'BBM': ['Egu002214', 'Egu023380', 'Egu002214'],\n",
       " 'ABI5': ['Egu019891', 'Egu019891', 'Egu023685'],\n",
       " 'ABI3': ['Egu019493', 'Egu019493', 'Egu006921'],\n",
       " 'LEC2': ['Egu019493', 'Egu006921', 'Egu020554'],\n",
       " 'FUS3': ['Egu006921', 'Egu020554', 'Egu019493'],\n",
       " 'AGL15': ['Egu000535', 'Egu000535', 'Egu003576'],\n",
       " 'WUS': ['Egu026907', 'Egu009827', 'Egu012330'],\n",
       " 'LEA': ['Egu005882', 'Egu024837']}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_pairs"
   ]
  },
  {
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}