3 Xenopus laevis samples.

1 Introduction

We received via a google drive link 3 samples of RNASeq data. They are dissociated_embryo_x6, whole_embryo_x1, and whole_embryo_x3.

I downloaded a copy of the xenopus laevis genome revision 9.1 along with the tropicalis 10.1 as genbank flat files from NCBI and converted them to the requisite gff/genome/CDS/amino acid entries.

Najib suggested we should apply our own trimming methods, though I fully expected the provided trimmed files are perfect.

2 Preprocessing

I will use the following ~ 6 shell blocks to invoke the preprocessing tools which will handle the data from trimming->quantification. Depending on those results, I may add some further analyses.

At the bottom of this document I am going to add some blocks into which I will copy/paste the actual scripts which were invoked on the computers in our cluster. I will just copy/paste one example each for an arbitrarily chosen sample.

2.1 Some re-used settings

I will put a few parameters as environment variables to be reused across invocations.

species="xenopus_laevis_v10.1"
xl_type="gene"
xl_tag="gene"
stranded="no"

2.2 Trimming and QC

I will use a few tools for quick trimming/QC metrics.

2.2.1 Trimomatic

The following for loop is how I do pretty much everything.

cd preprocessing
module add cyoa
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd "${start}/${i}"
    input=$(/bin/ls unprocessed/*.fastq.gz | tr '\n' ':')
    cyoa --method trim --input "${input}"
done
cd "${start}"

2.2.2 Fastqc

cd preprocessing
module add cyoa
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd "${start}/${i}"
    input=$(/bin/ls unprocessed/*.fastq.gz | tr '\n' ':')
    cyoa --method fastqc --input "${input}"
done
cd "${start}"

2.2.3 Fastp

A new tool which performs trimming and QC metrics at the same time is fastp.

cd preprocessing
module add cyoa
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd "${start}/${i}"
    input=$(/bin/ls unprocessed/*.fastq.gz | tr '\n' ':')
    cyoa --method fastp --input "${input}"
done
cd "${start}"

2.3 Gene/Tx quantification

2.3.1 Transcript quantification via salmon

cd preprocessing
module add cyoa
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd "${start}/${i}"
    input=$(/bin/ls outputs/01trimomatic/*-trimmed.fastq.xz | tr '\n' ':')
    cyoa --method salmon --input "${input}" --species "${species}"
done
cd "${start}"

2.3.2 Gene quantification via hisat2

cd preprocessing
module add cyoa
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd "${start}/${i}"
    input=$(/bin/ls outputs/01trimomatic/*-trimmed.fastq.xz | tr '\n' ':')
    cyoa --method hisat --species "${species}" --input "${input}" --gff_type "${xl_type}" \
         --gff_tag "${xl_tag}" --stranded "${stranded}"
done
cd "${start}"

2.4 rRNA Query

cd preprocessing
module add cyoa
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd "${start}/${i}"
    input=$(/bin/ls unprocessed/*.fastq.gz | tr '\n' ':')
    cyoa --method hisat --libtype rRNA --species "${species}"_ --input "${input}" \
         --stranded "${stranded}" --gff_type rRNA --gff_tag gene
done
cd "${start}"

3 Load gene annotations

I will load my annotations from a combination of ensembl and the NCBI gff file; I presume the gene IDs will match up, since the NCBI assembly and ensembl have the same source.

xl_mRNA_annot <- load_gff_annotations("reference/xenopus_laevis_v10.1.gff",
                                      type = "mRNA", id_col = "gene")

## Returning a df with 42 columns and 1007786 rows.

xl_CDS_annot <- load_gff_annotations("reference/xenopus_laevis_v10.1.gff",
                                      type = "CDS", id_col = "gene")

## Returning a df with 42 columns and 956598 rows.

## xl_ens_annot <- load_biomart_annotations(species = "xlaevis", archive = FALSE)

4 Load the salmon transcript counts

This will likely be tricky because salmon has some peculiar rules about gene/transcript IDs which I strongly suspect will not match up with the IDs produced when reading from the gff annotations or biomart. Let us find out! Here are the first few entries of the salmon quant file:

XP_041442323.1_LOC108704873_GeneID:108704873_pepsin_A_isoform_X1 1143 853.109 0.487714 9.000 XP_041442326.1_LOC108704873_GeneID:108704873_pepsin_A_isoform_X2 1035 743.000 0.000000 0.000 XP_018097054.2_LOC108704861_GeneID:108704861_peptidase_inhibitor_16 1656 1364.000 0.000000 0.000 XP_041436388.1_dok1.L_GeneID:444335_docking_protein_1_L_homeolog_isoform_X1 2127 1811.679 6.149833 241.000 NP_001085908.1_dok1.L_GeneID:444335_docking_protein_1_L_homeolog 1782 1476.220 0.000000 0.000 NP_001084645.1_mrps26.L_GeneID:414604_uncharacterized_protein_LOC414604 615 317.758 155.818955 1071.000 XP_018097055.1_LOC108704862_GeneID:108704862_meiosis_1_arrest_protein 1515 1223.000 0.000000 0.000 XP_018097050.2_dqx1.L_GeneID:108704857_LOW_QUALITY_PROTEIN:_ATP-dependent_RNA_helicase_DQX1 2197 1852.495 2.944778 118.000 XP_041442332.1_LOC108704874_GeneID:108704874_matrix_metalloproteinase-21 1857 1565.000 0.000000 0.000

So, yeah, that is going to be annoying.

I usually get all of my genomes/annotations from ensembl, I am now remembering why this is. I do not think it should prove difficult to cross reference my ncbi and ensembl IDs, but we shall see.

The IDs from salmon are:

“{protein_id}{gene}{db_xref}_{product}” with the caveat that all spaces are replaced by underscore.

salmon_annot <- xl_CDS_annot
rownames(salmon_annot) <- make.names(salmon_annot[["protein_id"]], unique = TRUE)
salmon_tx_map <- salmon_annot[, c("protein_id", "gene")]
salmon_tx_map[["protein_id"]] <- gsub(x = salmon_tx_map[["protein_id"]], pattern = "\\.[[:digit:]]+",
                                      replacement = "")

5 Create SummarizedExperiment

I will create one version by tx and one by gene.

tx_se <- create_se("sample_sheets/all_samples.xlsx", gene_info = salmon_annot,
                   file_column = "salmon_quant")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 3 rows(samples) and 4 columns(metadata fields).

## Matched 60687 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final summarized experiment has 60905 rows and 4 columns.

salmon_gene_annot <- salmon_annot
rownames(salmon_gene_annot) <- make.names(salmon_gene_annot[["gene"]], unique = TRUE)
gene_se <- create_se("sample_sheets/all_samples.xlsx", gene_info = salmon_gene_annot,
                  file_column = "salmon_quant", tx_gene_map = salmon_tx_map)

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 3 rows(samples) and 4 columns(metadata fields).

## In some cases, (notably salmon) the format of the IDs used by this can be tricky.
## It is likely to require the transcript ID followed by a '.' and the ensembl column:
## 'transcript_version', which is explicitly different than the gene version column.
## If this is not correctly performed, very few genes will be observed

## reading in files with read_tsv

## 1 2 3 
## removing duplicated transcript rows from tx2gene
## transcripts missing from tx2gene: 218
## summarizing abundance
## summarizing counts
## summarizing length
## Matched 34144 annotations and counts.
## Bringing together the count matrix and gene information.
## Some annotations were lost in merging, setting them to 'undefined'.
## The final summarized experiment has 34255 rows and 4 columns.

6 Look at the 3 samples

dim(assay(tx_se))

## [1] 60905     3

plot_libsize(tx_se)

## Library sizes of 3 samples, 
## ranging from 6,088,658 to 14,516,405.

dim(assay(gene_se))

## [1] 34255     3

plot_libsize(gene_se)

## Library sizes of 3 samples, 
## ranging from 6,088,627 to 14,516,358.

plot_nonzero(tx_se)

## The following samples have less than 39588.25 genes.

## [1] "dssctd_m_6" "whl_mbry_1" "whl_mbry_3"

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## A non-zero genes plot of 3 samples.
## These samples have an average 9.817 CPM coverage and 30216 genes observed, ranging from 28553 to
## 32847.

plot_nonzero(gene_se)

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## A non-zero genes plot of 3 samples.
## These samples have an average 9.817 CPM coverage and 23806 genes observed, ranging from 22726 to
## 25651.

tx_norm <- normalize(tx_se, convert = "cpm", filter = TRUE, transform = "log2")

## Running normalize_se.

## Removing 41018 low-count genes (19887 remaining).

## transform_counts: Found 2183 values equal to 0, adding 1 to the matrix.

plot_corheat(tx_norm)

## A heatmap of pairwise sample correlations ranging from: 
## 0.736448006840574 to 0.844499630181708.

plot_pca(tx_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by disociated, whole_one, whole_three
## Shapes are defined by six, one, three.

write_xlsx(excel = "excel/l2_cpm_tx.xlsx", data = assay(tx_norm))

## Deleting the file excel/l2_cpm_tx.xlsx before writing the tables.

## write_xlsx() wrote excel/l2_cpm_tx.xlsx.
## The cursor is on sheet first, row: 19890 column: 5.

7 Example Invocations

7.1 Trimomatic

#!/usr/bin/bash
#SBATCH --export=ALL --mail-type=NONE --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1
#SBATCH --job-name=01trim_XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1 --nice=10
#SBATCH --output=/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt.sbatch
#SBATCH --account=cbcb
#SBATCH --partition=cbcb
#SBATCH --qos=medium
#SBATCH --nodes=1 --ntasks=1 --cpus-per-task=4
#SBATCH --time=36:00:00
#SBATCH --mem=6G
startdir=$(pwd)
set -o errexit
set -o errtrace
set -o pipefail
echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/scripts/01trim_XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt
function get_sigterm {
  cd "${startdir}"
  echo "A SIGTERM was sent to 01trim_XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1: ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt
  exit 1
}
trap get_sigterm SIGTERM
function get_sigerr {
  cd "${startdir}"
  echo "A SIGERR was sent to 01trim_XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1: ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt
  exit 1
}
trap get_sigerr ERR
mod=$( { type -t module || true; } )
if [[ -z "${mod}" ]]; then
  module() {
  # shellcheck disable=SC2086
    { eval "$(/usr/bin/modulecmd bash $*)" || true; }
  }
  export -f module
fi
module add trimomatic 2>/dev/null 1>&2
## This call to trimomatic removes illumina and epicentre adapters from unprocessed/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1.fastq.gz:unprocessed/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2.fastq.gz.
## It also performs a sliding window removal of anything with quality <25;
## cutadapt provides an alternative to this tool.
## The original sequence data is recompressed and saved in the sequences/ directory.
mkdir -p outputs/01trimomatic
## Note that trimomatic prints all output and errors to STDERR, so send both to output
trimomatic PE \
  -threads 1 \
  -phred33 \
  unprocessed/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1.fastq.gz unprocessed/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2.fastq.gz \
  outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed_paired.fastq outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed_unpaired.fastq \
  outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed_paired.fastq outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed_unpaired.fastq \
   ILLUMINACLIP:/fs/cbcb-software/RedHat-8-x86_64/local/cyoa/202506/prefix/lib/perl5/auto/share/dist/Bio-Adventure/genome/adapters.fa:2:20:10:2:keepBothReads  \
  SLIDINGWINDOW:4:20 MINLEN:50 \
  1>outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003-trimomatic.stdout \
  2>outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003-trimomatic.stderr
excepted=$( { grep "Exception" "outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003-trimomatic.stdout" || test $? = 1; } )
## The following is in case the illumina clipping fails, which it does if this has already been run I think.
if [[ "${excepted}" != "" ]]; then
  trimomatic PE \
    -threads 1 \
    -phred33 \
    unprocessed/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1.fastq.gz unprocessed/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2.fastq.gz \
    outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed_paired.fastq outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed_unpaired.fastq \
    outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed_paired.fastq outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed_unpaired.fastq \
     SLIDINGWINDOW:4:25 MINLEN:50 \
    1>outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003-trimomatic.stdout \
    2>outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003-trimomatic.stderr
fi
sleep 10
mv outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed_paired.fastq outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed.fastq
mv outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed_paired.fastq outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed.fastq

## Recompress the unpaired reads, this should not take long.
xz -9e -f outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed_unpaired.fastq
xz -9e -f outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed_unpaired.fastq
## Recompress the paired reads.
xz -9e -f outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed.fastq
xz -9e -f outputs/01trimomatic/XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed.fastq
ln=$(cd outputs/01trimomatic && ln -sf XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1-trimmed.fastq.xz r1_trimmed.fastq.xz)
ln=$(cd outputs/01trimomatic && ln -sf XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R2-trimmed.fastq.xz r2_trimmed.fastq.xz)



## The following lines give status codes and some logging
## This might not work because it is a little circular.
cd ${startdir}
minutes_used=$(( SECONDS / 60 ))
echo "  $(hostname) Finished ${SLURM_JOBID} 01trim_XNEME_20250520_A00904_IL22367-001_N5UD-F09_L003_R1.sh at $(date), it took ${minutes_used} minutes." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt
if [[ -x "$(command -v sstat)" && -n "${SLURM_JOBID}" ]]; then
  echo "  walltime used by ${SLURM_JOBID} was: ${minutes_used:-null} minutes." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt
  echo "" >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x1/outputs/log.txt
fi
## Note, you can score a bunch more information by running cyoa --method checkjob from the working directory.

7.2 FastQC

#!/usr/bin/bash
#SBATCH --export=ALL --mail-type=NONE --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6
#SBATCH --job-name=01fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed --nice=10
#SBATCH --output=/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt.sbatch
#SBATCH --account=nexus
#SBATCH --partition=tron
#SBATCH --qos=medium
#SBATCH --nodes=1 --ntasks=1 --cpus-per-task=8
#SBATCH --time=3:00:00
#SBATCH --mem=20G
startdir=$(pwd)
set -o errexit
set -o errtrace
set -o pipefail
echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/scripts/01fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt
function get_sigterm {
  cd "${startdir}"
  echo "A SIGTERM was sent to 01fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed: ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt
  exit 1
}
trap get_sigterm SIGTERM
function get_sigerr {
  cd "${startdir}"
  echo "A SIGERR was sent to 01fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed: ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt
  exit 1
}
trap get_sigerr ERR
mod=$( { type -t module || true; } )
if [[ -z "${mod}" ]]; then
  module() {
  # shellcheck disable=SC2086
    { eval "$(/usr/bin/modulecmd bash $*)" || true; }
  }
  export -f module
fi
module add fastqc 2>/dev/null 1>&2
## This FastQC run is against unfiltered data and is used for
## an initial estimation of the overall sequencing quality.
mkdir -p outputs/01fastqc
which perl 2>outputs/01fastqc/fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed-unfiltered_fastqc.stderr 1>&2
which fastqc 2>>outputs/01fastqc/fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed-unfiltered_fastqc.stderr 1>&2
## Even if fastqc finishes happily using a subshell, it might exit with SIGERR
trap - ERR
finished=$(fastqc --extract \
  -o outputs/01fastqc \
   /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/unprocessed/XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1.fastq.gz  /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/unprocessed/XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R2.fastq.gz  \
  2>>outputs/01fastqc/fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed-unfiltered_fastqc.stderr \
  1>>outputs/01fastqc/fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed-unfiltered_fastqc.stdout || /bin/true)
echo "fastqc finished with $? andor ${finished}."



## The following lines give status codes and some logging
## This might not work because it is a little circular.
cd ${startdir}
minutes_used=$(( SECONDS / 60 ))
echo "  $(hostname) Finished ${SLURM_JOBID} 01fqc_XNEME_20250520_A00904_IL22369-001_N5UD-H09_L003_R1_unprocessed.sh at $(date), it took ${minutes_used} minutes." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt
if [[ -x "$(command -v sstat)" && -n "${SLURM_JOBID}" ]]; then
  echo "  walltime used by ${SLURM_JOBID} was: ${minutes_used:-null} minutes." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt
  echo "" >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/dissociated_embryo_x6/outputs/log.txt
fi
## Note, you can score a bunch more information by running cyoa --method checkjob from the working directory.

7.3 Fastp

#!/usr/bin/bash
#SBATCH --export=ALL --mail-type=NONE --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3
#SBATCH --job-name=12fastp_XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1 --nice=10
#SBATCH --output=/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt.sbatch
#SBATCH --account=cbcb
#SBATCH --partition=cbcb
#SBATCH --qos=medium
#SBATCH --nodes=1 --ntasks=1 --cpus-per-task=2
#SBATCH --time=8:00:00
#SBATCH --mem=12G
startdir=$(pwd)
set -o errexit
set -o errtrace
set -o pipefail
echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/scripts/12fastp_XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt
function get_sigterm {
  cd "${startdir}"
  echo "A SIGTERM was sent to 12fastp_XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1: ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt
  exit 1
}
trap get_sigterm SIGTERM
function get_sigerr {
  cd "${startdir}"
  echo "A SIGERR was sent to 12fastp_XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1: ${SLURM_JOBID}." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt
  exit 1
}
trap get_sigerr ERR
mod=$( { type -t module || true; } )
if [[ -z "${mod}" ]]; then
  module() {
  # shellcheck disable=SC2086
    { eval "$(/usr/bin/modulecmd bash $*)" || true; }
  }
  export -f module
fi
module add fastp 2>/dev/null 1>&2
## Run fastp on raw data


mkdir -p outputs/12fastp
trimmed=1
{
  /usr/bin/time -v -o outputs/12fastp/stdout.time -a
    fastp   -i /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/unprocessed/XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1.fastq.gz \
  -o outputs/12fastp/XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1-fastp.fastq \
  -I /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/unprocessed/XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R2.fastq.gz \
  -O outputs/12fastp/XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R2-fastp.fastq  \
      -h outputs/12fastp/fastp_report.html \
  -j outputs/12fastp/fastp_report.json  -D  -c  -y  \
      2>outputs/12fastp/stderr \
      1>outputs/12fastp/stdout

xz -9e -f outputs/12fastp/XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1-fastp.fastq
xz -9e -f outputs/12fastp/XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R2-fastp.fastq

} || {
  echo 'fastp failed.' >> outputs/12fastp/stderr
}


## The following lines give status codes and some logging
## This might not work because it is a little circular.
cd ${startdir}
minutes_used=$(( SECONDS / 60 ))
echo "  $(hostname) Finished ${SLURM_JOBID} 12fastp_XNEME_20250520_A00904_IL22368-001_N5UD-G09_L003_R1.sh at $(date), it took ${minutes_used} minutes." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt
if [[ -x "$(command -v sstat)" && -n "${SLURM_JOBID}" ]]; then
  echo "  walltime used by ${SLURM_JOBID} was: ${minutes_used:-null} minutes." >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt
  echo "" >> /fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/nemes_2025/preprocessing/whole_embryo_x3/outputs/log.txt
fi
## Note, you can score a bunch more information by running cyoa --method checkjob from the working directory.

pander::pander(sessionInfo())

R version 4.4.2 (2024-10-31)

Platform: x86_64-pc-linux-gnu

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: tidyr(v.1.3.1), hpgltools(v.1.0), glue(v.1.8.0), forcats(v.1.0.0) and dplyr(v.1.1.4)

loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), jsonlite(v.2.0.0), tximport(v.1.34.0), magrittr(v.2.0.3), rmarkdown(v.2.29), farver(v.2.1.2), fs(v.1.6.6), BiocIO(v.1.16.0), zlibbioc(v.1.52.0), vctrs(v.0.6.5), memoise(v.2.0.1), Rsamtools(v.2.22.0), DelayedMatrixStats(v.1.28.1), RCurl(v.1.98-1.17), htmltools(v.0.5.8.1), S4Arrays(v.1.6.0), curl(v.6.4.0), broom(v.1.0.8), Rhdf5lib(v.1.28.0), SparseArray(v.1.6.2), rhdf5(v.2.50.2), sass(v.0.4.10), bslib(v.0.9.0), htmlwidgets(v.1.6.4), bsseq(v.1.42.0), plyr(v.1.8.9), plotly(v.4.10.4), cachem(v.1.1.0), GenomicAlignments(v.1.42.0), mime(v.0.13), lifecycle(v.1.0.4), iterators(v.1.0.14), pkgconfig(v.2.0.3), Matrix(v.1.7-3), R6(v.2.6.1), fastmap(v.1.2.0), GenomeInfoDbData(v.1.2.13), MatrixGenerics(v.1.18.1), shiny(v.1.11.0), digest(v.0.6.37), colorspace(v.2.1-1), AnnotationDbi(v.1.68.0), S4Vectors(v.0.44.0), DESeq2(v.1.46.0), GenomicRanges(v.1.58.0), RSQLite(v.2.4.1), labeling(v.0.4.3), httr(v.1.4.7), abind(v.1.4-8), compiler(v.4.4.2), pander(v.0.6.6), withr(v.3.0.2), bit64(v.4.6.0-1), backports(v.1.5.0), BiocParallel(v.1.40.2), DBI(v.1.2.3), HDF5Array(v.1.34.0), R.utils(v.2.13.0), DelayedArray(v.0.32.0), corpcor(v.1.6.10), rjson(v.0.2.23), gtools(v.3.9.5), permute(v.0.9-7), tools(v.4.4.2), DSS(v.2.54.0), zip(v.2.3.3), httpuv(v.1.6.16), varhandle(v.2.0.6), restez(v.2.1.5), R.oo(v.1.27.1), restfulr(v.0.0.15), GOSemSim(v.2.32.0), rhdf5filters(v.1.18.1), promises(v.1.3.3), grid(v.4.4.2), reshape2(v.1.4.4), fgsea(v.1.32.4), generics(v.0.1.4), gtable(v.0.3.6), BSgenome(v.1.74.0), tzdb(v.0.5.0), R.methodsS3(v.1.8.2), hms(v.1.1.3), data.table(v.1.17.6), XVector(v.0.46.0), BiocGenerics(v.0.52.0), ggrepel(v.0.9.6), foreach(v.1.5.2), pillar(v.1.10.2), stringr(v.1.5.1), vroom(v.1.6.5), yulab.utils(v.0.2.0), limma(v.3.62.2), later(v.1.4.2), splines(v.4.4.2), lattice(v.0.22-7), rtracklayer(v.1.66.0), bit(v.4.6.0), annotate(v.1.84.0), tidyselect(v.1.2.1), GO.db(v.3.20.0), locfit(v.1.5-9.12), Biostrings(v.2.74.1), knitr(v.1.50.1), IRanges(v.2.40.1), edgeR(v.4.4.2), SummarizedExperiment(v.1.36.0), xfun(v.0.52), stats4(v.4.4.2), Biobase(v.2.66.0), statmod(v.1.5.0), matrixStats(v.1.5.0), stringi(v.1.8.7), UCSC.utils(v.1.2.0), lazyeval(v.0.2.2), yaml(v.2.3.10), evaluate(v.1.0.4), codetools(v.0.2-20), tibble(v.3.3.0), qvalue(v.2.38.0), graph(v.1.84.1), cli(v.3.6.5), xtable(v.1.8-4), jquerylib(v.0.1.4), dichromat(v.2.0-0.1), Rcpp(v.1.0.14), GenomeInfoDb(v.1.42.3), png(v.0.1-8), fastcluster(v.1.3.0), XML(v.3.99-0.18), parallel(v.4.4.2), readr(v.2.1.5), ggplot2(v.3.5.2), blob(v.1.2.4), DOSE(v.4.0.1), sparseMatrixStats(v.1.18.0), bitops(v.1.0-9), viridisLite(v.0.4.2), GSEABase(v.1.68.0), scales(v.1.4.0), openxlsx(v.4.2.8), purrr(v.1.0.4), crayon(v.1.5.3), rlang(v.1.1.6), cowplot(v.1.1.3), fastmatch(v.1.1-6) and KEGGREST(v.1.46.0)

message(paste0("This is hpgltools commit: ", get_git_commit()))

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset cdba1d3cf1448fd60c6a78e7b7b171e0c9a8e47b

## This is hpgltools commit: Thu Jun 5 12:57:42 2025 -0400: cdba1d3cf1448fd60c6a78e7b7b171e0c9a8e47b

message(paste0("Saving to ", savefile))

## Saving to index.rda.xz

##tmp <- sm(saveme(filename = savefile))

tmp <- loadme(filename = savefile)