helper_Transcriptome_EXP_JSD_scores.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69364a45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate transcriptome based EXP (normalized expression) and JSD (tissue specificity) scores\n",
    "# start with transcriptome (RNA-seq) TPM data file (generated with salmon package)\n",
    "# output are EXP (normailzed log2 of TPM+1) and JSD score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5599771",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.stats import entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a1244d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Jensen–Shannon divergence Score (tissue specificity)\n",
    "# https://en.wikipedia.org/wiki/Jensen%E2%80%93Shannon_divergence\n",
    "def JSD(p, q):\n",
    "    \"\"\" Given two array (p,q), calculate Jensen-Shannon Divergence \"\"\"\n",
    "    p = np.asarray(p)\n",
    "    q = np.asarray(q)\n",
    "   # normalize\n",
    "    p /= p.sum()\n",
    "    q /= q.sum()\n",
    "    m = (p + q) / 2\n",
    "    return (entropy(p, m) + entropy(q, m)) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c62f5e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input file: transcriptome (RNA-seq) data process with salmon tximport workflow\n",
    "# https://f1000research.com/articles/4-1521\n",
    "query_df=pd.read_csv(\"salmon_abundance_TPM_geneLevel.txt\",sep=\"\\t\",index_col=0)\n",
    "query_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04cb80df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get backgroud reference RNA-seq data\n",
    "# for mouse RNA-seq data\n",
    "# Option 1) Use Encode project and 7 Organs Devevelopment RNA-seq as reference\n",
    "# this is the prefered reference RNA-seq for embryonic developmental samples\n",
    "ref_df=pd.read_csv(\"data/mouse_references_rnaseq/mouse_rnaseq_encode_7Organs-Dev_tpm_delSimi.txt\",sep=\"\\t\",index_col=0)\n",
    "# Option 2) Use FANTOM5 project cells and tissues as reference\n",
    "# ref_df=pd.read_csv(\"data/mouse_references_rnaseq/mouse_rnaseq_f5cells_f5tissues_tpm_delSimi.txt\",sep=\"\\t\",index_col=0)\n",
    "# for human RNA-seq data\n",
    "# ref_df=pd.read_csv(\"data/human_references_rnaseq/human_encode_rnaseq_combined_tpm.tsv\",sep=\"\\t\",index_col=0)\n",
    "ref_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61f1f565",
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine query data with reference data\n",
    "query_ref_combined=pd.concat([query_df,ref_df], join='inner', axis=1)\n",
    "query_ref_combined.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eac4ea8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save to file\n",
    "query_ref_combined.to_csv(\"query_and_encode_7Organs_salmon_abundance_TPM.txt\",sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43b4647b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# quantile normalization in R, in terminal run:\n",
    "# Rscript log2_quantile_normalization.R -i query_and_encode_7Organs_salmon_abundance_TPM.txt -o query_and_encode_7Organs_salmon_abundance_TPM_log2norm.txt\n",
    "# or run directly in this notebook\n",
    "! Rscript log2_quantile_normalization.R -i query_and_encode_7Organs_salmon_abundance_TPM.txt -o query_and_encode_7Organs_salmon_abundance_TPM_log2norm.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e87ad356",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load from file\n",
    "query_ref_combined=pd.read_csv(\"query_and_encode_7Organs_salmon_abundance_TPM.txt\",sep=\"\\t\",index_col=0)\n",
    "query_ref_combined.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9a8ecec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove similar samples based on Pearson correlation coefficient\n",
    "def removeSimilarSamples(bkg_df, query_col, pcc_cutoff=0.85):\n",
    "    pcc_list=list(np.corrcoef(bkg_df.T, query_col.T)[-1])\n",
    "    remove_columns=[]\n",
    "    for columnName,pcc_val in zip(bkg_df.columns.values, pcc_list[:-1]):\n",
    "        if pcc_val>pcc_cutoff:remove_columns.append(columnName)\n",
    "    if len(remove_columns)>0:\n",
    "        return bkg_df.drop(labels=remove_columns,axis=1)\n",
    "    else:\n",
    "        return bkg_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "943923c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate jsd scores for all genes, loop through all genes\n",
    "def jsdLoop(bkg_df, query_col, pcc_cutoff=0.85):\n",
    "    \"\"\" \n",
    "    rows for bkg_df and query_df should be exactly the same (same genes, same order)\n",
    "    bkg_df have no Gene Symbole or name column, all columns are expression values\n",
    "    query_col should be a pandas series (one column)\n",
    "    \"\"\"\n",
    "    # remove samples in bkg_df that have PCC higher than pcc_cutoff with query sample\n",
    "    bkg_df=removeSimilarSamples(bkg_df, query_col, pcc_cutoff=pcc_cutoff)\n",
    "    query_gene_list=[]\n",
    "    query_jsd_list=[]\n",
    "    gene_val_list=query_col.iteritems()\n",
    "    for rowIdx, (gene,val) in enumerate(query_col.iteritems()):\n",
    "        bkg_arr=bkg_df.iloc[rowIdx,]\n",
    "        P_arr=list(bkg_arr)\n",
    "        P_arr.insert(0,val)\n",
    "        Q_arr=[1.0]+[0.0 for i in range(len(P_arr)-1)]\n",
    "        jsd_score=1.0-JSD(P_arr,Q_arr)\n",
    "        query_gene_list.append(gene)\n",
    "        if not math.isnan(jsd_score):\n",
    "            query_jsd_list.append(jsd_score)\n",
    "        else:\n",
    "            query_jsd_list.append(0)\n",
    "    #sorted_query_gene_jsd_list=sorted(query_gene_jsd_list, key = lambda x: x[1], reverse=True)\n",
    "    return query_gene_list, query_jsd_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eeb72159",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract all genes that are transcription factors\n",
    "# annotation of genome-wide transcription factors (TFs) download from AnimalTFDB database\n",
    "# http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/\n",
    "mouse_TF_df=pd.read_csv(\"data/AnimalTFDB_Mus_musculus_TF.txt\",sep=\"\\t\",index_col=0)\n",
    "query_ref_tf_df=query_ref_combined.loc[list(set(query_ref_combined.index).intersection(set(mouse_TF_df.index)))]\n",
    "query_ref_tf_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e59aa98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# seperate query and reference dataframe\n",
    "query_sample_number=len(query_df.columns)\n",
    "query_tf_df=query_ref_tf_df.iloc[:,:query_sample_number]\n",
    "ref_tf_df=query_ref_tf_df.iloc[:,query_sample_number:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed3f4893",
   "metadata": {},
   "outputs": [],
   "source": [
    "query_tf_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaf9d7b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_tf_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "160695e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate TFs' JSD score for each sample in query_df\n",
    "sample_JSDList_dict={}\n",
    "for sampleName in query_tf_df.columns:\n",
    "    #print(sampleName)\n",
    "    query_col=query_tf_df.loc[:,sampleName]\n",
    "    query_gene_list, query_jsd_list = jsdLoop(ref_tf_df, query_col, pcc_cutoff=0.85)\n",
    "    # apply log2 transform\n",
    "    query_jsd_list=np.log2(query_jsd_list + 0.01)\n",
    "    sample_JSDList_dict[sampleName]=query_jsd_list\n",
    "# add gene column\n",
    "sample_JSDList_dict[\"Gene\"]=query_gene_list\n",
    "# create dataframe\n",
    "sample_JSDList_df=pd.DataFrame(sample_JSDList_dict)\n",
    "sample_JSDList_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "258c7ce2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add TFs annotation information\n",
    "sample_JSDList_df=sample_JSDList_df.set_index(\"Gene\")\n",
    "sample_JSDList_wName_df=pd.concat([sample_JSDList_df,mouse_TF_df], join='inner', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "781aa071",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export to file\n",
    "sample_JSDList_wName_df.to_csv(\"salmon_abundance_TPM_geneLevel_log2norm_JSDscore.txt\",sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b71e8fe5",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}