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.

3 A suggestion from Nour

Nour made a good suggestion: gatk/picard provides a tool to mark/deduplicate sequence reads given a set of mapped reads. I applied that to this data set and noted some interesting differences. This is implemented as a method called ‘MarkDuplicates’ or some such, I will include a script stanza below showing its usage.

3.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.

3.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

3.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

3.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

3.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

4 Deduplication

Nour’s suggestion resulted in the following addition (which has since been added as a precursor to my freebayes/mpileup invocation):

mkdir -p ${freebayes_dir}
gatk MarkDuplicates \\
  -I $options->{input} \\
  -O ${deduplicated} \\
  -M ${marked} --REMOVE_DUPLICATES true --COMPRESSION_LEVEL 9 \\
  2>${gatk_stdout} \\
  1>${gatk_stderr}
samtools index ${deduplicated}

5 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

5.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!

6 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 "

6.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"]]

6.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 9 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")
pa01_cpm <- normalize_expt(pa01_expt, convert="cpm")

## 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

7 Gathering mutations numbers by gene

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

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:hpgltools':
## 
##     combine

## The following object is masked from 'package:testthat':
## 
##     matches

## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, setdiff, union

## The following object is masked from 'package:GenomeInfoDb':
## 
##     intersect

## The following objects are masked from 'package:IRanges':
## 
##     collapse, desc, intersect, setdiff, slice, union

## The following objects are masked from 'package:S4Vectors':
## 
##     first, intersect, rename, setdiff, setequal, union

## The following object is masked from 'package:matrixStats':
## 
##     count

## The following object is masked from 'package:Biobase':
## 
##     combine

## The following objects are masked from 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

wanted_v2 <- wanted_data

for (sample in rownames(pData(pa01_expt))) {
  full_file <- glue::glue("preprocessing/{sample}/outputs/50freebayes_paeruginosa_pa01/all_tags.txt.xz")
  input_file <- glue::glue("preprocessing/{sample}/outputs/50freebayes_paeruginosa_pa01/variants_by_gene.txt.xz")
  if (!file.exists(full_file)) {
    next
  }
  input_data <- readr::read_tsv(
                           input_file, show_col_types=FALSE, skip=1,
                           col_names=c("gene", "chr", "position", "nt_string", "aa_string")) %>%
    as.data.frame()
  input_data[["from"]] <- gsub(x=input_data[["nt_string"]], pattern="(\\w+)_(\\w+)",
                               replacement="\\1", perl=TRUE)
  input_data[["to"]] <- gsub(x=input_data[["nt_string"]], pattern="(\\w+)_(\\w+)",
                             replacement="\\2", perl=TRUE)
  input_data[["encoded"]] <- paste0("chr_", input_data[["chr"]], "_pos_",
                                    input_data[["position"]], "_ref_",
                                    input_data[["from"]], "_alt_",
                                    input_data[["to"]])
  penetrance <- readr::read_tsv(full_file, show_col_types=FALSE) %>%
    select(c("position", "PAIRED"))
  input_data <- merge(input_data, penetrance, by.x="encoded", by.y="position", all.x=TRUE)
  input_column <- input_data %>%
    dplyr::group_by(gene) %>%
    dplyr::summarize(n=n())
  colnames(input_column) <- c("gene", paste0(sample, "_n"))
  penetrance_column <- input_data %>%
    group_by(gene) %>%
    dplyr::mutate(aa_penetrance=paste0(aa_string, '_', PAIRED, collapse=",")) %>%
    select(gene, aa_penetrance) %>%
    distinct()
  colnames(penetrance_column) <- c("gene", paste0(sample, "_aa_penetrance"))
  input_columns <- merge(input_column, penetrance_column, by="gene")
  wanted_v2 <- merge(wanted_v2, input_columns, by.x="row.names", by.y="gene", all.x=TRUE)
  rownames(wanted_v2) <- wanted_v2[["Row.names"]]
  wanted_v2[["Row.names"]] <- NULL
}

## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...

## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...

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

8 Try visualizing some of this information

pa_snp_expt <- count_expt_snps(pa01_expt, annot_column="freebayestable", snp_column="PAIRED")

## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...
## New names:
## * DP -> DP...3
## * RO -> RO...8
## * AO -> AO...9
## * QR -> QR...12
## * QA -> QA...13
## * ...

plot_libsize(pa_snp_expt)$plot

test <- plot_corheat(pa_snp_expt)

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))