3_Biology_Example_AWS_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "51b77691",
   "metadata": {},
   "source": [
    "## Setting up AWS Athena\n",
    "\n",
    "\n",
    "### Step 1: Go to AWS Glue Service to Create and Run a \"Crawler\"\n",
    "\n",
    "* A crawler is an automated service in AWS that looks through a set of classifiers prioritized by the schema to create metadata tables.  \n",
    "\n",
    "* To get the most recent version of the metadata, which is updated by SRA daily, it is necessary to rerun the crawler.\n",
    "\n",
    "* This differs from BigQuery because BigQuery database and table names are static and the metadata tables are recreated daily by SRA using automated scripts.\n",
    "\n",
    "<img src=\"./img/crawler.png\" align=\"left\" width=\"700\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b428867",
   "metadata": {},
   "source": [
    "### Step 2: Confirm the Name of the Database and the Tables\n",
    "\n",
    "* When the crawler is created, the user sets the database name. In this workshop, \"sra\" was the name set by the crawler.\n",
    "\n",
    "* The level of the tables that show up in the database depends on which datastore is targeted by the crawler. In this case, all the SRA tables are visible because the crawler was set to run on the highest level directory **s3://sra-pub-metadata-us-east-1**\n",
    "\n",
    "* Note that the tables are all in one database called **\"sra\"**. This differs from BigQuery which has two sets of tables in its schema: **\"sra\"** and **\"sra_tax_analysis_tool\"**\n",
    "\n",
    "\n",
    "#### Athena Tables\n",
    "\n",
    "<img src=\"./img/table.png\" align=\"center\" width=\"700\"/>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#### BigQuery Tables\n",
    "\n",
    "<img src=\"./img/BigQuery.png\" align=\"center\" width=\"350\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6202f97",
   "metadata": {},
   "source": [
    "The Query workspace for Athena looks very similar to the BigQuery workspace.  The syntax is almost identical, but before starting, make note of how the Database and Tables are structured. The column names are exactly the same between the two providers.\n",
    "\n",
    "<img src=\"./img/athena.png\" align=\"center\" width=\"800\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bf24831",
   "metadata": {},
   "source": [
    "### Step 3: Create a \"Staging\" Bucket\n",
    "\n",
    "* Athena differs from BigQuery in that it requires a user to create a bucket to hold all the metadata and query outputs\n",
    "\n",
    "* This is true for both command line interface as well as direct console queries.\n",
    "\n",
    "### Run the following command to create your bucket for this workshop\n",
    "\n",
    "This script also applies credentials that allow this notebook (which is on GCP) to access AWS to create buckets and run Athena. It also creates a connection using sqlalchemy to the Athena database.\n",
    "\n",
    "*The script is unique to this workshop, so make sure to add your own AWS credentials when working outside of this event. Refer to the setup_AWS_credentials.ipynb notebook for tips and tricks on how to do this.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90a16a2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run Setup_AWS_Credentials.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59a464ff",
   "metadata": {},
   "source": [
    "### Test that the Credentials were Loaded Correctly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fea03593",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%sql\n",
    "\n",
    "SELECT acc, bioproject \n",
    "FROM sra.metadata\n",
    "WHERE organism = 'Homo sapiens'\n",
    "LIMIT 5\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "639ec2d8",
   "metadata": {},
   "source": [
    "# Challenge : To find RNA-Seq data from estrogen receptor negative and progesterone positive breast cancer cell lines.\n",
    "\n",
    "\n",
    "Hormone receptor positive (HR+) breast cancer cells may have either estrogen receptors (ER+) or progesterone receptors (PR+) or both.\n",
    "\n",
    "We have already found one search strategy to find a dataset that had ER+ breast cancer cell lines.\n",
    "\n",
    "This case study will look for HR+ breast cancer cell lines that are ER- and PR+\n",
    "\n",
    "This will requre combining search strategies to narrow our choices.\n",
    "\n",
    "These searchers are taking place on Athena, so the syntax is mostly the same with some subtle differences.\n",
    "\n",
    "Some of they syntax, such as capitalization or lowercase, is dicated by the tool we are using, so if a query in this notebook does not work exactly the same in Athena, try that first.\n",
    "\n",
    "### Start with a more general search that uses an \"or\" to specify receptor types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36e390ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%sql\n",
    "\n",
    "\n",
    "SELECT acc, bioproject, extracted.k, extracted.v\n",
    "FROM sra.metadata, unnest(attributes) as t(extracted)\n",
    "WHERE assay_type = 'RNA-Seq'\n",
    "    AND consent = 'public'\n",
    "    AND (extracted.v LIKE '%ER-%' or extracted.v LIKE '%HR+%')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "477e1c21",
   "metadata": {},
   "source": [
    "### A lot of the metadata we are picking up is not related to samples. So the next search specifies that the key value needs to end in \"sam\" by using a wildcard \"%\" before \"sam\" and not after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "641d445c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%sql\n",
    "\n",
    "\n",
    "SELECT acc, bioproject, extracted.k, extracted.v\n",
    "FROM sra.metadata, unnest(attributes) as t(extracted)\n",
    "WHERE assay_type = 'RNA-Seq'\n",
    "    AND consent = 'public'\n",
    "    AND (extracted.v LIKE '%ER-%' or extracted.v LIKE '%HR+%')\n",
    "    AND extracted.k LIKE '%sam'\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8719806",
   "metadata": {},
   "source": [
    "### Adding the word \"receptor\" in front of the wildcard here may limit the samples to a reasonable set for the next steps. \n",
    "\n",
    "There are several strategies to find interesting data sets, this is one that seems to work just within this case study. A good practice would be to spend some time trying different methods to get the best data set for your research question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c86acbf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%sql\n",
    "\n",
    "\n",
    "SELECT acc, bioproject, extracted.k, extracted.v\n",
    "FROM sra.metadata, unnest(attributes) as t(extracted)\n",
    "WHERE assay_type = 'RNA-Seq'\n",
    "    AND consent = 'public'\n",
    "    AND (extracted.v LIKE '%ER-%' or extracted.v LIKE '%HR+%')\n",
    "    AND extracted.k LIKE 'receptor%sam'\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11b981d4",
   "metadata": {},
   "source": [
    "#### This looks like a good set, let's pull out the accession list for downstream analysis.\n",
    "\n",
    "The syntax is slightly different because we are using sqlmagic - a different program than bigquery.\n",
    "\n",
    "The output of a query can be saved as a results object which is easily made into a dataframe for downstream analysis.  We are creating a results object called \"accessions\" and then saving it into a text file called \"cell_line_ER_HR.txt\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c43097c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%sql accessions <<\n",
    "\n",
    "\n",
    "SELECT acc\n",
    "FROM sra.metadata, unnest(attributes) as t(extracted)\n",
    "WHERE assay_type = 'RNA-Seq'\n",
    "    AND consent = 'public'\n",
    "    AND (extracted.v LIKE '%ER-%' or extracted.v LIKE '%HR+%')\n",
    "    AND extracted.k LIKE 'receptor%sam'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63bf8762",
   "metadata": {},
   "source": [
    "#### The output is different in SQL magic compared to BigQuery magic in the other notebook.\n",
    "\n",
    "It is not a dataframe yet, it would need to have accessions.DataFrame() to show the index as in the BigQuery results object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cad3e217",
   "metadata": {},
   "outputs": [],
   "source": [
    "accessions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2a6aa96",
   "metadata": {},
   "outputs": [],
   "source": [
    "accessions.DataFrame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ebcaf5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "file = open(\"cell_line_ER_HR.txt\", \"w\")\n",
    "file.write(accessions.DataFrame().to_string(index=False, header=False))\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "527f201f",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "## Example Bioinformatic Workflow\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c45a4d6",
   "metadata": {},
   "source": [
    "### We may not have time to go over this in detail, the code is here for your information ###\n",
    "\n",
    "<img src=\"./img/Workflow_demo.png\" align=\"left\" width=\"700\"/>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### We are going to skip down to the final analysis for this demo  ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1979c00",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "cat cell_line_ER_HR.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d336e521",
   "metadata": {},
   "source": [
    "### Configure SRA Toolkit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c59b4430",
   "metadata": {},
   "outputs": [],
   "source": [
    "GUID = !uuidgen\n",
    "GUID = str(GUID).strip('[]')\n",
    "!vdb-config --set LIBS/GUID={GUID}\n",
    "!vdb-config --report-cloud-identity yes\n",
    "!vdb-config --set /repository/user/main/public/cache-disabled=true\n",
    "#!cat $HOME/.ncbi/user-settings.mkfg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e069b782",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "cat cell_line_ER_HR.txt | xargs -P 12 -n 1 fastq-dump -X 1000000 -O cell_line_fastq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1db0ce8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "ls cell_line_fastq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "406de847",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "datasets download genome accession GCF_000001405.39 --chromosomes 17 --include-gtf --exclude-protein --exclude-rna\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e03cdf7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "unzip -o ncbi_dataset.zip -d GRCh38.p13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f1fb800",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "makeblastdb -in GRCh38.p13/ncbi_dataset/data/GCF_000001405.39/chr17.fna -dbtype nucl -parse_seqids -out GRCh38.p13/GRCh38.p13\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebd439ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "mkdir cell_line_bams"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee932939",
   "metadata": {},
   "source": [
    "### Making alignments\n",
    "\n",
    "This script can be automated, but for the purposes of this demo, we are choosing only two of the fastq files to align to Chromosome 17 - a hot spot of activity for HER genes.  \n",
    "\n",
    "* SRR12060806 is HR+\n",
    "\n",
    "* SRR12060824 is ER-/PR+\n",
    "\n",
    "Doing a little more digging into bioproject PRJNA640820: \"62 breast cell lines composed of 27 Triple negative breast cancers (TNBC), 5 Non malignant (NM) lines, 14 Hormone receptor positive (HR+) and 16 Her2 amplified lines (Her2amp) were profiled.\" This was a relatively recent dataset produced on a NextSeq here in Boston.\n",
    "\n",
    "<img src=\"./img/citation.png\" align=\"left\" width=\"700\"/>\n",
    "\n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152908\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c865101",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "magicblast -query cell_line_fastq/SRR12060806.fastq -db GRCh38.p13/GRCh38.p13 \\\n",
    "-splice T -no_unaligned -infmt fastq -num_threads 8 | \\\n",
    "samtools view -bS | samtools sort -@ 8 -o cell_line_bams/SRR12060806.bam -O BAM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbdde7f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "magicblast -query cell_line_fastq/SRR12060824.fastq -db GRCh38.p13/GRCh38.p13 \\\n",
    "-splice T -no_unaligned -infmt fastq -num_threads 8 | \\\n",
    "samtools view -bS | samtools sort -@ 8 -o cell_line_bams/SRR12060824.bam -O BAM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6af51010",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "samtools index cell_line_bams/SRR12060824.bam\n",
    "samtools index cell_line_bams/SRR12060806.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96103ece",
   "metadata": {},
   "source": [
    "## Load the Data in the NCBI Genome Data Viewer\n",
    "\n",
    "These bam files with their indexes need to be loaded somewhere they can be retrieved by the NCBI Genome Data Viewer.\n",
    "\n",
    "If you want to load your own remote files into GDV Viewer, it is possible using your Github repository for files that are **less than 1 GB**. Just use the following file format based on the location of your file in the repository.\n",
    "\n",
    "\n",
    "Follow this format:  https:///<i></i>github.com/YOUR-GITHUB-ACCOUNT/YOUR-GITHUB-REPOSITORY/raw/master/PATH-TO-BAM\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b577dff1",
   "metadata": {},
   "source": [
    "## See the Two Bams from this Demo\n",
    "\n",
    "The two aligned bams have been pre-loaded at the following link: \n",
    "    \n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?cfg=NCID_1_53497080_130.14.18.128_9146_1631700092_3532671673\n",
    "\n",
    "which will be good for the next 90 days.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26005f7",
   "metadata": {},
   "source": [
    "\n",
    "Here are the links to the two bams if you want to load them again, make sure to zoom to Chromosome 17, that is where the reads were aligned. The screenshot below indicates where to enter these URL's to have the tracks show up.\n",
    "\n",
    "\n",
    "<img src=\"./img/add_remote_view.png\" align=\"left\" width=\"600\"/>\n",
    "\n",
    "\n",
    "https://github.com/ncbi/ASHG-Workshop-2021/raw/main/bams/SRR12060806.bam\n",
    "\n",
    "\n",
    "https://github.com/ncbi/ASHG-Workshop-2021/raw/main/bams/SRR12060824.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2599d1fc",
   "metadata": {},
   "source": [
    "## What did we find???\n",
    "\n",
    "Actually, starting with a general idea and targeting a specific question but being general in our search terms, this demo shows the power of using SQL to rapidly search millions of records in SRA. \n",
    "\n",
    "The speed and specificity of using SQL for will help researchers quickly locate sequences of interest, bring these sequences directly into a workflow and compare the outcomes within the command line interface.\n",
    "\n",
    "\n",
    "\n",
    "***\n",
    "\n",
    "#### The GDV link leads to the TP53 gene, an important protein in the development of cancer.\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"./img/tp53_gdv.png\" align=\"left\" width=\"600\"/>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18813379",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### The mismatch in the pileups help narrow the scope of where interesting SNPs might be.\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"./img/pileup_view.png\" align=\"left\" width=\"600\"/>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1095909d",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### One of the cell lines has a SNP that indicates that cisplatin will be a less effective drug treatment.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"./img/drug_response_variant.png\" align=\"left\" width=\"600\"/>\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b83894a",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### Reading more into the records shows the importance of this variant for many cancer drugs.\n",
    "\n",
    "<img src=\"./img/clinvar.png\" align=\"left\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7157a9e9",
   "metadata": {},
   "source": [
    "The cell line with the SNP may be a good candidate for drug therapy studies to emulate individuals whose cancers do not respond to cisplatin. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dbc61d4",
   "metadata": {},
   "source": [
    "# Thank you!!!\n",
    "\n",
    "## Please Take our Survey!\n",
    "\n",
    "https://nlmenterprise.co1.qualtrics.com/jfe/form/SV_bx6evo2ghfCjkrk\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e8e4e21",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}