Examining the Pseudomonas ORFeome via deep sequencing of a plasmid library.

1 Introduction

This document seeks to lay out the strategy employed for an examination of the sequencing results from a group of plasmid libraries from Pseudmonas aeruginosa PA01.

These samples are in 5 main libraries, each of which was done twice, once with a topoisomerase during isolation and once without.

The immediate goal of this experiment is to learn what ORFs exist in the library and what mutations have appeared in them.

2 Preprocessing

The preprocessing is mostly a normal DNA sequencing pipeline:

1. Trim the reads.
2. Map them against the plasmid, pDONR201.
3. Map the unmapped reads from #2 against PA01.
4. Count the reads by locus_tag to identify missing genes.
5. Use freebayes to identify variant positions.
6. Parse the freebayes output down to an easy to interpret format.

2.1 Trimming

This is my default trimomatic invocation.

cd preprocessing
start=$(pwd)
for i in $(/bin/ls); do
    cd ${i}
    cyoa --task trim --method trim --input $(/bin/ls *.gz | tr '\n' ':' | sed 's/:$//g')
    cd $start
done

Here is the script which resulted for the poola sample:

I still have an old check in it from an earlier version of trimomatic which would sometimes fail.

2.1.1 Individual trimming script

#SBATCH --export=ALL
#SBATCH --mail-type=NONE
#SBATCH --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola
#SBATCH --partition=dpart
#SBATCH --qos=workstation --nice=10
#SBATCH --nodes=1 --requeue
#SBATCH --time=24:00:00
#SBATCH --job-name=01trim_ORFeome_Pool_A_S1_R1_001
#SBATCH --mem=24G
#SBATCH --cpus-per-task=3
#SBATCH --output=outputs/log.txt.sbatch
set -o errexit
set -o errtrace
set -o pipefail
export LESS='--buffers 0'
script="$(pwd)/$0"
err() {
    echo "Error occurred:"
    awk 'NR>L-4 && NR<L+4 { printf "%-5d%3s%s\n",NR,(NR==L?">>>":""),$script }' L=$1 $script
}
trap 'err $LINENO' ERR

echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/scripts/01trim_ORFeome_Pool_A_S1_R1_001.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> outputs/log.txt

module add trimomatic

## This call to trimomatic removes illumina and epicentre adapters from ORFeome_Pool_A_S1_R1_001.fastq.gz:ORFeome_Pool_A_S1_R2_001.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
trimmomatic PE \
  -threads 1 \
  -phred33 \
  <(less r1.fastq.gz) <(less r2.fastq.gz) \
  ORFeome_Pool_A_S1_R1_001-trimmed_paired.fastq ORFeome_Pool_A_S1_R1_001-trimmed_unpaired.fastq \
  ORFeome_Pool_A_S1_R2_001-trimmed_paired.fastq ORFeome_Pool_A_S1_R2_001-trimmed_unpaired.fastq \
   ILLUMINACLIP:/cbcb/sw/RedHat-7-x86_64/common/local/perl/5.34/lib/site_perl/5.34.0/auto/share/dist/Bio-Adventure/genome/adapters.fa:2:20:10:2:keepBothReads \
  SLIDINGWINDOW:4:20 MINLEN:40 \
  1>outputs/01trimomatic/ORFeome_Pool_A_S1_R1_001-trimomatic.stdout \
  2>outputs/01trimomatic/ORFeome_Pool_A_S1_R1_001-trimomatic.stderr
excepted=$( { grep "Exception" "outputs/01trimomatic/ORFeome_Pool_A_S1_R1_001-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
  trimmomatic PE \
    -threads 1 \
    -phred33 \
    <(less r1.fastq.gz) <(less r2.fastq.gz) \
    ORFeome_Pool_A_S1_R1_001-trimmed_paired.fastq ORFeome_Pool_A_S1_R1_001-trimmed_unpaired.fastq \
    ORFeome_Pool_A_S1_R2_001-trimmed_paired.fastq ORFeome_Pool_A_S1_R2_001-trimmed_unpaired.fastq \
     SLIDINGWINDOW:4:25 MINLEN:50\
    1>outputs/01trimomatic/ORFeome_Pool_A_S1_R1_001-trimomatic.stdout \
    2>outputs/01trimomatic/ORFeome_Pool_A_S1_R1_001-trimomatic.stderr
fi
sleep 10
mv ORFeome_Pool_A_S1_R1_001-trimmed_paired.fastq ORFeome_Pool_A_S1_R1_001-trimmed.fastq
mv ORFeome_Pool_A_S1_R2_001-trimmed_paired.fastq ORFeome_Pool_A_S1_R2_001-trimmed.fastq

## Recompress the unpaired reads, this should not take long.
xz -9e -f ORFeome_Pool_A_S1_R1_001-trimmed_unpaired.fastq
xz -9e -f ORFeome_Pool_A_S1_R2_001-trimmed_unpaired.fastq
## Recompress the paired reads.
xz -9e -f ORFeome_Pool_A_S1_R1_001-trimmed.fastq
xz -9e -f ORFeome_Pool_A_S1_R2_001-trimmed.fastq
ln -sf ORFeome_Pool_A_S1_R1_001-trimmed.fastq.xz r1_trimmed.fastq.xz
ln -sf ORFeome_Pool_A_S1_R2_001-trimmed.fastq.xz r2_trimmed.fastq.xz


## The following lines give status codes and some logging
echo "## Job status: $? " >> outputs/log.txt
echo "## $(hostname) Finished ${SLURM_JOBID} 01trim_ORFeome_Pool_A_S1_R1_001.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

if [[ -x "$(command -v sstat)" && ! -z "${SLURM_JOBID}" ]]; then
  walltime=$(scontrol show job "${SLURM_JOBID}" | grep RunTime | perl -F'/\s+|=/' -lane '{print $F[2]}' 2>/dev/null)
  echo "#### walltime used by ${SLURM_JOBID} was: ${walltime:-null}" >> outputs/log.txt
  maxmem=$(sstat --format=MaxVMSize -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### maximum memory used by ${SLURM_JOBID} was: ${maxmem:-null}" >> outputs/log.txt
  avecpu=$(sstat --format=AveCPU -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### average cpu used by ${SLURM_JOBID} was: ${avecpu:-null}" >> outputs/log.txt
fi
## Adding a little logic to have skip finished jobs.
touch outputs/logs/01trim_ORFeome_Pool_A_S1_R1_001.finished

2.2 Mapping against pDONR

This was followed by my default mapping against the plasmid via hisat2.

I previously downloaded a copy of pDONR201 and ran the following on it:

cyoa --task index --method indexhisat --input pDONR201.fasta --species pDONR201

I did not generate a gff file for this, so htseq is not run.

The following invocation queues 2 jobs on the cluster, 1 for the mapping and 1 to sort/compress/index the alignment.

cd preprocessing
start=$(pwd)
for i in $(/bin/ls); do
    cd ${i}
    cyoa --task map --method hisat --species pDONR201 --input r1_trimmed.fastq.xz:r2_trimmed.fastq.xz
    cd $start
done

2.2.1 Individual pDONR mapping script

#!/usr/bin/env bash
#SBATCH --export=ALL
#SBATCH --mail-type=NONE
#SBATCH --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola
#SBATCH --partition=dpart
#SBATCH --qos=workstation --nice=10
#SBATCH --nodes=1 --requeue
#SBATCH --time=10:00:00
#SBATCH --job-name=40hisat2_pDONR201_genome
#SBATCH --mem=24G
#SBATCH --cpus-per-task=4
#SBATCH --output=outputs/log.txt.sbatch
set -o errexit
set -o errtrace
set -o pipefail
export LESS='--buffers 0'
script="$(pwd)/$0"
err() {
    echo "Error occurred:"
    awk 'NR>L-4 && NR<L+4 { printf "%-5d%3s%s\n",NR,(NR==L?">>>":""),$script }' L=$1 $script
}
trap 'err $LINENO' ERR

echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/scripts/40hisat2_pDONR201_genome.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> outputs/log.txt

module add hisat2 samtools htseq bamtools

## This is a hisat2 alignment of  -1 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/r1_trimmed.fastq.xz) -2 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/r2_trimmed.fastq.xz)  against /cbcbhomes/abelew/libraries/genome/indexes/pDONR201

mkdir -p outputs/40hisat2_pDONR201
hisat2 -x /cbcbhomes/abelew/libraries/genome/indexes/pDONR201  \
  -p 4 \
  -q   -1 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/r1_trimmed.fastq.xz) -2 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/r2_trimmed.fastq.xz)  \
  --phred33 \
  --un outputs/40hisat2_pDONR201/poola_unaldis_pDONR201_genome.fastq \
  --al outputs/40hisat2_pDONR201/poola_aldis_pDONR201_genome.fastq \
  --un-conc outputs/40hisat2_pDONR201/poola_unalcon_pDONR201_genome.fastq \
  --al-conc outputs/40hisat2_pDONR201/poola_alcon_pDONR201_genome.fastq \
  -S outputs/40hisat2_pDONR201/poola_pDONR201_genome.sam \
  2>outputs/40hisat2_pDONR201/hisat2_pDONR201_genome_poola.stderr \
  1>outputs/40hisat2_pDONR201/hisat2_pDONR201_genome_poola.stdout


## The following lines give status codes and some logging
echo "## Job status: $? " >> outputs/log.txt
echo "## $(hostname) Finished ${SLURM_JOBID} 40hisat2_pDONR201_genome.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

if [[ -x "$(command -v sstat)" && ! -z "${SLURM_JOBID}" ]]; then
  walltime=$(scontrol show job "${SLURM_JOBID}" | grep RunTime | perl -F'/\s+|=/' -lane '{print $F[2]}' 2>/dev/null)
  echo "#### walltime used by ${SLURM_JOBID} was: ${walltime:-null}" >> outputs/log.txt
  maxmem=$(sstat --format=MaxVMSize -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### maximum memory used by ${SLURM_JOBID} was: ${maxmem:-null}" >> outputs/log.txt
  avecpu=$(sstat --format=AveCPU -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### average cpu used by ${SLURM_JOBID} was: ${avecpu:-null}" >> outputs/log.txt
fi
## Adding a little logic to have skip finished jobs.
touch outputs/logs/40hisat2_pDONR201_genome.finished

The cyoa invocation which generated the above script also creates the following script along with a group of compression jobs which aggressively compress the fastq output files. I leave that as an exercise to the reader.

#!/usr/bin/env bash
#SBATCH --export=ALL
#SBATCH --mail-type=NONE
#SBATCH --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola
#SBATCH --partition=dpart
#SBATCH --qos=throughput --nice=10
#SBATCH --nodes=1 --requeue
#SBATCH --time=12:00:00
#SBATCH --job-name=40_2s2b_hisat2_pDONR201
#SBATCH --mem=12G
#SBATCH --cpus-per-task=1
#SBATCH --output=outputs/log.txt.sbatch
set -o errexit
set -o errtrace
set -o pipefail
export LESS='--buffers 0'
script="$(pwd)/$0"
err() {
    echo "Error occurred:"
    awk 'NR>L-4 && NR<L+4 { printf "%-5d%3s%s\n",NR,(NR==L?">>>":""),$script }' L=$1 $script
}
trap 'err $LINENO' ERR

echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/scripts/40_2s2b_hisat2_pDONR201.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> outputs/log.txt

module add samtools bamtools

## Converting the text sam to a compressed, sorted, indexed bamfile.
## Also printing alignment statistics to outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam.stats
## This job depended on: 240110

echo "Starting samtools"
if [[ ! -f "outputs/40hisat2_pDONR201/poola_pDONR201_genome.sam" ]]; then
    echo "Could not find the samtools input file."
    exit 1
fi
samtools view -u -t /cbcbhomes/abelew/libraries/genome/pDONR201.fasta \
  -S outputs/40hisat2_pDONR201/poola_pDONR201_genome.sam -o outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam  \
  2>outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam.err 1>outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam.out && \
echo 0
  samtools sort -l 9 outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam -o outputs/40hisat2_pDONR201/poola_pDONR201_genome-sorted.bam \
  2>outputs/40hisat2_pDONR201/poola_pDONR201_genome-sorted.stderr \
  1>outputs/40hisat2_pDONR201/poola_pDONR201_genome-sorted.stdout && \
  rm outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam && \
  rm outputs/40hisat2_pDONR201/poola_pDONR201_genome.sam && \
  mv outputs/40hisat2_pDONR201/poola_pDONR201_genome-sorted.bam outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam &&  samtools index outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam
echo 0
bamtools stats -in outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam 2>outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam.stats 1>&2
echo 0

## The following will fail if this is single-ended.
samtools view -b -f 2 -o outputs/40hisat2_pDONR201/poola_pDONR201_genome-paired.bam outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam && samtools index outputs/40hisat2_pDONR201/poola_pDONR201_genome-paired.bam
bamtools stats -in outputs/40hisat2_pDONR201/poola_pDONR201_genome-paired.bam 2>outputs/40hisat2_pDONR201/poola_pDONR201_genome-paired.stats 1>&2
bamtools filter -tag XM:0 -in outputs/40hisat2_pDONR201/poola_pDONR201_genome.bam -out outputs/40hisat2_pDONR201/poola_pDONR201_genome-sorted_nomismatch.bam &&
  samtools index outputs/40hisat2_pDONR201/poola_pDONR201_genome-sorted_nomismatch.bam

## The following lines give status codes and some logging
echo "## Job status: $? " >> outputs/log.txt
echo "## $(hostname) Finished ${SLURM_JOBID} 40_2s2b_hisat2_pDONR201.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

if [[ -x "$(command -v sstat)" && ! -z "${SLURM_JOBID}" ]]; then
  walltime=$(scontrol show job "${SLURM_JOBID}" | grep RunTime | perl -F'/\s+|=/' -lane '{print $F[2]}' 2>/dev/null)
  echo "#### walltime used by ${SLURM_JOBID} was: ${walltime:-null}" >> outputs/log.txt
  maxmem=$(sstat --format=MaxVMSize -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### maximum memory used by ${SLURM_JOBID} was: ${maxmem:-null}" >> outputs/log.txt
  avecpu=$(sstat --format=AveCPU -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### average cpu used by ${SLURM_JOBID} was: ${avecpu:-null}" >> outputs/log.txt
fi
## Adding a little logic to have skip finished jobs.
touch outputs/logs/40_2s2b_hisat2_pDONR201.finished

2.3 Mapping against PA01

Before mapping against the PA01 strain, I previously downloaded the reference assembly from NCBI as a complete genbank file and put it in my “${LIBDIR}/genome/” directory as ‘paeruginosa_pa01.gb’.

I then performed the following to convert it to fasta/gff and index it (note I did this quite a while ago, I am just writing it here for completeness), also here is an extra salmon index operation for fun:

start=$(pwd)
cd ~/libraries/genome
cyoa --task convert --method gb2gff --input paeruginosa_pa01.gb
cyoa --task index --method indexhisat --input paeruginosa_pa01.fasta --species paeruginosa_pa01
cyoa --task index --method indexsalmon --input paeruginosa_pa01_cds.fasta --species paeruginosa_pa01
cd $start

This is nearly identical to the previous step, except it uses the unaligned reads from the previous step, and I explicitly state that I want htseq to count the genes by the locus_tag tags.

cd preprocessing
start=$(pwd)
for i in $(/bin/ls); do
    cd ${i}
    cyoa --task map --method hisat --species paeruginosa_pa01 \
         --htseq_type gene --htseq_id locus_tag \
         --input $(/bin/ls outputs/40hisat2_pDONR201/*unalcon*.fastq.xz | tr '\n' ':' | sed 's/:$//g')
    cd $start
done

The resulting scripts are nearly identical to those for pDONR, except of course the species used and inputs:

#!/usr/bin/env bash
#SBATCH --export=ALL
#SBATCH --mail-type=NONE
#SBATCH --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola
#SBATCH --partition=dpart
#SBATCH --qos=workstation --nice=10
#SBATCH --nodes=1 --requeue
#SBATCH --time=10:00:00
#SBATCH --job-name=40hisat2_paeruginosa_pa01_genome
#SBATCH --mem=24G
#SBATCH --cpus-per-task=4
#SBATCH --output=outputs/log.txt.sbatch
set -o errexit
set -o errtrace
set -o pipefail
export LESS='--buffers 0'
script="$(pwd)/$0"
err() {
    echo "Error occurred:"
    awk 'NR>L-4 && NR<L+4 { printf "%-5d%3s%s\n",NR,(NR==L?">>>":""),$script }' L=$1 $script
}
trap 'err $LINENO' ERR

echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/scripts/40hisat2_paeruginosa_pa01_genome.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> outputs/log.txt

module add hisat2 samtools htseq bamtools

## This is a hisat2 alignment of  -1 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/outputs/40hisat2_pDONR201/poola_unalcon_pDONR201_genome.1.fastq.xz) -2 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/outputs/40hisat2_pDONR201/poola_unalcon_pDONR201_genome.2.fastq.xz)  against /cbcbhomes/abelew/libraries/genome/indexes/paeruginosa_pa01

mkdir -p outputs/40hisat2_paeruginosa_pa01
hisat2 -x /cbcbhomes/abelew/libraries/genome/indexes/paeruginosa_pa01  \
  -p 4 \
  -q   -1 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/outputs/40hisat2_pDONR201/poola_unalcon_pDONR201_genome.1.fastq.xz) -2 <(less /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/outputs/40hisat2_pDONR201/poola_unalcon_pDONR201_genome.2.fastq.xz)  \
  --phred33 \
  --un outputs/40hisat2_paeruginosa_pa01/poola_unaldis_paeruginosa_pa01_genome.fastq \
  --al outputs/40hisat2_paeruginosa_pa01/poola_aldis_paeruginosa_pa01_genome.fastq \
  --un-conc outputs/40hisat2_paeruginosa_pa01/poola_unalcon_paeruginosa_pa01_genome.fastq \
  --al-conc outputs/40hisat2_paeruginosa_pa01/poola_alcon_paeruginosa_pa01_genome.fastq \
  -S outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome.sam \
  2>outputs/40hisat2_paeruginosa_pa01/hisat2_paeruginosa_pa01_genome_poola.stderr \
  1>outputs/40hisat2_paeruginosa_pa01/hisat2_paeruginosa_pa01_genome_poola.stdout


## The following lines give status codes and some logging
echo "## Job status: $? " >> outputs/log.txt
echo "## $(hostname) Finished ${SLURM_JOBID} 40hisat2_paeruginosa_pa01_genome.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

if [[ -x "$(command -v sstat)" && ! -z "${SLURM_JOBID}" ]]; then
  walltime=$(scontrol show job "${SLURM_JOBID}" | grep RunTime | perl -F'/\s+|=/' -lane '{print $F[2]}' 2>/dev/null)
  echo "#### walltime used by ${SLURM_JOBID} was: ${walltime:-null}" >> outputs/log.txt
  maxmem=$(sstat --format=MaxVMSize -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### maximum memory used by ${SLURM_JOBID} was: ${maxmem:-null}" >> outputs/log.txt
  avecpu=$(sstat --format=AveCPU -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### average cpu used by ${SLURM_JOBID} was: ${avecpu:-null}" >> outputs/log.txt
fi
## Adding a little logic to have skip finished jobs.
touch outputs/logs/40hisat2_paeruginosa_pa01_genome.finished

I will leave the nearly identical conversion of the sam alignments to compressed/sorted/indexed bam as an exercise to the reader, but here is the resulting htseq script:

#!/usr/bin/env bash
#SBATCH --export=ALL
#SBATCH --mail-type=NONE
#SBATCH --open-mode=append
#SBATCH --chdir=/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola
#SBATCH --partition=dpart
#SBATCH --qos=throughput --nice=10
#SBATCH --nodes=1 --requeue
#SBATCH --time=12:00:00
#SBATCH --job-name=40_3hts_poola_all_hisat2_paeruginosa_pa01_sno_gene_Alias
#SBATCH --mem=6G
#SBATCH --cpus-per-task=1
#SBATCH --output=outputs/log.txt.sbatch
set -o errexit
set -o errtrace
set -o pipefail
export LESS='--buffers 0'
script="$(pwd)/$0"
err() {
    echo "Error occurred:"
    awk 'NR>L-4 && NR<L+4 { printf "%-5d%3s%s\n",NR,(NR==L?">>>":""),$script }' L=$1 $script
}
trap 'err $LINENO' ERR

echo "## Started /fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/scripts/40_3hts_poola_all_hisat2_paeruginosa_pa01_sno_gene_Alias.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> outputs/log.txt

module add htseq

## Counting the number of hits in outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired.bam for each feature found in /cbcbhomes/abelew/libraries/genome/paeruginosa_pa01.gff
htseq-count  --help 2>&1 | tail -n 3
htseq-count \
  -q -f bam -s no  --type gene  --idattr locus_tag \
  outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired.bam \
  /cbcbhomes/abelew/libraries_fs/genome/paeruginosa_pa01.gff \
  2>outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired_all_sno_gene_Alias.stderr \
  1>outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired_all_sno_gene_Alias.count && \
    xz -f -9e outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired_all_sno_gene_Alias.count 2>outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired_all_sno_gene_Alias.stderr.xz 1>outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome-paired_all_sno_gene_Alias.count.xz


## The following lines give status codes and some logging
echo "## Job status: $? " >> outputs/log.txt
echo "## $(hostname) Finished ${SLURM_JOBID} 40_3hts_poola_all_hisat2_paeruginosa_pa01_sno_gene_Alias.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

if [[ -x "$(command -v sstat)" && ! -z "${SLURM_JOBID}" ]]; then
  walltime=$(scontrol show job "${SLURM_JOBID}" | grep RunTime | perl -F'/\s+|=/' -lane '{print $F[2]}' 2>/dev/null)
  echo "#### walltime used by ${SLURM_JOBID} was: ${walltime:-null}" >> outputs/log.txt
  maxmem=$(sstat --format=MaxVMSize -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### maximum memory used by ${SLURM_JOBID} was: ${maxmem:-null}" >> outputs/log.txt
  avecpu=$(sstat --format=AveCPU -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "#### average cpu used by ${SLURM_JOBID} was: ${avecpu:-null}" >> outputs/log.txt
fi
## Adding a little logic to have skip finished jobs.
touch outputs/logs/40_3hts_poola_all_hisat2_paeruginosa_pa01_sno_gene_locus_tag.finished

3 Variant searching

I have two methods of variant searching in my repetoire, mpileup and freebayes. I only learned freebayes recently and so chose it. My current implementation also spins off a series of scripts to parse the freebayes output down into a simpler tsv format. Also, unlike the mapping script above, the freebayes invocation itself includes the conversion of the text vcf format to a sorted/compressed/index bcf format.

Note also that the previous step also creates 3 bam files:

1. \({sample}_\){species}_genome.bam: All alignments
2. \({sample}_\){species}_genome-paired.bam: All properly paired alignments, this is the default used by htseq when available.
3. \({sample}_\){species}_genome-sorted_nomismatch.bam: This is filtered for the set of alignments with no mismatches.

Since running these samples, I added another suffix to the hisat output directories, depending on whether the -k parameter is used to limit the number of alignments reported. That is not particularly relevant here, but worth noting.

cd preprocessing
start=$(pwd)
for i in $(/bin/ls); do
    cd ${i}
    cyoa --task snp --method freebayes --species paeruginosa_pa01 \
         --input $(/bin/ls outputs/40hisat2_paeruginosa_pa01/*_genome.bam)
    cd $start
done

3.0.1 Single freebayes invocation

#!/usr/bin/env bash
#SBATCH --export=ALL
#SBATCH --mail-type=NONE
#SBATCH --open-mode=append
#SBATCH --chdir=/mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola
#SBATCH --partition=dpart
#SBATCH --qos=workstation --nice=10
#SBATCH --nodes=1 --requeue
#SBATCH --time=10:00:00
#SBATCH --job-name=40freebayes_paeruginosa_pa01
#SBATCH --mem=48G
#SBATCH --cpus-per-task=4
#SBATCH --output=outputs/log.txt.sbatch
set -o errexit
set -o errtrace
set -o pipefail
export LESS='--buffers 0'
script="$(pwd)/$0"
err() {
    echo "Error occurred:"
    awk 'NR>L-4 && NR<L+4 { printf "%-5d%3s%s\n",NR,(NR==L?">>>":""),$script }' L=$1 $script
}
trap 'err $LINENO' ERR

echo "## Started /mnt/cbcb/fs01_abelew/cbcb-lab/nelsayed/scratch/atb/dnaseq/paeruginosa_plasmids_2022/preprocessing/poola/scripts/40freebayes_paeruginosa_pa01.sh at $(date) on $(hostname) with id ${SLURM_JOBID}." >> outputs/log.txt

module add freebayes libgsl/2.7.1 libhts/1.13 samtools/1.13 bcftools vcftools

## Use freebayes, bcftools, and vcfutils to get some idea about how many variant positions are in the data.

mkdir -p outputs/40freebayes_paeruginosa_pa01
freebayes -f /home/trey/libraries_fs/genome/paeruginosa_pa01.fasta \
  -v outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.vcf \
  outputs/40hisat2_paeruginosa_pa01/poola_paeruginosa_pa01_genome.bam \
  1>outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.stdout \
  2>outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.stderr
bcftools convert outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.vcf \
  -Ob -o outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.bcf \
  2>>outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.stderr \
  1>outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.stdout
bcftools index outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.bcf \
  2>>outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.stderr \
  1>>outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.stdout
rm outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01.vcf


## The following lines give status codes and some logging
echo "Job status: $? " >> outputs/log.txt
echo "  $(hostname) Finished ${SLURM_JOBID} 40freebayes_paeruginosa_pa01.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

if [[ -x "$(command -v sstat)" && ! -z "${SLURM_JOBID}" ]]; then
  walltime=$(scontrol show job "${SLURM_JOBID}" | grep RunTime | perl -F'/\s+|=/' -lane '{print $F[2]}' | head -n 1 2>/dev/null)
  echo "  walltime used by ${SLURM_JOBID} was: ${walltime:-null}" >> outputs/log.txt
  maxmem=$(sstat --format=MaxVMSize -n "${SLURM_JOBID}.batch" 2>/dev/null)
  echo "  maximum memory used by ${SLURM_JOBID} was: ${maxmem:-null}" >> outputs/log.txt
  echo "" >> outputs/log.txt
fi
## Adding a little logic to have skip finished jobs.
touch outputs/logs/40freebayes_paeruginosa_pa01.finished

Upon completion we have a series of files which in theory should provide everything of interest to Vince. Lets poke at them in R!

4 Vince’s Requests

Here is my email from Vince with his initial, pared down queries:

" Based on the data you have analyzed, can you please share with me a Excel spreadsheet (or CSV file) with 2 columns for each library and each PA gene:

1. Coverage (RPKM) or 0 for no reads
2. Number of SNPs "

4.1 Annotation

I always grab my favorite annotation sources, even if I will not really use them. I usually copy the gff and fasta for the version of the genome I used to the local directory, but I did not bother for this.

pa01_gff <- load_gff_annotations("~/libraries_fs/genome/paeruginosa_pa01.gff", id_col="locus_tag")

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = FALSE)

## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = FALSE)

## Returning a df with 15 columns and 5697 rows.

rownames(pa01_gff) <- pa01_gff[["locus_tag"]]

4.2 Expressionset

Let us now pull out the htseq count tables in order to make an rpkm table. This uses the metadata sample sheet I created in the sample_sheets/ directory.

I bolded the two most important columns for the purposes of these tasks.

pa01_expt <- create_expt("sample_sheets/all_samples.xlsx", gene_info=pa01_gff)

## Reading the sample metadata.

## The sample definitions comprises: 10 rows(samples) and 7 columns(metadata fields).

## Matched 5697 annotations and counts.

## Bringing together the count matrix and gene information.

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

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 5697 rows and 10 columns.

pa01_rpkm <- normalize_expt(pa01_expt, convert="rpkm")

## Plot a log-scale density of rpkm values.
plot_boxplot(pa01_expt)

## 12590 entries are 0.  We are on a log scale, adding 1 to the data.

plot_density(normalize_expt(pa01_rpkm, transform="log2", filter=TRUE))

## Removing 255 low-count genes (5442 remaining).

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

## $plot

## 
## $condition_summary
##    condition min   1st median  mean   3rd    max
## 1:         a   0 6.468  7.348 6.833 7.767  9.282
## 2:         b   0 6.174  7.332 6.796 7.811  9.268
## 3:         c   0 6.566  7.413 6.871 7.773  9.105
## 4:         d   0 6.535  7.335 6.756 7.766  9.166
## 5:       pmg   0 0.000  0.000 1.420 0.000 11.833
## 
## $batch_summary
##     batch min   1st median  mean   3rd   max
## 1: notopo   0 5.083  7.111 5.740 7.747 11.58
## 2:   topo   0 5.121  7.105 5.731 7.719 11.83
## 
## $sample_summary
##         sample min   1st median  mean   3rd    max
##  1:      poola   0 6.499  7.375 6.848 7.779  9.282
##  2:      poolb   0 6.119  7.318 6.783 7.827  9.268
##  3:      poolc   0 6.572  7.422 6.876 7.786  9.087
##  4:      poold   0 6.543  7.356 6.766 7.789  9.085
##  5:    poolpmg   0 0.000  0.000 1.425 0.000 11.584
##  6: poola_topo   0 6.441  7.319 6.818 7.756  9.165
##  7: poolb_topo   0 6.233  7.345 6.809 7.798  9.149
##  8: poolc_topo   0 6.560  7.405 6.866 7.762  9.105
##  9: poold_topo   0 6.529  7.308 6.746 7.739  9.166
## 10: poolpmg_tp   0 0.000  0.000 1.414 0.000 11.833
## 
## $table
##            id     sample counts condition  batch
##     1: PA0001      poola  7.797         a notopo
##     2: PA0002      poola  7.957         a notopo
##     3: PA0003      poola  8.387         a notopo
##     4: PA0004      poola  7.627         a notopo
##     5: PA0005      poola  7.700         a notopo
##    ---                                          
## 54416: PA5566 poolpmg_tp  0.000       pmg   topo
## 54417: PA5567 poolpmg_tp  0.000       pmg   topo
## 54418: PA5568 poolpmg_tp  0.000       pmg   topo
## 54419: PA5569 poolpmg_tp  0.000       pmg   topo
## 54420: PA5570 poolpmg_tp  0.000       pmg   topo

plot_nonzero(pa01_expt)$plot

rpkm_df <- as.data.frame(exprs(pa01_rpkm)) %>%
  round(2)
wanted_annotations <- fData(pa01_expt)[, c("gene", "locus_tag")]
wanted_data <- merge(wanted_annotations, rpkm_df,
                     by="row.names")
rownames(wanted_data) <- wanted_data[["Row.names"]]
wanted_data[["Row.names"]] <- NULL

5 Gathering mutations numbers by gene

I do not have an instant function for this task, but it should not prove difficult.

wanted_v2 <- wanted_data
for (sample in rownames(pData(pa01_expt))) {
  input_file <- glue::glue("preprocessing/{sample}/outputs/40freebayes_paeruginosa_pa01/paeruginosa_pa01_hits_by_gene.txt.xz")
  input_data <- readr::read_tsv(input_file, col_names=c("gene", "chr", "postion", "nt_string", "aa_string"))
  input_column <- input_data %>%
    dplyr::group_by(gene) %>%
    dplyr::summarize(n=n())
  colnames(input_column) <- c("gene", paste0(sample, "_n"))
  wanted_v2 <- merge(wanted_v2, input_column, by.x="row.names", by.y="gene")
  rownames(wanted_v2) <- wanted_v2[["Row.names"]]
  wanted_v2[["Row.names"]] <- NULL
}

## Rows: 2586 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2390 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2596 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2665 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 245 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2776 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2734 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2816 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 2921 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Rows: 280 Columns: 5

## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): gene, chr, nt_string, aa_string
## dbl (1): postion

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Write it out
written <- write_xlsx(data=wanted_v2, excel="excel/rpkm_mutations_per_gene.xlsx")

pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))
tmp <- sm(saveme(filename=this_save))