L.braziliensis 2017: RNASeq sample preprocessing.

1 Copy reads to the working tree

I prefer to concatenate the multiple sequence files into a single file for each the forward(R1) and reverse(R2) sequences. Interestingly, for this data I decided to handle all steps in one shot with a single script which is invoked with 1 argument for each sample. The text of that script ‘copy_process.sh’ follows:

cd preprocessing
for i in $(/bin/ls hpgl*); do
  ./copy_process.sh ${i}
done

## startdir=$(pwd)
## for i in "$@"; do
##   cd $startdir
##   cyoa --task rnaseq --method copyraw --raw_dir "/cbcb/lab/nelsayed/raw_data/ade" --hpgl ${i}
##   cd $i
##   inputs="${i}_forward.fastq.gz:${i}_reverse.fastq.gz"
##   trimmed="${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz"
##   cyoa --task rnaseq --method fastqc --input ${inputs}
##   cyoa --task rnaseq --method trimomatic --input ${inputs}
##   trim_id=$(cat last_job.txt)
##   cyoa --task rnaseq --method tophat --input ${trimmed} -s lpanamensis --depends "${trim_id}"
##   lp_id=$(cat last_job.txt)
##   cyoa --task rnaseq --method tophat --input ${trimmed} -s lbraziliensis --depends "${trim_id}"
##   lb_id=$(cat last_job.txt)
##   cyoa --task rnaseq --method tophat --input ${trimmed} -s hsapiens --depends "${trim_id}"
##   hs_id=$(cat last_job.txt)
##   rm last_job.txt
##   cd $startdir
## done

This process is therefore done in 6 steps: 1. Copy the raw data from Najib’s raw_data directory 2. Perform fastqc on the raw reads 3. Perform trimomatic 4. Align with bowtie2 the L.braziliensis genome

2 Fastqc

The cyoa tool writes out a script which looks like the following:

## This FastQC run is against rnaseq data and is used for
## an initial estimation of the overall sequencing quality.
mkdir -p outputs/fastqc && \
  fastqc --extract -o outputs/fastqc hpgl0241_forward-trimmed.fastq hpgl0241_reverse-trimmed.fastq \
  2>outputs/fastqc.out 1>&2

A relevant plot from fastqc is (from sample hpgl0241):

Forward qualities: Reverse qualities:

3 Trimomatic

The trimomatic script looks like: It is worth noting that in some cases, the illumina adapter removal results in a java exception. If that happens, rather than lose the data I made my cyoa script re-invoke trimomatic without the illumina clipping. In addition, if things need to be redone, then the sequences may in fact have been moved into the sequences/ directory and compressed with xz -9e, thus the if [[]] at the beginning testing for those files.

## Trimomatic_Pairwise: In case a trimming needs to be redone...
if [[ ! -r "axa1_forward.fastq.gz" ]]; then
  if [[ -r "sequences/axa1_forward.fastq.xz" ]]; then
    mv sequences/axa1_forward.fastq.xz . && pxz -d axa1_forward.fastq.xz && pigz axa1_forward.fastq && mv sequences/axa1_reverse.fastq.xz . &&
 pxz -d axa1_reverse.fastq.xz && pigz axa1_reverse.fastq
  else
    echo "Missing files. Did not find axa1_forward.fastq.gz nor sequences/axa1_forward.fastq.xz"
    exit 1
  fi
fi
trimomatic PE -threads 1 -phred33 axa1_forward.fastq.gz axa1_reverse.fastq.gz axa1_forward-trimmed_paired.fastq.gz axa1_forward-trimmed_unpair
ed.fastq.gz axa1_reverse-trimmed_paired.fastq.gz axa1_reverse-trimmed_unpaired.fastq.gz ILLUMINACLIP:/cbcbhomes/abelew/libraries/adapters.fa:2:20:4 SL
IDINGWINDOW:4:25 1>outputs/axa1-trimomatic.out 2>&1
excepted=$(grep "Exception" outputs/axa1-trimomatic.out)
## 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 axa1_forward.fastq.gz axa1_reverse.fastq.gz axa1_forward-trimmed_paired.fastq.gz axa1_forward-trimmed_unpa
ired.fastq.gz axa1_reverse-trimmed_paired.fastq.gz axa1_reverse-trimmed_unpaired.fastq.gz SLIDINGWINDOW:4:25 1>>outputs/axa1-trimomatic.out 2>&1
fi
sleep 10
mv axa1_forward-trimmed_paired.fastq.gz axa1_forward-trimmed.fastq.gz && mv axa1_reverse-trimmed_paired.fastq.gz axa1_reverse-trimmed.fastq.gz

4 Bowtie 2

The end result should be a set of files with trimmed, paired sequences separated from unpaired sequences.

5 Tophat

The tophat process is actually a couple of steps rolled into one, the trimmed/paired reads are fed to tophat, then the accepted_hits.bam is tested for properly paired reads and those reads are separated into accepted_paired.bam. These 2 bam files are sorted and indexed, and finally passed to htseq. Thus the tophat scripts look like:

These scripts however, do not include the additional processing step which split the accepted_hits.bam into accepted_paired.bam I am not sure why, but these scripts are missing the following step which was performed but not properly logged:

if [ -r "${tophat_dir}/accepted_hits.bam" ]; then
  samtools view -b -f 2 ${tophat_dir}/accepted_hits.bam > ${tophat_dir}/accepted_paired.bam && samtools index ${tophat_dir}/accepted_paired.bam
fi

I know this step was performed, because the files accepted_paired.bam exists for each directory, but the option -f 2 says that the ‘2’ flag must be set, which is the flag for an alignment which is paired.

## scripts/th_lbraziliensis-axa1.sh
mkdir -p outputs/tophat_lbraziliensis && tophat  -g 1  \
  -G /cbcbhomes/abelew/libraries/genome/lbraziliensis.gff \
  --b2-very-sensitive -p 4 -o outputs/tophat_lbraziliensis \
/cbcbhomes/abelew/libraries/genome/indexes/lbraziliensis \
  axa1_forward-trimmed.fastq.gz axa1_reverse-trimmed.fastq.gz && \
  samtools sort -l 9 -n outputs/tophat_lbraziliensis/accepted_hits.bam outputs/tophat_lbraziliensis/accepted_sorted && \
  mv outputs/tophat_lbraziliensis/accepted_sorted.bam outputs/tophat_lbraziliensis/accepted_hits.bam && \
  samtools index outputs/tophat_lbraziliensis/accepted_hits.bam && \
  samtools sort -l 9 -n outputs/tophat_lbraziliensis/unmapped.bam outputs/tophat_lbraziliensis/unmapped_sorted && \
  mv outputs/tophat_lbraziliensis/unmapped_sorted.bam outputs/tophat_lbraziliensis/unmapped.bam && \
  samtools index outputs/tophat_lbraziliensis/unmapped.bam

6 Run htseq

The above scripts create sorted accepted_hits.bam files, I think that I manually made the accepted_paired.bam files and added those steps to the cyoa script post-facto, but the end result is the same and so I needed to run htseq on both:

## scripts/th_lbraziliensis-hpgl0241.sh
## Counting the number of hits in outputs/tophat_lbraziliensis/accepted_hits.bam for each feature found in /cbcbhomes/abelew/libraries/genome/lbraziliensis.gff
## Is this stranded? no.  The defaults of htseq are:
##  --order=name --idattr=gene_id --minaqual=10 --type=exon --stranded=yes --mode=union
htseq-count -q -f bam -s no  -i ID  \
  outputs/tophat_lbraziliensis/accepted_hits.bam /cbcbhomes/abelew/libraries/genome/lbraziliensis.gff \
  1>outputs/tophat_lbraziliensis/accepted_hits.count 2>outputs/tophat_lbraziliensis/accepted_hits.error && \
    xz -9e outputs/tophat_lbraziliensis/accepted_hits.count
## I am not going to belabor the point and print the htseq commands for all species.

if (!isTRUE(get0("skip_load"))) {
  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))
}

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

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

## > git reset 2a0661d6e37f8a3d8831eb3bbd6347c0d9c4b3b7

## R> packrat::restore()

## This is hpgltools commit: Thu Mar 29 16:59:07 2018 -0400: 2a0661d6e37f8a3d8831eb3bbd6347c0d9c4b3b7

## Saving to index-v20170201.rda.xz