---
title: "12-Peve-RNAseq-kallisto"
author: "Kathleen Durkin"
date: "2024-01-25"
always_allow_html: true
output: 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true
---

```{r setup, include=FALSE}
library(knitr)
knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = FALSE,        # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  comment = ""         # Prevents appending '##' to beginning of lines in code output
)
```

**Description:**

Use kallisto to quantify *P. evermanni* RNAseq transcript abundances

**Inputs:**

Trimmed RNAseq reads (e.g. `*.fastq.gz`).  *P. evermanni* reads found [here](https://gannet.fish.washington.edu/Atumefaciens/20230519-E5_coral-fastqc-fastp-multiqc-RNAseq/P_evermanni/trimmed/).

Coding sequence gff file (e.g. `*.gff`) and scaffold genome (e.g. `*.fa`), found [here](https://www.genoscope.cns.fr/corals/data/).

**Outputs:**

kalisto counts matrix (`kallisto.isoform.counts.matrix`)


---


# Create a Bash variables file

This allows usage of Bash variables (e.g. paths to common directories) across R Markdown chunks.
```{r save-bash-variables-to-rvars-file, engine='bash', eval=TRUE}
{
echo "#### Assign Variables ####"
echo ""

echo "# Data directories"
echo 'export Peve_dir=/home/shared/8TB_HDD_02/shedurkin/deep-dive/E-Peve'
echo 'export output_dir_top=${Peve_dir}/output/12-Peve-RNAseq-kallisto'
echo 'export trimmed_reads_dir=${output_dir_top}/trimmed-reads'
echo 'export kallisto_output_dir=${output_dir_top}/kallisto'
echo ""

echo "# Input/Output files"
echo 'export transcriptome_dir=${Peve_dir}/data'
echo 'export transcriptome_gff_name="Porites_evermanni_v1.annot.gff"'
echo 'export transcriptome_gff=${transcriptome_dir}/${transcriptome_gff_name}'
echo 'export transcriptome_gff_filtered_name="Porites_evermanni_v1_CDS.annot.gff"'
echo 'export transcriptome_gff_filtered=${transcriptome_dir}/${transcriptome_gff_filtered_name}'
echo 'export transcriptome_bed_name="Porites_evermanni_v1_CDS.annot.bed"'
echo 'export transcriptome_bed=${transcriptome_dir}/${transcriptome_bed_name}'
echo 'export genome_fasta_name="Porites_evermanni_v1.fa"'
echo 'export genome_fasta=${transcriptome_dir}/${genome_fasta_name}'
echo 'export transcriptome_fasta_name="Porites_evermanni_CDS.fasta"'
echo 'export transcriptome_fasta=${transcriptome_dir}/${transcriptome_fasta_name}'
echo 'export kallisto_index_name="Peve_kallisto_index.idx"'

echo "# External data URLs"
echo 'export trimmed_reads_url="https://gannet.fish.washington.edu/Atumefaciens/20230519-E5_coral-fastqc-fastp-multiqc-RNAseq/P_evermanni/trimmed/"'
echo 'export transcriptome_url="https://www.genoscope.cns.fr/corals/data/Porites_evermanni_v1.annot.gff"'
echo 'export genome_url="https://www.genoscope.cns.fr/corals/data/Porites_evermanni_v1.fa"'
echo ""

echo "# Set filename patterns"
echo "export fastq_pattern='*.fastq.gz'"
echo "export R1_fastq_pattern='*_R1.fastq.gz'"
echo "export R2_fastq_pattern='*_R2.fastq.gz'"
echo ""

echo "# Paths to programs"
echo 'export kallisto=/home/shared/kallisto_linux-v0.50.1/kallisto'
echo 'export trinity_abund_to_matrix=/home/shared/trinityrnaseq-v2.12.0/util/abundance_estimates_to_matrix.pl'
echo 'export bedtools=/home/shared/bedtools2/bin/bedtools'
echo 'export bedops=/home/shared/bedops_linux_x86_64-v2.4.41/bin'
echo ""

echo "# Set number of CPUs to use"
echo 'export threads=20'
echo ""

echo "# Programs associative array"
echo "declare -A programs_array"
echo "programs_array=("
echo '[kallisto]="${kallisto}" \'
echo '[trinity_abund_to_matrix]="${trinity_abund_to_matrix}" \'
echo '[bedtools]="${bedtools}" \'
echo '[bedops]="${bedops}" \'
echo ")"

} > .bashvars

cat .bashvars
```


# Download trimmed RNAseq reads

Reads are downloaded from: https://gannet.fish.washington.edu/Atumefaciens/20230519-E5_coral-fastqc-fastp-multiqc-RNAseq/P_evermanni/trimmed/

The `--cut-dirs 4` command cuts the preceding directory structure (i.e. `Atumefaciens/20230519-E5_coral-fastqc-fastp-multiqc-RNAseq/P_evermanni/trimmed/`) so that we just end up with the reads.

```{bash download-trimmed-reads, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${trimmed_reads_dir} \
--recursive \
--no-check-certificate \
--continue \
--cut-dirs 4 \
--no-host-directories \
--no-parent \
--quiet \
--accept "RNA-*.fastq.gz" ${trimmed_reads_url}
```

```{r view-reads-directory, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

ls -lh "${trimmed_reads_dir}"
```

# FastQC/MultiQC on trimmed reads

Already performed, can view multiqc report at https://gannet.fish.washington.edu/Atumefaciens/20230519-E5_coral-fastqc-fastp-multiqc-RNAseq/P_evermanni/trimmed/multiqc_report.html

# Retrieve the reference transcriptome and genome

Provided by https://www.genoscope.cns.fr/corals/genomes.html

Download the gene CDS (coding sequence) gff file
```{r download-transcriptome-gff, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${transcriptome_dir} \
--no-check-certificate \
--continue \
--no-host-directories \
--no-directories \
--no-parent \
--quiet \
--execute robots=off \
--accept "${transcriptome_gff_name}" ${transcriptome_url}
```
Note that this is a CDS (coding sequence) gff file, not a FASTA, so can't input directly into kallisto. We'll need the reference genome as well to convert gff to FASTA

Download the genome scaffolds FASTA file
```{r download-genome-fasta, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${transcriptome_dir} \
--no-check-certificate \
--continue \
--no-host-directories \
--no-directories \
--no-parent \
--quiet \
--execute robots=off \
--accept "${genome_fasta_name}" ${genome_url}
```

```{r view-data-directory, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

ls -lh "${transcriptome_dir}"
```

## Verify transcriptome/genome FastA MD5 checksum

No checksum file(s) provided with download, so skipping this step

# Convert gff to a genes FASTA

## Extract only CDS from gff

We only want the sequences classified as "CDS"
```{r extract-CDS-from-transcriptome-gff, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

# Extract only the CDS sequence lines from the gff
grep -w 'CDS' ${transcriptome_gff} > ${transcriptome_gff_filtered}

head -n 5 ${transcriptome_gff_filtered}
```

## Convert gff to bed file

To extract FASTAs for each of the CDS sequences we just extracted we'll be using bedtools, which takes bed files as input, so we need to convert our CDS gff file to a bed file. (Note that bedtools does accept gff files, but since gff and bed files have slightly different coordinate systems we're going to convert bed just in case)
```{r convert-gff-to-bed, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

# Ensure bedops can find its dependencies when running
export PATH=/home/shared/bedops_linux_x86_64-v2.4.41/bin:$PATH

${bedops}/gff2bed \
--do-not-sort \
< ${transcriptome_gff_filtered} \
> ${transcriptome_bed}

head -n 3 ${transcriptome_gff_filtered}
echo ""
head -n 3 ${transcriptome_bed}
```

## Generate transcriptome fasta

The below script will take as input a bed file containing information on CDS sequences, where multiple CDS sequences may originate from the same parent mRNA. The script will extract FASTAs for each sequence, and concatenate and label by parent. This should output a gene FASTA that we can use for kallisto pseudoalignment and abundance quantification!
**Warning**: This script will take a while to run -- for our bed file of 231,320 CDS sequences, the script took ~4hours to output a complete gene fasta.

```{r generate-transcriptome-fasta, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

# Navigate to correct directory and make output file
cd ${transcriptome_dir}
echo > ${transcriptome_fasta_name}

# Helper list for processing all parent IDs
processed_ids=()

######################################################

# Helper function to concatenate and format several bedtools output sequences 
# into a single, appropriately named contig
concatenate_helper() {
    local input_bedtools_fastas="$1"
    local parent_ID="$2"
    local reference_name=""
    local positions=""
    local concatenated_sequences=""

    # Read the input line by line
    while IFS= read -r line; do
        # Check if the line starts with ">"
        if [[ "$line" == ">"* ]]; then
            # Extract reference name and position from the line
            reference_position="${line:1}"  # Remove ">"
            reference_name=$(echo "$reference_position" | cut -d: -f1)
            position=$(echo "$reference_position" | cut -d: -f2)

            # Append position to the positions variable
            positions+="$position,"
        else
            # Concatenate sequences
            concatenated_sequences+="$line"
        fi
    done <<< "$input_bedtools_fastas"

    # Remove trailing comma from positions
    positions="${positions%,}"

    # Output the reformatted result
    echo ">$parent_ID $reference_name:$positions"
    echo "$concatenated_sequences"
}

######################################################

# Process your input bed file
while IFS= read -r line; do

    # pull the parent ID number for the current line of the bed
    parentID=$(echo "$line" | grep -o 'Parent=Peve_[0-9]\+')
    
    # Only continue if you haven't already processed the CDS sequences associated with this parent ID
    if [[ ! " ${processed_ids[@]} " =~ " $parentID " ]]; then
 
        # Add the current parentID to the processed list
        processed_ids+=("$parentID")

        # Create temporary files to store intermediate results
        temp_CDS_bed_file=$(mktemp)
        temp_bedtools_fasta_file=$(mktemp)

        # Grab all of the CDS sequences with the same parent ID and write to temporary file
        grep "$parentID" ${transcriptome_bed} > "$temp_CDS_bed_file"

        # Use bedtools to extract corresponding FASTAs and write to temporary file
        ${programs_array[bedtools]} getfasta -fi ${genome_fasta} -bed "$temp_CDS_bed_file" -fo "$temp_bedtools_fasta_file"

        # Use our helper function to concatenate and format all of these CDS fastas into a single contig
        concatenated_fasta=$(concatenate_helper "$(cat "$temp_bedtools_fasta_file")" "$parentID")
 
        # Add the concatenated CDS fasta to our output file on a new line
        echo "$concatenated_fasta" >> ${transcriptome_fasta}

        # Remove the temporary files
        rm "$temp_CDS_bed_file" "$temp_bedtools_fasta_file"
    fi
done < ${transcriptome_bed}

# The output file ends up having a blank first line before the data, so delete that unwanted empty first line
sed -i '1{/^$/d}' ${transcriptome_fasta}

head -n 4 ${transcriptome_fasta}
```

## Check transcriptome fasta

Let's do a quick check to see whether the output file contains all the CDS sequences we want and is grouping them appropriately
```{r check-transcriptome-fasta-head, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

# Output the first five fasta names in our transcriptome fasta (i.e. first five odd lines)
sed -nu '1~2p' ${transcriptome_fasta} | head -n 5

echo ""

# Output the first twenty CDS sequences listed in the CDS gff
head -n 20 ${transcriptome_bed}

```

```{r check-transcriptome-fasta-tail, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

# Output the first five fasta names in our transcriptome fasta (i.e. first five odd lines)
sed -nu '1~2p' ${transcriptome_fasta} | tail -n 5

echo ""

# Output the first twenty CDS sequences listed in the CDS gff
tail -n 20 ${transcriptome_bed}

```

It looks like each grouped/concatenated FASTA in our output contains the correct number of original gff sequences from the correct parent and coordinates, and the output also contains sequences for all of the parents listed in the original gff!

# Align to reference transcriptome (Kallisto pseudoalignment)

## Building Index

```{r kallisto-indexing, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

cd "${kallisto_output_dir}"

${programs_array[kallisto]} index \
--threads=${threads} \
--index="${kallisto_index_name}" \
"${transcriptome_fasta}"

ls -lh ${kallisto_output_dir}
```
Note that, when building an index, kallisto warns us that it "replaced 43264 non-ACGUT characters in the input sequence with pseudorandom nucleotides." This high number of identified "non-ACGUT" characters is related to the type of reference sequences we used to build the index. We obtained a coding sequence (CDS) gff file for P.evermanni and the associated scaffold genome fasta, filtered the gff to retain only coding sequences, and then used bedtools to extract the fasta sequences of every CDS from the scaffold fasta. Notably, scaffolds are basically fragments of known DNA sequences "stitched" together by stretches of Ns to approximate the full sequence structure without complete sequence data. This means some of our mRNA sequences contain long, relatively-meaningless stretches of Ns. I'm not sure how/to what extent this will interfere with the kallisto pseudoallignment process, since it differs from standard alignment of a full read to reference, but we'll continue for now

## Sample Quantification

Kallisto can run quantification on either single- or paired-end reads. The default option is paired-end, which requires the input of an even number of paired fastq files (e.g., pairA_R1.fastq, pairA_R2.fastq). 
To use single-end mode, include the --single flag, as well as -l (--fragment-length=DOUBLE, estimated avg. fragment length) and -s (--sd=DOUBLE, estimates stand. dev. of fragment length), and a number of fastq files.
Again, gzipped files are acceptable.

Kallisto quant is rather finicky about how you input sets of paired reads, and you can only input a single pair at a time. To circumvent, I'll create symlinks to each of the input files with simplified names, create a quantification function, and apply it iteratively to each pair using a loop.

```{r rename-trimmed-reads, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

# Create sym links to each of the trimmed read files with simplified names

for file in "${trimmed_reads_dir}"/*.fastp-trim.20230519.fastq.gz; do
    # Extract sample ID and read number from the file name
    sample_id=$(echo "$file" | grep -oP 'RNA-POR-\K\d+')
    read_number=$(echo "$file" | grep -oP '_R\K\d+')

    # Create the shortened name
    shortened_name="sample${sample_id}_R${read_number}.fastq.gz"

    # Create symbolic link
    ln -s "$file" "${trimmed_reads_dir}/${shortened_name}"

done

ls -lh ${trimmed_reads_dir}
```

```{r kallisto-quantification, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

# Function to run kallisto quant. Takes two (paired) reads as input, outputs to sample-associated directory
run_kallisto_quant() {
    source .bashvars  # Source .bashvars inside the function to make its variables accessible
    local R1_fastq=${1}
    local R2_fastq=${2}
    
    cd ${kallisto_output_dir}
    sample_num=$(basename "${R1_fastq}" "_R1.fastq.gz")
    mkdir kallisto_quant_${sample_num}

    ${programs_array[kallisto]} quant \
        --threads=${threads} \
        --index="${kallisto_output_dir}/${kallisto_index_name}" \
        --output-dir="${kallisto_output_dir}/kallisto_quant_${sample_num}" \
        --bootstrap-samples=100 \
        ${trimmed_reads_dir}/${R1_fastq} ${trimmed_reads_dir}/${R2_fastq}
}



# Iteratively apply run_kallisto_quant on each pair of input reads
for file_r1 in "${trimmed_reads_dir}"/*_R1.fastq.gz; do
    # Extract the sample name from the file name
    sample_name=$(basename "${file_r1}" "_R1.fastq.gz")

    # Form the file names (function takes input file names, not paths)
    file_r1_name="${sample_name}_R1.fastq.gz"
    file_r2_name="${sample_name}_R2.fastq.gz"

    # Check that the sample hasn't already been quantified
    if [ ! -d "${kallisto_output_dir}/kallisto_quant_${sample_name}" ]; then
    
        # Check if the corresponding R2 file exists
        if [ -e "${trimmed_reads_dir}/${file_r2}" ]; then
            # Run kallisto quant on the file pair
            run_kallisto_quant "${file_r1_name}" "${file_r2_name}" 

            echo "Processed sample: ${sample_name}"
        fi
    else
        echo "Sample already processed: ${sample_name}"
    fi
done

ls -lh ${kallisto_output_dir}
```


## Trinity Matrix with Kallisto Output

```{r kallisto-trinity-matrix, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

cd ${kallisto_output_dir}

${programs_array[trinity_abund_to_matrix]} \
--est_method 'kallisto' \
--gene_trans_map 'none' \
--out_prefix 'kallisto' \
--name_sample_by_basedir ${kallisto_output_dir}/kallisto_quant_*/abundance.tsv

ls -lh ${kallisto_output_dir}
```


# Summary

