data_processing/AnVIL/Generate_Terra_Tables_HiFi_QC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Data Tables For HiFi Data QC<a class=\"tocSkip\">\n",
    "\n",
    "**This notebook automatically reads in data stored in the AnVIL_HPRC workspace's bucket and generates data table so the data can be run through QC (NTSM and ReadStat's WDLs)**\n",
    "\n",
    "**Jeltje changed this in Feb 2023 to output per-file tables instead of per-sample**\n",
    "\n",
    "**Below are the steps taken in this notebook:**\n",
    "1. Import Statements & Global Variable Definitions\n",
    "2. Define Functions\n",
    "3. Read In Sample Names\n",
    "4. Create Dataframe Of Files\n",
    "5. Write data frame to data tables\n",
    "\n",
    "**Prior to running set the following variables:**\n",
    "1. sample_info_fn (line 112)\n",
    "2. subm_key (line 113)\n",
    "3. file_type_ls (line 226)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import Statements & Global Variable Definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%pip install gcsfs\n",
    "## capture CANNOT have comments above it\n",
    "## For reading CSVs stored in Google Cloud (without downloading them first)\n",
    "## May need to restart kernel after install"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%pip install --upgrade --no-cache-dir --force-reinstall terra-pandas\n",
    "%pip install --upgrade --no-cache-dir  --force-reinstall git+https://github.com/DataBiosphere/terra-notebook-utils\n",
    "## For reading/writing data tables into pandas data frames\n",
    "## May need to restart kernel after install\n",
    "## If terra_pandas module not found, try installing w/o upgrade option:\n",
    "## %pip install terra-pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from firecloud import fiss\n",
    "import pandas as pd         \n",
    "import os                 \n",
    "import subprocess       \n",
    "import re                 \n",
    "import io\n",
    "import gcsfs\n",
    "\n",
    "from typing import Any, Callable, List, Optional\n",
    "from terra_notebook_utils import table, WORKSPACE_NAME, WORKSPACE_GOOGLE_PROJECT\n",
    "from terra_pandas import dataframe_to_table, table_to_dataframe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Global Variable Declarations\n",
    "Update this with project specific info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# AnVIL_HPRC WorkspaceBucket\n",
    "anvil_hprc_bucket       = \"gs://fc-4310e737-a388-4a10-8c9e-babe06aaf0cf/\"\n",
    "\n",
    "sample_info_fn = 'UW_HPRC_HiFi_samples.tsv' # a list of IDs,e.g. HG01243\n",
    "subm_key       = '48f9717b-23c6-4045-a33f-10a3df8fc924--UW_HPRC_HiFi'\n",
    "\n",
    "# Get the Google billing project name and workspace name for current workspace\n",
    "PROJECT = os.environ['WORKSPACE_NAMESPACE']\n",
    "WORKSPACE =os.path.basename(os.path.dirname(os.getcwd()))\n",
    "bucket = os.environ['WORKSPACE_BUCKET'] + \"/\"\n",
    "\n",
    "# Verify that we've captured the environment variables\n",
    "print(\"Billing project: \" + PROJECT)\n",
    "print(\"Workspace: \" + WORKSPACE)\n",
    "print(\"Workspace storage bucket: \" + bucket)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function Definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def rtn_datatype_ls_for_subm(bucket_url, submission_key, data_type_subdir, file_type_ls):\n",
    "    '''Takes in:\n",
    "            * bucket_url (string): url of bucket to search (should be the AnVIL_HPRC bucket)\n",
    "                ex: \"gs://fc-4310e737-a388-4a10-8c9e-babe06aaf0cf/\"\n",
    "            * submission_key (string): the UUID plus the submission name to search\n",
    "                ex: \"5c68b972-8534-402f-9861-11c93558765f--UW_HPRC_HiFi_Y3\"\n",
    "            * data_type_subdir (string): name of the subfolder to search \n",
    "              (used when a submission has more than one data type.)\n",
    "                ex: \"PacBio_HiFi\" if the data is not in subfolders, pass \".\"\n",
    "            * file_type_ls (list of strings): file extensions to search for. Often a submission will \n",
    "              have more than one type of file that represents the same dataset. Be careful to not \n",
    "              include replicate data, however.\n",
    "                ex: \".hifi_reads.bam\"\n",
    "     then performs a list of the bucket, then returns a filtered list files.'''\n",
    "    \n",
    "    rtn_file_ls = []\n",
    "    \n",
    "    submission_path = str(bucket_url + \"submissions/\" + submission_key)\n",
    "    file_list_byt   = subprocess.run(['gsutil', '-u', 'human-pangenome-ucsc', 'ls', '-r', submission_path], \n",
    "                                     stdout=subprocess.PIPE)\n",
    "\n",
    "    file_list_str   = file_list_byt.stdout.decode('utf-8')\n",
    "    file_list       = file_list_str.split('\\n')  ## Pull out \"\\n\"\n",
    "   \n",
    "    ## filter out empty strings\n",
    "    file_list = [ elem for elem in file_list if elem != '']\n",
    "    \n",
    "    ## Pull from subdir type we are targeting\n",
    "    file_list = list(filter(lambda x:re.search(rf\"{data_type_subdir}\", x), file_list))\n",
    "    \n",
    "    for file_type in file_type_ls:\n",
    "    \n",
    "        ## Pull files of correct type (ex: ccs.bam)\n",
    "        file_list_by_type = list(filter(lambda x:re.search(rf\"{file_type}$\", x), file_list))\n",
    "\n",
    "        ## Add to list of files to return\n",
    "        rtn_file_ls += file_list_by_type\n",
    "\n",
    "    return rtn_file_ls    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read In Sample Names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## This file (a list of sample IDs) should be stored in the Terra workspace (in the bucket) that is being used for \n",
    "## the submission you  are wrangling...\n",
    "sample_info_fp = os.path.join(bucket + sample_info_fn) \n",
    "\n",
    "sample_df = pd.read_csv(sample_info_fp, header=None)\n",
    "\n",
    "sample_df.rename(columns = {0:'sample_id'}, inplace = True)\n",
    "sample_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Dataframe Of Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_type_subdir = \".\"\n",
    "file_type_ls     = [\".ccs.bam\"]\n",
    "\n",
    "## get a list of the files in submission (subm_key) that match the \n",
    "## data_type_subdir and file_type_ls\n",
    "file_ls  = rtn_datatype_ls_for_subm(anvil_hprc_bucket, subm_key, \n",
    "                                      data_type_subdir, file_type_ls)\n",
    "\n",
    "## Check that the number of files matches what we expect\n",
    "len(file_ls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add Sample Data To File Data Frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filedict = dict()\n",
    "\n",
    "for f in file_ls:\n",
    "    for s in sample_df.sample_id:\n",
    "        if re.search(s,f):\n",
    "            filedict[f] = s\n",
    "            break            \n",
    "print(len(filedict))\n",
    "file_df = pd.DataFrame.from_dict(filedict.items())\n",
    "file_df.columns = ['hifi', 'sample']\n",
    "file_df.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add 1000G data (for NTSM run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Read in 1000G data from another workspace\n",
    "## We will be using this to compare against out submission to make sure that\n",
    "## the data comes from the same samples\n",
    "thousand_genomes_df = table_to_dataframe(\"sample\", \n",
    "                                        workspace=\"1000G-high-coverage-2019\", \n",
    "                                        workspace_namespace=\"anvil-datastorage\")\n",
    "## Use the 1kg library name (i.e. HG00621) as the index (that matches our sample name)\n",
    "thousand_genomes_df.set_index('library_name', inplace=True)\n",
    "\n",
    "## We only need the cram file (represents 30X Ilmn dataset)\n",
    "thousand_genomes_df = thousand_genomes_df[['cram']]\n",
    "\n",
    "## name the column to be a bit more descriptive\n",
    "thousand_genomes_df.rename(columns = {'cram':'1000g_cram'}, inplace = True)\n",
    "\n",
    "## Sometimes some of the samples are not present in the 1000g file\n",
    "addSamples = list(set(file_df['sample']) - set(thousand_genomes_df.index.tolist()))\n",
    "print('missing:', addSamples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## If there are any missing samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## For the missing files we're going to get genome info from anvil_HPRC\n",
    "## the data comes from the same samples\n",
    "anvil_hprc_df = table_to_dataframe(\"sample\", \n",
    "                                        workspace=\"AnVIL_HPRC\", \n",
    "                                        workspace_namespace=\"anvil-datastorage\")\n",
    "## sample_id is the index\n",
    "anvil_hprc_df = anvil_hprc_df[['child_ilmn', 'hic']]\n",
    "\n",
    "## we want the child_ilnm value for HG002 instead of the hic value\n",
    "anvil_hprc_df.loc[['HG002']]['hic'] = anvil_hprc_df.loc[['HG002']]['child_ilmn']\n",
    "\n",
    "## now only retain the one column\n",
    "anvil_hprc_df = anvil_hprc_df[['hic']]\n",
    "\n",
    "## name the column to match 1000genomes\n",
    "anvil_hprc_df.rename(columns = {'hic':'1000g_cram'}, inplace = True)\n",
    "\n",
    "# and concatenate the two\n",
    "thousand_genomes_df = pd.concat([anvil_hprc_df.loc[addSamples], thousand_genomes_df], ignore_index=False)\n",
    "\n",
    "# check\n",
    "list(set(file_df['sample']) - set(thousand_genomes_df.index.tolist()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## merge the two dataframes\n",
    "ntsm_df = pd.merge(\n",
    "    file_df,\n",
    "    thousand_genomes_df,\n",
    "    left_on='sample',\n",
    "    right_index=True)\n",
    "ntsm_df.head()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "## Read in 1000G data from another workspace\n",
    "## instead of cram files (as before) we can use fastq file so we're not dependent\n",
    "## on the public hg18.fa file, which has been giving problems.\n",
    "## We will be using this to compare against out submission to make sure that\n",
    "## the data comes from the same samples\n",
    "t2t_reads_df = table_to_dataframe(\"participant\", \n",
    "                                 workspace_namespace=\"anvil-datastorage\", \n",
    "                                 workspace=\"AnVIL_T2T\")\n",
    "t2t_reads_df = t2t_reads_df[['read_1_fastq', 'read_2_fastq']]\n",
    "\n",
    "t2t_reads_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ntsm_df.shape, file_df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Upload To Tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create tables for running NTSM and ReadStats\n",
    "dataframe_to_table(\"ntsm\",      ntsm_df, WORKSPACE, PROJECT)\n",
    "dataframe_to_table(\"readstats\", file_df, WORKSPACE, PROJECT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}