data_processing/AnVIL/mergeSRAMetaReadstats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create a table with SRA IDs, readstats QC data, and submitter metadata\n",
    "\n",
    "**Below are the steps taken in this notebook:**\n",
    "1. Installs, import Statements & Global Variable Definitions\n",
    "2. Read in tables as dataframes\n",
    "3. Merge dataframes\n",
    "4. Write merged dataframe to data table and file\n",
    "\n",
    "**Important Note**  \n",
    "This may not be the best way to aggregate the data, because SRA submissions are done after QC, which means the workspace needs to be kept and run with newly uploaded info. It might be easier to just output the Readstats table and merge outside Terra/Anvil\n"
   ]
  },
  {
   "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 "
   ]
  },
  {
   "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"
   ]
  },
  {
   "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",
    "# table filenames (expected in the workspace bucket, see below)\n",
    "submitter_metadata = '201112_UW_HPRC_PacBio_HiFi_Metadata_Submission_v0.2_kmmod.txt'\n",
    "sra_metadata = 'metadata-8323615-processed-ok.tsv'\n",
    "\n",
    "# submission id\n",
    "submission_id = 'UW_HPRC_HiFi_Y1'\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": [
    "# Read In tables as dataframes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Readstats table\n",
    "Find the individual readstats output files in the table and concatenate their contents.\n",
    "Add the hifi filename and sample name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# readstats\n",
    "readstats_df = table_to_dataframe(\"readstats\", workspace=WORKSPACE, workspace_namespace=PROJECT)\n",
    "qc_list = []\n",
    "\n",
    "for index, row in readstats_df.iterrows():\n",
    "\n",
    "        sample_readstats_fp = row['ReadStatsReport']\n",
    "        sample_readstats_fn = os.path.basename(sample_readstats_fp)\n",
    "        df = pd.read_csv(sample_readstats_fp, header=None, sep='\\t')\n",
    "        df = df[df[0]=='sample.fastq']\n",
    "        # add sample and file name\n",
    "        #df = df.append({0:'undef', 1:'sample', 2:row['sample']}, ignore_index=True)\n",
    "        df = df.append({0:'undef', 1:'filename', 2:os.path.basename(row['hifi'])}, ignore_index=True)\n",
    "        # remove all but the variables (total_bp, quartile_25 etc)\n",
    "        df2 = df.drop(df.columns[[0, 1]], axis=1)\n",
    "        # make rownames ('total_bp', 'quartile_25' etc)\n",
    "        df2.index = df.iloc[:,1]\n",
    "        # remove the now meaningless 'sample.fastq' filename\n",
    "        df2 = df2.drop(index=['file'])\n",
    "        qc_list.append(df2)\n",
    "\n",
    "# merge\n",
    "readstats_df = pd.concat(qc_list, axis=1).transpose()\n",
    "readstats_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Submitter metadata table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_fp = os.path.join(bucket, submitter_metadata)\n",
    "metadata_df = pd.read_csv(meta_fp, sep='\\t')\n",
    "metadata_df = metadata_df[['filename', 'sample_ID', 'instrument_model', 'shear_method', 'size_selection', 'ccs_algorithm', \n",
    "    'polymerase_version', 'seq_plate_chemistry_version', 'generator_facility', 'generator_contact']]\n",
    "metadata_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SRA table\n",
    "Must download after finishing submission.\n",
    "The downloaded table needs to be split in one row per file (instead of filename, filename2 <...> filenameN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Read in SRA file and split by filename\n",
    "\n",
    "sra_file = os.path.join(bucket + sra_metadata)\n",
    "sample_df = pd.read_csv(sra_file, sep='\\t')\n",
    "\n",
    "fnames = [x for x in sample_df.columns if x.startswith('filename')]\n",
    "sra_df = pd.melt(sample_df, id_vars =['sample_name', 'accession'], value_vars = fnames, value_name='file')\n",
    "sra_df= sra_df.dropna()\n",
    "sra_df = sra_df[['file', 'accession', 'sample_name']]\n",
    "sra_df.columns = ['filename', 'accession', 'sample']\n",
    "sra_df.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge all dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sra_df.shape, readstats_df.shape, metadata_df.shape)\n",
    "# two samples of HG02572 were not uploaded to SRA with the rest of this batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(readstats_df['filename']) - set(sra_df['filename'])\n",
    "#m54329U_201103_231616.ccs.bam (HG002 sample) is not present in the submitter metadata or on SRA, but it is in the submissions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df = pd.merge(\n",
    "    sra_df,\n",
    "    readstats_df,\n",
    "    on='filename')\n",
    "merged_df = pd.merge(\n",
    "    merged_df,\n",
    "    metadata_df,\n",
    "    on='filename')\n",
    "merged_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write output files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create table\n",
    "#dataframe_to_table(\"post_sra_metadata\", merged_df, WORKSPACE, PROJECT)\n",
    "outf = os.path.join(bucket, submission_id + '_post_sra_metadata.tsv')\n",
    "merged_df.to_csv(outf, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df.shape"
   ]
  },
  {
   "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
}