Preprocessing RNAseq samples from Leishmania panamensis.

index.html

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 tophat against the L.panamensis genome 5. Align with tophat against the L.braziliensis genome 6. Align with tophat against the H.sapiens genome

One small caveat, these sequences were created in two separate groups, hpgl0241-hpgl0322 and hpgl0630-hpgl0663.

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

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

4 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-hpgl0663.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 \
  hpgl0663_forward-trimmed.fastq.gz hpgl0663_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

## scripts/th_lpanamensis-hpgl0663.sh
mkdir -p outputs/tophat_lpanamensis && tophat  -g 1  \
  -G /cbcbhomes/abelew/libraries/genome/lpanamensis.gff \
  --b2-very-sensitive -p 4 -o outputs/tophat_lpanamensis \
/cbcbhomes/abelew/libraries/genome/indexes/lpanamensis \
  hpgl0663_forward-trimmed.fastq.gz hpgl0663_reverse-trimmed.fastq.gz && \
  samtools sort -l 9 -n outputs/tophat_lpanamensis/accepted_hits.bam outputs/tophat_lpanamensis/accepted_sorted && \
  mv outputs/tophat_lpanamensis/accepted_sorted.bam outputs/tophat_lpanamensis/accepted_hits.bam && \
  samtools index outputs/tophat_lpanamensis/accepted_hits.bam && \
  samtools sort -l 9 -n outputs/tophat_lpanamensis/unmapped.bam outputs/tophat_lpanamensis/unmapped_sorted && \
  mv outputs/tophat_lpanamensis/unmapped_sorted.bam outputs/tophat_lpanamensis/unmapped.bam && \
  samtools index outputs/tophat_lpanamensis/unmapped.bam

## scripts/th_hsapiens-hpgl0633.sh
mkdir -p outputs/tophat_hsapiens && tophat  -g 1  \
  -G /cbcbhomes/abelew/libraries/genome/hsapiens.gtf \
  --b2-very-sensitive -p 4 -o outputs/tophat_hsapiens \
/cbcbhomes/abelew/libraries/genome/indexes/hsapiens \
  hpgl0241_forward-trimmed.fastq hpgl0241_reverse-trimmed.fastq && \
  samtools sort -l 9 -n outputs/tophat_hsapiens/accepted_hits.bam outputs/tophat_hsapiens/accepted_sorted && \
  mv outputs/tophat_hsapiens/accepted_sorted.bam outputs/tophat_hsapiens/accepted_hits.bam && \
  samtools index outputs/tophat_hsapiens/accepted_hits.bam

5 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_lpanamensis-hpgl0241.sh
## Counting the number of hits in outputs/tophat_lpanamensis/accepted_hits.bam for each feature found in /cbcbhomes/abelew/libraries/genome/lpanamensis.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_lpanamensis/accepted_hits.bam /cbcbhomes/abelew/libraries/genome/lpanamensis.gff \
  1>outputs/tophat_lpanamensis/accepted_hits.count 2>outputs/tophat_lpanamensis/accepted_hits.error && \
    xz -9e outputs/tophat_lpanamensis/accepted_hits.count
## I am not going to belabor the point and print the htseq commands for all species.

6 Run Kallisto

A new alignment tool is kallisto, which skips all the intermediate steps and generates a set of pseudo-counts using a database of all possible transcripts as the index. The following describes one such invocation:

## scripts/kall_hsapiens-hpgl0241.sh
## This is a kallisto pseudoalignment of HPGL0241_R1_filtered_adaptertrim.fastq.xz HPGL0241_R2_filtered_adaptertrim.fastq against
## /cbcbhomes/abelew/libraries/genome/indexes/hsapiens.idx.
## This jobs depended on: 685883.cbcbsrv.umiacs.umd.edu
## I am not sure why I have that sleep...
mkdir -p outputs/kallisto_hsapiens && sleep 10 && \
kallisto quant --plaintext -t 4 -b 100 -o outputs/kallisto_hsapiens -i /cbcbhomes/abelew/libraries/genome/indexes/hsapiens.idx \
  HPGL0241_R1_filtered_adaptertrim.fastq.xz HPGL0241_R2_filtered_adaptertrim.fastq 2>outputs/kallisto/HPGL0241_R1_filtered_adaptertrim.fastq.xz:HPGL0241_R2_filtered_adaptertrim-kallisto.err 1>outputs/kallisto/HPGL0241_R1_filtered_adaptertrim.fastq.xz:HPGL0241_R2_filtered_adaptertrim-kallisto.out

index.html