DisGeNET_similarity_network.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "confused-joining",
   "metadata": {},
   "source": [
    "# **USECASE: Building a disease-disease similairty netwrok**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unlimited-bahrain",
   "metadata": {},
   "source": [
    "* DBRetina is an efficent tool for  building a similarity network for a set of items by pairwaise calcuation of their shared features using a linear-time algorithm.\n",
    "* [DisGeNET](https://www.disgenet.org/) has one of the largest collections of genes associated to human diseases.\n",
    "* In this tutorial, we will use DBRetina to build a disease-disease similairty netwrok based on the shared number of genes asscoiated with them in the DisGeNET database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "handmade-latest",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## Download the most recent version of disease gene associations from DisGeNET\n",
    "if [ ! -f all_gene_disease_associations.tsv ];then\n",
    "  wget --quiet http://www.disgenet.org/static/disgenet_ap1/files/downloads/all_gene_disease_associations.tsv.gz\n",
    "  gunzip all_gene_disease_associations.tsv.gz\n",
    "else echo \"all_gene_disease_associations.tsv file exists in the disgenet DB\";fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "scientific-chair",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "head -n3 all_gene_disease_associations.tsv | tr '\\t' ',' | sed 's/ //g'\n",
    "## Here are some quaility matrices\n",
    "# Score: estimate for the curation and size of evidence \n",
    "# DSI (Disease Specificity Index): Some genes are associated wiht multiple diseases (e.g. TNF) while others are more specific\n",
    "# DPI (Disease Pleiotropy Index): similar but consider the disease similarity based on MeSH classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "elect-direction",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "## transform the data table into the DBRetina format\n",
    "## DBRetina expects 2 files. Both are tab-separated files with two columns. Files must have header lines\n",
    "## 1) Associations file: The 1st column for \"items\" and the 2nd for their asscoiated \"features\".\n",
    "## 2) Super-association file: The 1st column for \"items\" and the 2nd for their \"aliases\". You can use this column to update the item name or if you want to pool multiple items together as one super item otherwise the 2nd column should be a copy of the 1st column\n",
    "## In addition, we will filter the input list to keep trusted disease-gene associations only (DisGeNET score > 0.3)\n",
    "cat all_gene_disease_associations.tsv | sed -e 's/^[ \\t]*//' | awk 'BEGIN{FS=OFS=\"\\t\";}{if($10>0.3)print $6,$2}' > disgenet.asc\n",
    "echo \"item alias\" | tr ' ' '\\t' > disgenet.names\n",
    "tail -n+2 disgenet.asc | awk 'BEGIN{FS=OFS=\"\\t\";}{print $1,$1}' | sort | uniq >> disgenet.names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "other-causing",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "## Let us explore the format of the prepared filed\n",
    "echo \"DisGeNET input file\"\n",
    "wc -l all_gene_disease_associations.tsv\n",
    "echo \"===================\"\n",
    "\n",
    "echo \"Associations file\"\n",
    "wc -l disgenet.asc\n",
    "echo \"-------------------\"\n",
    "head -n3 disgenet.asc\n",
    "echo \"===================\"\n",
    "\n",
    "echo \"Super-associations file\"\n",
    "wc -l disgenet.names\n",
    "echo \"-------------------\"\n",
    "head -n3 disgenet.names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "injured-chorus",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "## Now we can run DBRetina\n",
    "kPro_index=\"disgenetDBR\"\n",
    "DBRetina items_indexing -i disgenet.asc -n disgenet.names -p ${kPro_index}\n",
    "DBRetina pairwise -i ${kPro_index}\n",
    "\n",
    "# DBRetina has 3 output files:\n",
    "# 1) The table of pairwise combinations \"${kPro_index}_kSpider_pairwise.tsv\": The items (i.e. diseases) are encoded as numerical IDs\n",
    "# 2) The names map \"${kPro_index}.namesMap\": It maps each item to its ID\n",
    "# 3) The nodes size \"${kPro_index}_nodes_size.tsv\": It reports the number of features associated with each item\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "single-rotation",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Have a look on the output files\n",
    "kPro_index=\"disgenetDBR\"\n",
    "echo \"The table of pairwise combinations\"\n",
    "echo \"----------------------------------\"\n",
    "head -n3 ${kPro_index}_kSpider_pairwise.tsv\n",
    "echo \"----------------------------------\"\n",
    "echo \"The table of names map\"\n",
    "echo \"----------------------\"\n",
    "head -n5 ${kPro_index}.namesMap\n",
    "echo \"----------------------\"\n",
    "echo \"The table of nodes sizes\"\n",
    "echo \"----------------------\"\n",
    "head -n3 ${kPro_index}_kSpider_seqToKmersNo.tsv\n",
    "echo \"----------------------\"\n",
    "# How many pairwise combinations did we do?\n",
    "echo \"The number of pairwise combinations =\" $(wc -l ${kPro_index}_kSpider_pairwise.tsv | cut -d\" \" -f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accurate-success",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Let us pool the output files into one final output \n",
    "# 1. merge the names map and nodes sizes in one output file\n",
    "kPro_index=\"disgenetDBR\"\n",
    "paste <(tail -n+2 ${kPro_index}.namesMap |cut -d\" \" -f1)  <(tail -n+2 ${kPro_index}.namesMap |cut -d\" \" -f2-) > ${kPro_index}.namesMap.tmp\n",
    "echo \"node_id node_name size\" | tr ' ' '\\t' > ${kPro_index}_nodes_size.tsv\n",
    "awk 'BEGIN{FS=OFS=\"\\t\";}FNR==NR{a[$2]=$3;next;}{if(a[$1]!=\"\")print $0,a[$1]}' ${kPro_index}_kSpider_seqToKmersNo.tsv ${kPro_index}.namesMap.tmp >> ${kPro_index}_nodes_size.tsv\n",
    "rm ${kPro_index}.namesMap.tmp*\n",
    "\n",
    "# Let us have a look\n",
    "head -n5 ${kPro_index}_nodes_size.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rapid-encoding",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# 2. Annotation of the numerically coded association file\n",
    "# a) Add items names and no of associated items \n",
    "# b) calc jaccard distance and containment ratio for each pair\n",
    "# c) Additionally, we can filter out those with minimal similarities \n",
    "kPro_index=\"disgenetDBR\"\n",
    "md=1   ## minimum jaccard distance (as a percentage) to keep the record \n",
    "mc=10  ## minimum containment ratio (as a percentage) to keep the record \n",
    "echo \":START_ID|START_name|START_size|shared_count:int|jDist:float|smPerc:float|END_name|END_size|:END_ID\" > ${kPro_index}_relations.csv\n",
    "awk -v md=$md -v mc=$mc 'BEGIN{FS=\"\\t\";S=\"|\";}FNR==NR{a[$1]=$3;b[$1]=$2S$3;next;}{\n",
    "   g1=a[$2]; g2=a[$3]; min=g1;min=(min < g2 ? min : g2); \n",
    "   jDist=$4*100/(g1+g2-$4); smPerc=$4*100/min; \n",
    "   if(jDist>md || smPerc>mc)\n",
    "     printf(\"%s%s%s%s%s%s%.1f%s%.1f%s%s%s%s\\n\", $2,S,b[$2],S,$4,S,jDist,S,smPerc,S,b[$3],S,$3)}' \\\n",
    "   ${kPro_index}_nodes_size.tsv <(tail -n+2 ${kPro_index}_kSpider_pairwise.tsv) >> ${kPro_index}_relations.csv\n",
    "\n",
    "# Let us have a look\n",
    "wc -l ${kPro_index}_relations.csv \n",
    "echo \"----------------------\"\n",
    "head -n5 ${kPro_index}_relations.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "stock-blackjack",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Query for diseases with high similarity to any Alzheimer disease\n",
    "kPro_index=\"disgenetDBR\"\n",
    "cat ${kPro_index}_relations.csv | grep -i \"Alzheimer\" | awk -F\"|\" '{if($5>10 && $6>25)print}' > ${kPro_index}_relations_Alzheimer.csv\n",
    "# Check the impact of filtration\n",
    "wc -l ${kPro_index}_relations.csv\n",
    "echo \"----------------------\"\n",
    "wc -l ${kPro_index}_relations_Alzheimer.csv\n",
    "echo \"----------------------\"\n",
    "# update the IDs to be from 1 to n where n is the number of diseases in this small file\n",
    "tail -n+2 ${kPro_index}_relations_Alzheimer.csv | awk 'BEGIN{FS=\"|\";OFS=\"\\n\";}{print $1\"|\"$2,$9\"|\"$7}' | sort -n | uniq |\\\n",
    "awk 'BEGIN{FS=OFS=\"|\";}{print $1,NR,$2}' > id.map\n",
    "wc -l id.map\n",
    "echo \"----------------------\"\n",
    "head -n5 id.map\n",
    "echo \"----------------------\"\n",
    "# Collect all pairwise distances of the diseases with high similarity to Alzheimer\n",
    "head -n1 ${kPro_index}_relations.csv > ${kPro_index}_relations_newIds.csv\n",
    "awk 'BEGIN{FS=\"|\";OFS=\"|\";}FNR==NR{a[$1]=$2;next;}{if(a[$1]&&a[$9]){$1=a[$1];$9=a[$9];print $0;}}' id.map ${kPro_index}_relations.csv >> ${kPro_index}_relations_newIds.csv\n",
    "wc -l ${kPro_index}_relations_newIds.csv\n",
    "echo \"----------------------\"\n",
    "head -n5 ${kPro_index}_relations_newIds.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alpine-copying",
   "metadata": {},
   "outputs": [],
   "source": [
    "# enables the %%R magic\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dramatic-channel",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R \n",
    "if (!requireNamespace(\"RFLPtools\", quietly = TRUE))\n",
    "    install.packages(\"RFLPtools\")\n",
    "if (!requireNamespace(\"pheatmap\", quietly = TRUE))\n",
    "    install.packages(\"pheatmap\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adequate-concentrate",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R \n",
    "kPro_index=\"disgenetDBR\";\n",
    "target=\"smPerc\"\n",
    "\n",
    "rels=paste(kPro_index,\"relations_newIds.csv\",sep=\"_\");\n",
    "dbretina_rel=read.table(file=rels ,skip=1, header=FALSE, sep = \"|\", quote = \"\", dec = \".\"); \n",
    "colnames(dbretina_rel)=c(\"START_ID\",\"START_name\",\"START_size\",\"count\",\"jDist\",\"smPerc\",\"END_name\",\"END_size\",\"END_ID\")\n",
    "\n",
    "dbretina_nodes=read.table(file=\"id.map\", header=FALSE, sep = \"|\", quote = \"\");\n",
    "size=dim(dbretina_nodes)[1]\n",
    "dbret_count_matrix = matrix(0, nrow = size, ncol = size)\n",
    "for (i in 1:dim(dbretina_rel)[1]){\n",
    "  dbret_count_matrix[dbretina_rel$START_ID[i],dbretina_rel$END_ID[i]]=log2(dbretina_rel[[target]][i])\n",
    "  dbret_count_matrix[dbretina_rel$END_ID[i],dbretina_rel$START_ID[i]]=log2(dbretina_rel[[target]][i])\n",
    "}\n",
    "\n",
    "## transform the similarity matrix into distance object (measure the dissimilarity)\n",
    "\n",
    "library(RFLPtools)\n",
    "dbret_count_dist = sim2dist(dbret_count_matrix, log2(100))\n",
    "## Transform the dist object back into into matrix\n",
    "dbret_count_matrix = as.matrix(dbret_count_dist)\n",
    "\n",
    "rownames(dbret_count_matrix) <- dbretina_nodes$V3\n",
    "colnames(dbret_count_matrix) <- NULL\n",
    "\n",
    "library(pheatmap)\n",
    "#library(repr)\n",
    "\n",
    "colors <- colorRampPalette( rev(brewer.pal(9, \"Blues\")) )(255) ## specify a blue color palette for painting\n",
    "#options(repr.plot.width = 4, repr.plot.height = 2)\n",
    "#jpeg(\"outfig\")\n",
    "pheatmap(dbret_count_matrix,\n",
    "         clustering_distance_rows = dbret_count_dist,\n",
    "         clustering_distance_cols = dbret_count_dist,\n",
    "         col = colors, cellwidth=5,cellheight=8, fontsize_row=7.5,fontsize_col=0)\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "traditional-massage",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "kPro_index=\"disgenetDBR\"\n",
    "#grep \"ACNE\" ${kPro_index}_nodes_size.tsv\n",
    "cat ${kPro_index}_relations_newIds.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hawaiian-corps",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}