jupyter/20190130_Cg_Giga_cont_AA.fa_BLASTP_uniprot_swprot2019.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2019-01-30 13:03:30--  ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz\n",
      "           => ‘/Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta.gz’\n",
      "Resolving ftp.uniprot.org... 141.161.180.197\n",
      "Connecting to ftp.uniprot.org|141.161.180.197|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /pub/databases/uniprot/current_release/knowledgebase/complete ... done.\n",
      "==> SIZE uniprot_sprot.fasta.gz ... 88385361\n",
      "==> PASV ... done.    ==> RETR uniprot_sprot.fasta.gz ... done.\n",
      "Length: 88385361 (84M) (unauthoritative)\n",
      "\n",
      "uniprot_sprot.fasta 100%[===================>]  84.29M  6.85MB/s    in 11s     \n",
      "\n",
      "2019-01-30 13:03:43 (7.44 MB/s) - ‘/Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta.gz’ saved [88385361]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz -O /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gunzip -k /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Building a new DB, current time: 01/30/2019 13:11:04\n",
      "New DB name:   /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_db\n",
      "New DB title:  /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta\n",
      "Sequence type: Protein\n",
      "Keep MBits: T\n",
      "Maximum file size: 1000000000B\n",
      "Adding sequences from FASTA; added 559077 sequences in 55.0913 seconds.\n"
     ]
    }
   ],
   "source": [
    "!/Users/Shared/bioinformatics/ncbi-blast-2.5.0+/bin/makeblastdb \\\n",
    "-in /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta \\\n",
    "-out /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_db \\\n",
    "-dbtype prot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Confirm Steven's .fa file is the same as Sean's\n",
    "!diff /Volumes/web/halfshell/bu-git-repos/nb-2017/C_gigas/data/Cg_Giga_cont_AA.fa \\\n",
    "/Volumes/web/scaphapoda/Sean/Rhonda-2016-Oyster-Intermediates/Cg-Giga_cont_AA.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#can sort database before aligning to change the order of entries. We don't care that much about taxanomic group, we car about protein functionality. So it doesn't matter for us if human gets called over mouse or vice-versa. So I am not worrying about changing the max_target_seqs parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "!/Users/Shared/bioinformatics/ncbi-blast-2.5.0+/bin/blastp \\\n",
    "-query /Volumes/web/halfshell/bu-git-repos/nb-2017/C_gigas/data/Cg_Giga_cont_AA.fa \\\n",
    "-db /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_db \\\n",
    "-out /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab \\\n",
    "-num_threads 8 \\\n",
    "-max_hsps 1 \\\n",
    "-max_target_seqs 1 \\\n",
    "-outfmt 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   38294 /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab\r\n"
     ]
    }
   ],
   "source": [
    "!wc -l /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   40751\r\n"
     ]
    }
   ],
   "source": [
    "!awk '{if($1~/>/) print $1}' /Volumes/web/halfshell/bu-git-repos/nb-2017/C_gigas/data/Cg_Giga_cont_AA.fa | \\\n",
    "wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#check how many query sequences are not in the Blast output\n",
    "#first extract all sequence names from query file with first awk command\n",
    "#sed command removes the carrot \n",
    "#compare query file sequence names to names blast output\n",
    "#this awk command creates an array with the first file (blast output)\n",
    "#the array creates a key out of the first field 'a[$1]' and stores a value \n",
    "#which in this case is the same as the key 'a[$1]=$1'\n",
    "#;next tells it do this for each line in the first file (blast output)\n",
    "#!($1 in a) means if the first field in the second file is not in the array we made\n",
    "#print the first field of the second file\n",
    "#The second file is denoted as '-' because we are passing it the file we made in the previous sed and awk commands\n",
    "#add a column to the unmapped sequence names that specifies the sequences were unmapped\n",
    "\n",
    "!awk '{if($1~/>/) print $1}' /Volumes/web/halfshell/bu-git-repos/nb-2017/C_gigas/data/Cg_Giga_cont_AA.fa | \\\n",
    "sed 's/>//g'|\\\n",
    "awk 'NR==FNR{a[$1]=$1;next}!($1 in a){print $1\"\\t\"\"unmapped\"}' \\\n",
    "/Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab - \\\n",
    "| tr '|' '.' > /Users/strigg/Desktop/gigas_blastP_uniprot/query_seq_names_NOTin_Blast_out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    2457 /Users/strigg/Desktop/gigas_blastP_uniprot/query_seq_names_NOTin_Blast_out\r\n"
     ]
    }
   ],
   "source": [
    "#count the number of sequence names in the query file that are NOT in the blast output\n",
    "!wc -l /Users/strigg/Desktop/gigas_blastP_uniprot/query_seq_names_NOTin_Blast_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "559077 sp|\r\n"
     ]
    }
   ],
   "source": [
    "#formatting output so that uniprot accessions are in one column\n",
    "#first test to see if all uniprot db entries start with the same pattern \"sp|\"\n",
    "#This awk command finds every line containing a '>' and prints the first three characters after the carrot\n",
    "#Awk command literally: if first field contains a carrot, print a substring of the first field string, starting at the second character in the string, that is 3 characters long\n",
    "# sort the list of substrings and print only unique substrings and their frequencies\n",
    "#if there is only one substring, we can conclude all uniprot DB entries start with the same pattern \"sp|\"\n",
    "!awk '{if($1~/>/)print substr($1,2,3)}' /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta | sort | uniq -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "559077 1\r\n"
     ]
    }
   ],
   "source": [
    "#Next check that there is only 1 pipe in each uniprot DB name after removing the \"sp|\" pattern\n",
    "#we tell awk that our file is tab delimited, then print only the sequence name\n",
    "#then sed to remove the 'sp|' pattern\n",
    "#then sed to replace all non '|' characters with nothing and awk to count the number of characters remaining in each line\n",
    "#https://unix.stackexchange.com/questions/18736/how-to-count-the-number-of-a-specific-character-in-each-line\n",
    "#then sort and unique to see if there are different number of pipes in each line\n",
    "!awk -F\"\\t\" '{if($1~/>/)print $1}' \\\n",
    "/Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta | \\\n",
    "sed 's/sp\\|//g' | sed 's/[^|]//g' | awk '{print length}' | \\\n",
    "sort | uniq -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Next check that the uniprot DB does not contain the pattern \"|m\"\n",
    "!grep '\\|m' /Users/strigg/Desktop/gigas_blastP_uniprot/uniprot_swprot2019_01.fasta |\\\n",
    "head\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">CHOYP_043R.1.5|m.16874\r\n",
      ">CHOYP_043R.5.5|m.64252\r\n",
      ">CHOYP_14332.1.2|m.5643\r\n",
      ">CHOYP_14332.1.2|m.5644\r\n",
      ">CHOYP_14332.2.2|m.61737\r\n",
      ">CHOYP_14332.2.2|m.61738\r\n",
      ">CHOYP_1433E.1.2|m.3639\r\n",
      ">CHOYP_1433E.1.2|m.3638\r\n",
      ">CHOYP_1433E.2.2|m.63376\r\n",
      ">CHOYP_1433E.2.2|m.63378\r\n"
     ]
    }
   ],
   "source": [
    "#And that only the C_gigas sequence names contain \"|m\"\n",
    "#We want to check this so that we can replace \"\\m\" with \".m\"\n",
    "!grep '\\|m' /Volumes/web/halfshell/bu-git-repos/nb-2017/C_gigas/data/Cg_Giga_cont_AA.fa | \\\n",
    "head\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now in our BLAST output file, we can replace 'sp|' with nothing, \"|m\" with \".m\", and then finally we can replace the only remaining pipe with tabs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CHOYP_043R.1.5.m.16874\tQ06852\tSLAP1_CLOTH\t56.944\t216\t61\t18\t10\t197\t1388\t1599\t5.20e-08\t55.8\r\n",
      "CHOYP_043R.5.5.m.64252\tQ06852\tSLAP1_CLOTH\t52.381\t294\t80\t24\t575\t816\t1351\t1636\t2.02e-16\t88.2\r\n",
      "CHOYP_14332.1.2.m.5643\tQ2F637\t1433Z_BOMMO\t66.031\t262\t74\t2\t19\t280\t1\t247\t2.78e-119\t344\r\n",
      "CHOYP_14332.1.2.m.5644\tP62325\tBTG1_MOUSE\t47.205\t161\t80\t2\t1\t156\t11\t171\t2.22e-47\t155\r\n",
      "CHOYP_14332.2.2.m.61737\tQ2F637\t1433Z_BOMMO\t67.331\t251\t78\t1\t1\t251\t1\t247\t4.22e-119\t342\r\n",
      "CHOYP_14332.2.2.m.61738\tQ8IDX6\tRBP2A_PLAF7\t27.778\t108\t76\t1\t31\t136\t2766\t2873\t5.98e-04\t43.9\r\n",
      "CHOYP_1433E.1.2.m.3639\tQ9CWP8\tDPOD4_MOUSE\t38.235\t102\t61\t2\t30\t130\t7\t107\t9.73e-19\t78.6\r\n",
      "CHOYP_1433E.1.2.m.3638\tP92177\t1433E_DROME\t77.692\t260\t53\t1\t1\t255\t1\t260\t1.95e-149\t420\r\n",
      "CHOYP_1433E.2.2.m.63376\tP92177\t1433E_DROME\t77.692\t260\t53\t1\t1\t255\t1\t260\t1.95e-149\t420\r\n",
      "CHOYP_1433E.2.2.m.63378\tQ6MG82\tPRRT1_RAT\t44.444\t99\t53\t1\t41\t137\t196\t294\t1.84e-18\t82.8\r\n"
     ]
    }
   ],
   "source": [
    "#Because all uniprot DB entries start with the same pattern, we can replace it with nothing so that each entry starts with the uniprot accession\n",
    "!sed 's/sp\\|//g' /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab | \\\n",
    "awk -F\"\\t\" '{sub(/\\|/,\".\",$1);print $0}'  \\\n",
    "| tr ' ' '\\t' | tr '|' '\\t' | \\\n",
    "head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CYC_HORSE.\tP00004\tCYC_HORSE\t100.000\t105\t0\t0\t1\t105\t1\t105\t3.35e-73\t215\r\n",
      "ALDOA_RABIT.\tP00883\tALDOA_RABIT\t100.000\t360\t0\t0\t1\t360\t1\t360\t0.0\t745\r\n",
      "LYSC_CHICK.\tP00698\tLYSC_CHICK\t100.000\t147\t0\t0\t1\t147\t1\t147\t1.30e-105\t300\r\n",
      "MYG_HORSE.\tP68082\tMYG_HORSE\t100.000\t154\t0\t0\t1\t154\t1\t154\t6.48e-109\t309\r\n",
      "OVAL_CHICK.\tP01012\tOVAL_CHICK\t100.000\t386\t0\t0\t1\t386\t1\t386\t0.0\t802\r\n",
      "GFP_AEQVI.\tP42212\tGFP_AEQVI\t100.000\t238\t0\t0\t1\t238\t1\t238\t1.26e-179\t494\r\n",
      "SRPP_HEVBR.\tO82803\tSRPP_HEVBR\t100.000\t204\t0\t0\t1\t204\t1\t204\t2.12e-144\t403\r\n",
      "REF_HEVBR.\tP15252\tREF_HEVBR\t100.000\t138\t0\t0\t1\t138\t1\t138\t1.14e-96\t277\r\n",
      "PLMP_GRIFR.\tP81054\tPLMP_GRIFR\t100.000\t348\t0\t0\t1\t348\t1\t348\t0.0\t715\r\n",
      "KKA1_ECOLX\tP00551\tKKA1_ECOLX\t100.000\t271\t0\t0\t1\t271\t1\t271\t0.0\t565\r\n"
     ]
    }
   ],
   "source": [
    "#check end of file too\n",
    "!sed 's/sp\\|//g' /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab | \\\n",
    "awk -F\"\\t\" '{sub(/\\|/,\".\",$1);print $0}'  \\\n",
    "| tr ' ' '\\t'| tr '|' '\\t' | tail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed 's/sp\\|//g' /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.tab | \\\n",
    "awk -F\"\\t\" '{sub(/\\|/,\".\",$1);print $0}'  \\\n",
    "| tr ' ' '\\t' | tr '|' '\\t' \\\n",
    "> /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.nopipe.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat /Users/strigg/Desktop/gigas_blastP_uniprot/giga-uniprot-blastP-out.nopipe.tab \\\n",
    "/Users/strigg/Desktop/gigas_blastP_uniprot/query_seq_names_NOTin_Blast_out \\\n",
    "> /Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   40751 /Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.tab\r\n"
     ]
    }
   ],
   "source": [
    "!wc -l /Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.tab\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#put place holders in fields for unmapped IDs\n",
    "!awk -F\"\\t\" '{if($2==\"unmapped\")print $1FS$2FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\"FS\"@\";else print $0}' \\\n",
    "/Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.tab \\\n",
    "> /Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.placeholder.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#combine Uniprot annotations with BLAST'd proteins\n",
    "!awk -F\"\\t\" 'NR==FNR{a[$1]=$0;next}{print $0\"\\t\"a[$2]}' \\\n",
    "/Users/strigg/Desktop/gigas_blastP_uniprot/uniprot-reviewed%3Ayes.tab \\\n",
    "/Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.placeholder.tab \\\n",
    "> /Users/strigg/Desktop/gigas_blastP_uniprot/all_giga-uniprot-blastP-out.nopipe.annotations.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mv directory onto Metacarcinus on Gannet\n",
    "!scp -r /Users/strigg/Desktop/gigas_blastP_uniprot/ /Volumes/web-1/metacarcinus/Cgigas/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 1916362\r\n",
      "-rwx------  1 strigg  staff    40M Feb  4 14:46 \u001b[31mall_giga-uniprot-blastP-out.nopipe.annotations.tab\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   3.4M Feb  4 14:46 \u001b[31mall_giga-uniprot-blastP-out.nopipe.placeholder.tab\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   3.4M Feb  4 14:46 \u001b[31mall_giga-uniprot-blastP-out.nopipe.tab\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   3.3M Feb  4 14:46 \u001b[31mgiga-uniprot-blastP-out.nopipe.tab\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   3.4M Feb  4 14:46 \u001b[31mgiga-uniprot-blastP-out.tab\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff    90K Feb  4 14:46 \u001b[31mquery_seq_names_NOTin_Blast_out\u001b[m\u001b[m*\r\n",
      "-rwx------@ 1 strigg  staff   321M Feb  5 16:27 \u001b[31muniprot-reviewed%3Ayes.tab\u001b[m\u001b[m*\r\n",
      "-rwx------@ 1 strigg  staff   6.7K Feb  4 14:46 \u001b[31muniprot_annotations_first10.org\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   263M Feb  5 16:27 \u001b[31muniprot_swprot2019_01.fasta\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   102M Feb  4 14:46 \u001b[31muniprot_swprot2019_db.phr\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   4.3M Feb  4 14:46 \u001b[31muniprot_swprot2019_db.pin\u001b[m\u001b[m*\r\n",
      "-rwx------  1 strigg  staff   192M Feb  5 16:27 \u001b[31muniprot_swprot2019_db.psq\u001b[m\u001b[m*\r\n"
     ]
    }
   ],
   "source": [
    "ls -lh /Volumes/web-1/metacarcinus/Cgigas/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}