{ "metadata": { "name": "", "signature": "sha256:3ccdae9211882f40cbfde52a1055f818726b0cb2dec0d59bf882ee5fb2a87116" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## IPython Notebook and data supplemental to the manuscript: \"_Up in arms: Immune and nervous system response to sea star wasting disease_\"\n", "\n", "\n", "\n", "---\n", "\n", "The code in this IPython notebook is meant to serve as a platform for anyone to explore the data associated with this research effort. In order to explore data interactively, one should download the GitHub Repository https://github.com/sr320/eimd-sswd locally, navigate to local copy of eimd-sswd and launch IPython.\n", "\n", "\n", "---\n", "To execute the IPython Notebook in its entirety you will need: \n", "\n", "* IPython - [install instructions](http://ipython.org/install.html) \n", "\n", "\n", "Note this notebook was originally developed on MacOSx and thus certain assumptions are made regarding \"built\" in software including `perl`.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Examining _Pycnopodia helianthoides_ coelomocyte _de novo_ transcriptome" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#Look at first 2 lines\n", "!head -2 ./data/Phel_transcriptome.fasta" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ ">Phel_contig_1\r\n", "CAAATATATGAACGGTTGATTGTCAACGATTAGTACATGTTTTCATTGTTCCCCACGCCCGCCCCCCCCCACTCAAACATTTAAAGTGTGAAATATTATTTATCCACAAATTTCCTTAAACCTGCAAACTTGTCTGCTGTCTCTTATTGGAAGTTATGAAAAAGAACAACGGGTTTTCTTTAAAGGGTCTGCGTGCGATTTTCAACCTTTTGAGTAATAGCAGTTATTTTGATAACCGATTTTTTTCAAAGCTCAACAGCTTTTTAAAATAAGGAATCCTATAATGGCCAAACGAATACTATAAAAATAAGGGTTCTCTTAATTGTATAAAACGTATAATTTTATCAATTTTGGGACCGTGTAATTTTTTAAAGACCACAAGAATGTTACATACAACAAATAGACGAAACTCGTAGCTTTGGAAACTACGTCATGGGCGTTTGGTCAAAAGCTGGAGAGAAAGAGAGGTGGGGTGCCAGACTTAAGTAGTCACGTGATCTGACCAACGCACATCGGAAGCTCGATCGGATGAAATCTTCTCTATCGTTCTTGCGTCTATACGTGCTACGAAGAGCTGACAGAAGTTTGGACTTGTTTACTTCTTGCACCTGTTGATGGAACGGCCACGGACCTTGTCGCACGCACACCTGGAGCCAGTGCTCGGATCGACGCAACGGATGTACTGTCTTCCCCTTCCGCGTTTCTCAAGTAGGTACTCAAAGTCGTCCGCGTCGAAGTTGGCCTCGGCGTCCCTCTTCTCCAGCTCCTCCATGTCCTCCTCTGTGTAGTACGGGGTGACGAGCACCACCAGGGCGGCCACAATGGCCAGTGCTAGAAGACACTTCGTATTCATTCTGCTGGTGGTTGGATGTGCGCAAACAAGACAGGAGAGACTTATTAGAATC\r\n" ] } ], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "#Run perl script cont_fasta.pl - Perl script\n", "#Author: Joseph Fass (modified from script by Brad Sickler) \n", "#last revised: November 2010 - http://bioinformatics.ucdavis.edu.\n", "#-i can be modified to change bin size\n", "!perl ./scripts/count_fasta.pl -i 1000 ./data/Phel_transcriptome.fasta" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\r\n", "0:999 \t15505\r\n", "1000:1999 \t8559\r\n", "2000:2999 \t3003\r\n", "3000:3999 \t1240\r\n", "4000:4999 \t580\r\n", "5000:5999 \t312\r\n", "6000:6999 \t129\r\n", "7000:7999 \t72\r\n", "8000:8999 \t31\r\n", "9000:9999 \t19\r\n", "10000:10999 \t9\r\n", "11000:11999 \t7\r\n", "12000:12999 \t2\r\n", "13000:13999 \t4\r\n", "14000:14999 \t3\r\n", "15000:15999 \t0\r\n", "16000:16999 \t1\r\n", "\r\n", "Total length of sequence:\t40747496 bp\r\n", "Total number of sequences:\t29476\r\n", "N25 stats:\t\t\t25% of total sequence length is contained in the 2260 sequences >= 3085 bp\r\n", "N50 stats:\t\t\t50% of total sequence length is contained in the 6715 sequences >= 1757 bp\r\n", "N75 stats:\t\t\t75% of total sequence length is contained in the 14612 sequences >= 959 bp\r\n", "Total GC count:\t\t\t16459121 bp\r\n", "GC %:\t\t\t\t40.39 %\r\n", "\r\n" ] } ], "prompt_number": 7 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Annotation of transcriptome" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to annoate the P. hel transcriptome `blastx` was used to compare sequences to the UniProtKB/Swiss_Prot database. \n", "The specific code is \n", "\n", "```\n", "blastx \\\n", "-query Phel_transcriptome.fasta \\\n", "-db /databaselocation/uniprot_sprot \\\n", "-out Phel_uniprot_sprot.tab \\\n", "-evalue 1E-05 \\\n", "-max_target_seqs 1 \\\n", "-max_hsps 1 \\\n", "-outfmt 6 \\\n", "-num_threads 16\n", "```\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# blastx output\n", "!head -4 ./wd/Phel_uniprot_sprot.tab" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Phel_contig_100\tsp|Q16513|PKN2_HUMAN\t81.33\t332\t61\t1\t7935\t6940\t653\t983\t5e-162\t 537\r\n", "Phel_contig_1000\tsp|Q8R4U2|PDIA1_CRIGR\t53.62\t442\t201\t2\t199\t1512\t31\t472\t5e-146\t 464\r\n", "Phel_contig_10006\tsp|Q9Y2H9|MAST1_HUMAN\t70.93\t289\t82\t1\t861\t1\t434\t722\t6e-132\t 415\r\n", "Phel_contig_10021\tsp|Q96MU7|YTDC1_HUMAN\t60.85\t212\t82\t1\t1115\t1750\t294\t504\t1e-73\t 258\r\n" ] } ], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "#determining how many sequences have blast hits\n", "!wc -l ./wd/Phel_uniprot_sprot.tab" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 10513 ./wd/Phel_uniprot_sprot.tab\r\n" ] } ], "prompt_number": 13 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Joining with UniProt information in SQLShare" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "!python {sqls}fetchdata.py \\\n", "-s \"SELECT * FROM [{usr}].[Phel_uniprot_sprot_sql.tab]phel left join [sr320@washington.edu].[uniprot-reviewed_wGO_010714]des on phel.Column3=des.Entry\" \\\n", "-f tsv \\\n", "-o ./wd/Phel_uniprot_info.tab\n", "```\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "!head -3 ./wd/Phel_uniprot_ID.tab" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Phel_contig_100 Q16513 Serine/threonine-protein kinase N2 (EC 2.7.11.13) (PKN gamma) (Protein kinase C-like 2) (Protein-kinase C-related kinase 2) 5E-162 Homo sapiens (Human)\r\n", "Phel_contig_1000 Q8R4U2 Protein disulfide-isomerase (PDI) (EC 5.3.4.1) (Prolyl 4-hydroxylase subunit beta) (p58) 5E-146 Cricetulus griseus (Chinese hamster) (Cricetulus barabensis griseus)\r\n", "Phel_contig_10006 Q9Y2H9 Microtubule-associated serine/threonine-protein kinase 1 (EC 2.7.11.1) (Syntrophin-associated serine/threonine-protein kinase) 6E-132 Homo sapiens (Human)\r\n" ] } ], "prompt_number": 36 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Differential Expression Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Count data from CLC..." ] }, { "cell_type": "code", "collapsed": false, "input": [ "!head ./data/Phel_countdata.txt" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Feature ID\tTreated_FHL - Total gene reads\tTreated_PH - Total gene reads\tTreated_L - Total gene reads\tControl_FHL - Total gene reads\tControl_DB - Total gene reads\tControl_PH - Total gene reads\r\n", "Phel_contig_1\t168\t37\t8\t89\t28\t38\r\n", "Phel_contig_10\t9518\t2752\t839\t22\t42\t180\r\n", "Phel_contig_100\t260\t565\t413\t616\t1234\t6104\r\n", "Phel_contig_1000\t2043\t3842\t3070\t4311\t8527\t31946\r\n", "Phel_contig_10000\t9\t12\t13\t32\t21\t211\r\n", "Phel_contig_10001\t44\t225\t89\t90\t54\t365\r\n", "Phel_contig_10002\t38\t61\t80\t185\t478\t1267\r\n", "Phel_contig_10003\t9\t29\t20\t17\t29\t186\r\n", "Phel_contig_10004\t8\t25\t6\t4\t19\t92\r\n" ] } ], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "!wc -l ./data/Phel_countdata.txt" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 29476 ./data/Phel_countdata.txt\r\n" ] } ], "prompt_number": 51 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### DESeq2 Code\n", "\n", "```R\n", "%%R\n", "data <- read.table(\"./data/Phel_countdata.txt\", header = T, sep = \"\\t\")\n", "rownames(data) <- data$Feature\n", "data <- data[,-1]\n", "\n", "# Build Objects\n", "# Specify which columns are in which groups\n", "deseq2.colData <- data.frame(condition=factor(c(rep(\"Treated\", 3), rep(\"Control\", 3))), \n", " type=factor(rep(\"single-read\", 6)))\n", "rownames(deseq2.colData) <- colnames(data)\n", "deseq2.dds <- DESeqDataSetFromMatrix(countData = data,\n", " colData = deseq2.colData, \n", " design = ~ condition)\n", "# Count number of hits with adjusted p-value less then 0.05\n", "dim(deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ])\n", "\n", "tmp <- deseq2.res\n", "# The main plot\n", "plot(tmp$baseMean, tmp$log2FoldChange, pch=20, cex=0.45, ylim=c(-3, 3), log=\"x\", col=\"darkgray\",\n", " main=\"DEG Virus Exposure (pval <= 0.05)\",\n", " xlab=\"mean of normalized counts\",\n", " ylab=\"Log2 Fold Change\")\n", "# Getting the significant points and plotting them again so they're a different color\n", "tmp.sig <- deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ]\n", "points(tmp.sig$baseMean, tmp.sig$log2FoldChange, pch=20, cex=0.45, col=\"red\")\n", "# 2 FC lines\n", "abline(h=c(-1,1), col=\"blue\")\n", "\n", "# Run Analysis\n", "deseq2.dds <- DESeq(deseq2.dds)\n", "deseq2.res <- results(deseq2.dds)\n", "deseq2.res <- deseq2.res[order(rownames(deseq2.res)), ]\n", "\n", "# Count number of hits with adjusted p-value less then 0.05\n", "dim(deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ])\n", "\n", "tmp <- deseq2.res\n", "# The main plot\n", "plot(tmp$baseMean, tmp$log2FoldChange, pch=20, cex=0.45, ylim=c(-3, 3), log=\"x\", col=\"darkgray\",\n", " main=\"DEG Virus Exposure (pval <= 0.05)\",\n", " xlab=\"mean of normalized counts\",\n", " ylab=\"Log2 Fold Change\")\n", "# Getting the significant points and plotting them again so they're a different color\n", "tmp.sig <- deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ]\n", "points(tmp.sig$baseMean, tmp.sig$log2FoldChange, pch=20, cex=0.45, col=\"red\")\n", "# 2 FC lines\n", "abline(h=c(-1,1), col=\"blue\")\n", "\n", "write.table(tmp.sig, \"./wd/Phel_DEGlist.tab\", row.names = T)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![plot](./precompiled_wd/deseq2.png)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "!head -3 ./wd/Phel_DEGlist.tab" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\"baseMean\" \"log2FoldChange\" \"lfcSE\" \"stat\" \"pvalue\" \"padj\"\r\n", "\"Phel_contig_10\" 5492.70960072414 7.62981963801218 0.791209689281713 9.64323331901912 5.250750567205e-22 2.77316470200627e-19\r\n", "\"Phel_contig_10001\" 136.224412089297 1.6018041057117 0.596619525301856 2.68480000700827 0.00725732180374353 0.0435439308224612\r\n" ] } ], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }