---
title: "Cicada Data"
subtitle: "Data analysis reproduction concerning the cicada's genome."

author: 
    - name: "Eric Mossotti"
      corresponding: true
      email: ecmossotti@gmail.com
      roles:
        - Data analysis

keywords: 
    - Genome Assembly
    - ETL Pipeline
    - Python
    - R
    - Bash

date: 07/15/2024
date-modified: last-modified
date-format: "MMM D, YYYY"

bibliography: references.bib

repo: https://github.com/ericMossotti/Cicada_Data
source: index.qmd


abstract-title: "Objective"
abstract: |
      To reproduce DNA Zoo's summary table on the 17-year cicada.
      
description-meta: "Reproducible genomics analysis."

license: CC BY-SA
funding: "The author(s) received no specific funding for this work."

csl: csl/apa.csl
citation-location: margin
citations-hover: true
link-citations: true
citeproc: true

toc-expand: true
toc: true
number-sections: true
toc-float: true

link-external-newwindow: true
smooth-scroll: true
fig-responsive: true
echo: true
callout-appearance: simple
code-fold: true
code-overflow: wrap

filters:
     - nutshell

# For favicon support on most platforms
# Generated the image at https://svg.io/
# Converted to favicon package at https://realfavicongenerator.net/
include-in-header:
  text: |
    <link rel="apple-touch-icon" sizes="180x180" href="favicon/apple-touch-icon.png">
    <link rel="icon" type="image/png" sizes="32x32" href="favicon/favicon-32x32.png">
    <link rel="icon" type="image/png" sizes="16x16" href="favicon/favicon-16x16.png">
    <link rel="manifest" href="favicon/site.webmanifest">
    <link rel="mask-icon" href="favicon/safari-pinned-tab.svg" color="#5bbad5">
    <meta name="msapplication-TileColor" content="#da532c">
    <meta name="theme-color" content="#ffffff">
---

## Introduction {.unnumbered}

[The full repo for this project with instructions on how to replicate and run with near verbatim output.](https://github.com/ericMossotti/Cicada_Data)

### Problem {.unnumbered}

The steps involved in reproducing data can be unclear.

### Purpose {.unnumbered}

To elaborate on the objective stated at the top of this document, I seek to supplement DNA Zoo's report with an accessible data analysis (DA) pipeline. To accomplish this, I independently reproduce the original article's table while documenting every data processing step. Although there's nothing wrong with the original works, things can always be taken further. [@dnazoo], [@magicica]

### Stakeholders {.unnumbered}

This might be of interest to the original authors of the article. More generally, the spirit of this work could transfer to other domains of data intensive research and analytics.

### Source {.unnumbered}

All data used within this report was freely available from a public database hosted by DNA Zoo. [@dnazoo.s]

### Pipeline {.unnumbered .column-screen-right}

```{mermaid}
%%| code-fold: true
%%| code-summary: Create a Mermaid diagram from code.

flowchart TB
    A((1)):::circle --> B((2)):::circle
    B --> C((3)):::circle
    C --> D((4)):::circle

    subgraph Extract ["1. Extract"]
        direction LR
        A1["directorize.py"] --> A2["importer.py"]
    end
    subgraph Transform ["2. Transform"]
      direction TB
      B1{"decompress.py"} -.->|.fasta| B2["assembly_stats"] -.->|summary_stats.txt| B3["assemblyFramer.py"]
      B1 -.->|.fasta| B4["assemblyDictionary.py"]
    end
    subgraph Load ["3. Load"]
        direction TB
        C1{"strint.py"}
    end
    subgraph Present ["4. Present"]
        direction TB
        D1["DNA Zoo's Table, Reproduced"]
    end
	
    A ~~~ Extract
    B ~~~ Transform
    C ~~~ Load
    D ~~~ Present
	
    A2 -.->|fasta.gz| B1
    B3 -.->|dataframe| C1
    B4 -.->|dict| C1
    C1 -.->|strings| Present
    C1 -.->|strings| Present
    
```

## Extract

This would be the data extraction phase of the DA pipeline.

### Create Project Directory

```{python}
#| label: code_directorize
#| code-summary: A function used further below to establish the project's DA directory.

# Automates creation of the DA pipeline directory needed for this project to avoid confusion for those replicating this project, if desired.

import os

def directorize(base_path, structure):
    # Nested directories
    for dir_name, subdirs in structure.items():
        dir_path = os.path.join(base_path, dir_name)
        os.makedirs(dir_path, exist_ok = True)

        for subdir in subdirs:
            subdir_path = os.path.join(dir_path, subdir)
            os.makedirs(subdir_path, exist_ok = True)
```

```{python}
#| label: run_directorize
#| code-summary: If the directory folders already exist, it shouldn't cause any issues to allow this code to evaluate during render.

# Define the directory structure
structure = {
    "01_Extract/": ["data/", "scripts/"],
    "02_Transform/": ["data/", "scripts/"],
    "03_Load/": ["data/", "scripts/"],
    "04_Present/": ["data/", "scripts/"]
}

# Create the analysis folder structure in a preferred base directory
# "" = project's working directory
directorize("", structure)
```

### Download to Local Machine

If you want the exact verbatim output, then delete the cloned folder structure (01_Extract through the 04_Present and run this code).

```{python}
#| label: code_importer
#| code-summary: Importing all scripts and data needed for the project from the web to make things easier for replicating this project. If you want the exact verbatim output, then delete the cloned folder structure (01_Extract through the 04_Present and run this code). Otherwise, it's ok to re-run this with the scripts already there. 

# ---- Import data from the web with wget
import os
import sys
import wget

def importer (fileMap):
    # Download from URL to path and notify when complete
    for url, file_path in fileMap.items():
        # Checking file existence
        if not os.path.exists(file_path):
            wget.download(url, file_path)
            print(f"{file_path} written")
        else:
            print(f"{file_path} already exists.")
```

```{python}
#| label: run_importer
#| code-summary: Run `importer()`, passing in the mapped URLs keys to their local directory path destination values.

# Set the url
fasta_URL = "https://dnazoo.s3.wasabisys.com/Magicicada_septendecula/magicicada_hifiasm.asm.bp.p_ctg_HiC.fasta.gz"
# Set the desired local file path
fasta_PATH = "01_Extract/data/magicicada.fasta.gz"

# Map the url to the file path
fileMap = {
  fasta_URL: fasta_PATH,
  }

importer(fileMap)

```

::: callout-note
### The specific link used to download all data from

<https://dnazoo.s3.wasabisys.com/Magicicada_septendecula/magicicada_hifiasm.asm.bp.p_ctg_HiC.fasta.gz>
:::

## Transform

The data transformation phase of the pipeline.

### Decompress .GZ

```{r}
#| label: src_decompress.py
#| code-summary: Source the `decompress.py` script.

reticulate::source_python("02_Transform/scripts/decompress.py")
```

```{python}
#| label: code_decompress.py
#| code-summary: Transform the compressed file to its decompressed form.
#| file: "02_Transform/scripts/decompress.py"
#| eval: false 
```

```{python}
#| label: run_decompress.py
#| code-summary: Run `decompress()`, using the compressed path and the desired path to the decompressed file. 

# Set the compressed fasta.gz file variable
gzFasta = "01_Extract/data/magicicada.fasta.gz"

# Set the decompressed fasta file variable
fasta = "02_Transform/data/magicicada.fasta"

# Pass file paths to the function
decompress(gzFasta, fasta)
```

### FASTA to Text, to DataFrame

Sourcing the dataframe formatting script here, as there is some issue with knitr engine and Python with these parts. Seems like some kind of conflict with switching between dataframe and dictionary data types.

```{r}
#| label: src_formaFrame_df.py
#| code-summary: Source the `formaFrame.py` script.

reticulate::source_python("04_Present/scripts/formaFrame.py")
```

This chunk should be ran locally instead of with `quarto render`. When working with the source file, change the code-chunk language specifier from `{.bash}` back to `{bash}`. You might have to add the `{bash}` tag entirely back to the div. Not sure how else to go about accomplishing this within my current Quarto project setup. [@trizna2020]

```{.bash}
# BASH SCRIPT

# The uncompressed fasta file variable
fasta=02_Transform/data/magicicada.fasta

# The text file path variable generated by the script
summary_stats=02_Transform/data/summary_stats.txt

assembly_stats $fasta > $summary_stats
```

I need to do this for the rendering part to display the code-cell outputs how I wanted them to. However, you can use the bash code chunk further below and comment this code out if you want. This makes it so you can deactivate the bash code chunk for render but still have everything else work without a hitch if you go simply from empty directory to render with this document.

Now, transform the text file into a Python dataframe. I am opting not to blanket change data-types as output format could vary by user preference.

```{r}
#| label: src_assemblyFramer.py
#| code-summary: Source the `assemblyFramer.py` script.

# Import external python script to local library environment
reticulate::source_python("02_Transform/scripts/assemblyFramer.py")
```

::: scroll-code
```{python}
#| label: assemblyFramer_code_view
#| code-summary: Transform the text file data into a mult-indexed dataframe. Multi-indexing simplifies query syntax later on.
#| file: "02_Transform/scripts/assemblyFramer.py"
#| eval: false 
```
:::

```{python}
#| label: run_assemblyFramer
#| code-summary: Run `assemblyFramer()`, passing in the path to the text file generated by running the Bash script from earlier. 

# Set the local text file path
statsPath = "02_Transform/data/summary_stats.txt"
# Run to yield a multi-indexed dataframe
df = assemblyFramer(statsPath)
```

### FASTA to Dictionary

```{r}
#| label: src_assemblyDictionary.py
#| code-summary: Source the external `assemblyDictionary.py` script.

reticulate::source_python("02_Transform/scripts/assemblyDictionary.py")
```

::: scroll-code
```{python}
#| label: code_assemblyDictionary.py
#| code-summary: Return a stats dictionary.
#| file: "02_Transform/scripts/assemblyDictionary.py"
#| eval: false 
```
:::

```{r}
#| label: src_strint.py
#| code-summary: Source the `strint.py` script.

reticulate::source_python("03_Load/scripts/strint.py")
```

```{python}
#| label: run_assemblyDictionary.py
#| code-summary: Instead of a text file, this skips several time-consuming steps. 
#| eval: false

statsDict = assemblyDictionary("02_Transform/data/magicicada.fasta")
```

## Load

This is the data loading phase. Following completion of this stage, querying the data should be more intuitive than before.

### Assign Variables with `strint.py`

::: scroll-code
```{python}
#| label: code_strint.py
#| code-summary: This code transforms the strings in the value column of the dataframe to a more visually appealing, thousands-separated form.
#| file: "03_Load/scripts/strint.py"
#| eval: false 
```
:::

### DataFrame Input

```{python}
#| label: run_strint.py
#| code-summary: Call `strint()`, with the simpflified query syntax noted earlier. Then simply specify the literal 'Category' and 'Label' keys to return the desired value.
#| eval: true

#---- Contigs
ctig_len = strint(df, "Contigs", "total_bps")
ctig_count = strint(df, "Contigs", "sequence_count")
ctig_n50 = strint(df, "Contigs", "N50")
ctig_max = strint(df, "Contigs", "longest")

#---- Scaffolds
sfld_len = strint(df, "Scaffolds", "total_bps")
sfld_count = strint(df, "Scaffolds", "sequence_count")
sfld_n50 = strint(df, "Scaffolds", "N50")
sfld_max = strint(df, "Scaffolds", "longest")
```

:::: {.callout-tip appearance="simple"}
#### Multi-index dataframe query syntax

::: scroll-code
```{python}
#| eval: false
#| code-fold: show
#| echo: true
#| code-summary: An example of the syntax used to query the multi-indexed dataframe. This code is not evaluated, but only for illustrative purposes.

"""

strint(dataframe, category, label)

category options
----------------
  Contigs
  Scaffolds


label options
----------------
  L10
  L20
  L30
  L40
  L50
  N10
  N20
  N30
  N40
  N50
  gc_content
  longest
  mean
  median
  sequence_count
  shortest
  total_bps
  
"""

ctig_len = strint(df, "Contigs", "N50")

# -->> Looking inside the strint(dataframe, category, label) function ---->>

# Which then finds the desired value or 'Value'
stat = df.loc[(category, label), "Value"]

```
:::
::::

### Dictionary Input

```{python}
#| label: assign_strings
#| code-summary: Assigns the dictionary value as string objects.
#| eval: false

#---- Contigs
ctig_len = strint(statsDict, "Contig Stats", "total_bps")
ctig_count = strint(statsDict, "Contig Stats", "sequence_count")
ctig_n50 = strint(statsDict, "Contig Stats", "N50")
ctig_max = strint(statsDict, "Contig Stats", "longest")

#---- Scaffolds
sfld_len = strint(statsDict, "Scaffold Stats", "total_bps")
sfld_count = strint(statsDict, "Scaffold Stats", "sequence_count")
sfld_n50 = strint(statsDict, "Scaffold Stats", "N50")
sfld_max = strint(statsDict, "Scaffold Stats", "longest")
```

:::: {.callout-tip appearance="simple"}
#### Dictionary query syntax

::: scroll-code
```{python}
#| eval: false
#| code-fold: show
#| echo: true
#| code-summary: An example of the syntax used to query the dictionary This code is not evaluated, but only for illustrative purposes.

"""

strint(dataframe, category, label)

category options
----------------
  Contig Stats
  Scaffold Stats


label options
----------------
  L10
  L20
  L30
  L40
  L50
  N10
  N20
  N30
  N40
  N50
  gc_content
  longest
  mean
  median
  sequence_count
  shortest
  total_bps
  
"""

# A look inside assemblyDictionary.py where stat_output is returned
stat_output = {'Contig Stats': contig_stats,
               'Scaffold Stats': scaffold_stats}

ctig_len = statsDict["Contig Stats"]["total_bps"]
```
:::
::::

## Present

### The Pandas Table

This is a slightly formatted view of the Pandas table designed to be more easily queried to return the desired statistic. If, however, you'd like to treat the Styler object as the unchanged, dataframe object, use the `forma_df.data` syntax.

[:The original dataframe output:](#NutFrame)

```{python}
#| label: styled_dataframe
#| code-summary: The Pandas dataframe as a Styler object for better viewing.

# Display with Style 
forma_df = df.loc[:].style.pipe(formaFrame)

forma_df
```

### DNA Zoo's Table, Reproduced

```{r}
#| label: load_reticulate
#| code-summary: To make the variable calls simpler in the table below.

library(reticulate)
```

Importing the library makes it simpler for inserting the values into the table below. For example, I would have had to type `{r} reticulate::py$ctig_n50`, but now, I only need to type `{r} py$ctig_n50` into the individual cells. I needed to convert the Python into R objects, as the knitr engine used in rendering this document does not seem to display output from execution of inline Python code directly.

::: {#ZooTable}
|  |  |  |  |
|------------------|------------------|------------------|-------------------|
| **Contig length (bp)** | **Number of contigs** | **Contig N50 (bp)** | **Longest contig (bp)** |
| `{r} py$ctig_len` | `{r} py$ctig_count` | `{r} py$ctig_n50` | `{r} py$ctig_max` |
| **Scaffold length (bp)** | **Number of scaffolds** | **Scaffold N50 (bp)** | **Longest scaffold (bp)** |
| `{r} py$sfld_len` | `{r} py$sfld_count` | `{r} py$sfld_n50` | `{r} py$ sfld_max` |
:::

## :x NutFrame {.unnumbered .unlisted}

```{python}
#| label: nutshFrame
#| code-summary: The dataframe object on display, in all its glory.

# The un-styled dataframe output
forma_df.data
```
