L.major 2018: Mouse PMN preprocessing.

1 Preprocessing some peculiar samples.

2 Separating samples by hpgl identifier

Najib provided a series of his traditional hpgl identifiers for this data set. I therefore manually copied out the samples by sequence file names in the run0395 and run0396 directories into the appropriate directories. This is therefore an obvious potential source of error. I therefore did this with a second person watching (thanks, Yoann!) and checked each sample 2 times before moving forward.

2.1 New samples

20180402: The folks at the Sacks lab kindly provided a new set of samples, for which Najib provided a fresh set of IDs (HPGL1065-1088). I did the same tasks using the sample sheet he provided me. The source of human error remains, I did double-check my work, but human error in copying remains a possibility.

For each hpgl identifier, I did the following:

1. Check the sample sheet and copy the files in lane1/lane2 corresponding to the 132* identifier to the corresponding HPGL identifier directory.
2. cd into each HPGL directory / unprocessed and:

1. cp 132*_L001_R1* hpglxxxx_forward.fastq.gz
2. cat 132R1 >> hpglxxxx_forward.fastq.gz (Note I did not do 132R1 as that would be redundant, but instead used tab completion in bash to pick up the correct filenames. I did human error the second directory but happily noticed it immediately and do not think I messed up any thereafter.
3. cp 132L001_R2 hpglxxxx_reverse.fastq.gz
4. cat 132R2 >> hpglxxxx_reverse.fastq.gz

Once I finished this, invoked cyoa with options appropriate for trimming/fastqc:

cyoa --task rnaseq --method fastqc --input hpglxxxx_forward.fastq.gz:hpglxxxx_reverse.fastq.gz
cyoa --task rnaseq --method trim --input hpglxxxx_forward.fastq.gz:hpglxxxx_reverse.fastq.gz

3 Preprocessing by sample

The preprocessing of this data happened in two parts and was performed via the cyoa script, which attempts to simplify the creation of jobs on our computer cluster.

I split the preprocessing tasks into two parts, the first part consists of the initial sequencing quality assurance and adapter trimming.

3.1 Part 1

The following is contained in a shell script in preprocessing/ called process_part1.sh

#!/usr/bin/env bash
start=$(pwd)
for i in $(/bin/ls -d hpgl*); do
  cd ${i}
  echo "Starting ${i}"
  cyoa --task rnaseq --method trimom --input ${i}_forward.fastq.gz:${i}_reverse.fastq.gz
  cyoa --task rnaseq --method fastq --input ${i}_forward.fastq.gz:${i}_reverse.fastq.gz
  echo "Finishing ${i}"
  cd ${start}
done

3.2 Part 2

Once the sequence adapters were removed, I perform the set of tasks which depend on the trimmed data. It is worth noting that I tried two (pseudo)aligners for the mouse genome, tophat and kallisto; and two aligners for the parasite genomes, tophat and bowtie2.

In addition, I attempted alignments against L.major and L.amazonensis. My reasoning will become clear soon I think.

#!/usr/bin/env bash
start=$(pwd)
for i in $(/bin/ls -d hpgl084[0-2]); do
  cd "${i}" || exit
  echo "Starting ${i}"
  cyoa --task rnaseq --method biopiece --input "${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz"
  cyoa --task rnaseq --method kallisto --input "${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz" --species mmusculus
  cyoa --task rnaseq --method tophat --input "${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz" --species mmusculus --htseq_type exon --htseq_stranded yes
  cyoa --task rnaseq --method bt2 --input "${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz" --species lamazonensis --htseq_type exon --htseq_stranded yes
  cyoa --task rnaseq --method bt2 --input "${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz" --species lmajor --htseq_type exon --htseq_stranded yes
  cyoa --task rnaseq --method tophat --input "${i}_forward-trimmed.fastq.gz:${i}_reverse-trimmed.fastq.gz" --species lmajor --htseq_type exon --htseq_stranded yes
  echo "Finishing ${i}"
  cd "${start}" || exit
done

4 What happened in each step?

The cyoa script in turn creates and executes jobs on our torque cluster, I am copy/pasting them below:

Two caveats:

1. There is a common header and footer for all of these scripts. I will include it for this script, but skip it for the following:
2. I am going to print here a copy of my script for each step for the first sample and leave as an exercise for the reader the rest.

4.1 Fastqc

The first step in most(all) of our data sets is to make sure that the sequencer provided high-quality data at the expected length(s).

There is only one interesting point in this script, the –extract argument was provided so that a dependent job may be run in order to extract some statistics from the fastqc results.

The following is preprocessing/hpgl0480/scripts/00fqc-hpgl0840.sh

#!/usr/bin/env bash
#PBS -V -S /usr/bin/bash -q workstation
#PBS -d /cbcb/nelsayed-scratch/atb/rnaseq/lmajor_sax/preprocessing/hpgl0840
#PBS -N fqc-hpgl0840 -l mem=6gb -l walltime=10:00:00 -l ncpus=8
#PBS -o /cbcb/nelsayed-scratch/atb/rnaseq/lmajor_sax/preprocessing/hpgl0840/outputs/qsub/fqc-hpgl0840.qsubout -j oe -V -m n
echo "####Started /cbcb/nelsayed-scratch/atb/rnaseq/lmajor_sax/preprocessing/hpgl0840/scripts/00fqc-hpgl0840.sh at $(date)" >> outputs/log.txt
cd /cbcb/nelsayed-scratch/atb/rnaseq/lmajor_sax/preprocessing/hpgl0840 || exit

## 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 hpgl0840_forward.fastq.gz hpgl0840_reverse.fastq.gz \
  2>outputs/fastqc.out 1>&2

## The following lines give status codes and some logging
echo $? > outputs/status/fqc-hpgl0840.status
echo "###Finished ${PBS_JOBID} 00fqc-hpgl0840.sh at $(date), it took $(( SECONDS / 60 )) minutes." >> outputs/log.txt

walltime=$(qstat -f -t "${PBS_JOBID}" | grep 'resources_used.walltime' | awk -F ' = ' '{print $2}')
echo "####PBS walltime used by ${PBS_JOBID} was: ${walltime:-null}" >> outputs/log.txt
mem=$(qstat -f -t | grep "${PBS_JOBID}" | grep 'resources_used.mem' | awk -F ' = ' '{print $2}')
echo "####PBS memory used by ${PBS_JOBID} was: ${mem:-null}" >> outputs/log.txt
vmmemory=$(qstat -f -t "${PBS_JOBID}" | grep 'resources_used.vmem' | awk -F ' = ' '{print $2}')
echo "####PBS vmemory used by ${PBS_JOBID} was: ${vmmemory:-null}" >> outputs/log.txt
cputime=$(qstat -f -t "${PBS_JOBID}" | grep 'resources_used.cput' | awk -F ' = ' '{print $2}')
echo "####PBS cputime used by ${PBS_JOBID} was: ${cputime:-null}" >> outputs/log.txt
##qstat -f -t ${PBS_JOBID} >> outputs/log.txt

4.1.1 Fastqc statistics

The fastqc statistics are generated separately for the forward/reverse reads.

The following is preprocessing/hpgl0480/scripts/01fqcst_forward-hpgl0840.sh

## This is a stupidly simple job to collect alignment statistics.

if [ ! -r outputs/fastqc_forward_stats.csv ]; then
  echo "name,total_reads,poor_quality,per_quality,per_base_content,per_sequence_gc,per_base_n,per_seq_length,over_rep,adapter_content,kmer_content" > outputs/fastqc_forward_stats.csv
fi
total_reads_tmp=$(grep "^Total Sequences" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
total_reads=${total_reads_tmp:-0}
poor_quality_tmp=$(grep "^Sequences flagged as poor quality" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
poor_quality=${poor_quality_tmp:-0}
per_quality_tmp=$(grep "Per base sequence quality" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
per_quality=${per_quality_tmp:-0}
per_base_content_tmp=$(grep "Per base sequence content" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
per_base_content=${per_base_content_tmp:-0}
per_sequence_gc_tmp=$(grep "Per sequence GC content" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
per_sequence_gc=${per_sequence_gc_tmp:-0}
per_base_n_tmp=$(grep "Per base N content" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
per_base_n=${per_base_n_tmp:-0}
per_seq_length_tmp=$(grep "Sequence Length Distribution" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
per_seq_length=${per_seq_length_tmp:-0}
over_rep_tmp=$(grep "Overrepresented sequences" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
over_rep=${over_rep_tmp:-0}
adapter_content_tmp=$(grep "Adapter Content" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
adapter_content=${adapter_content_tmp:-0}
kmer_content_tmp=$(grep "Kmer Content" outputs/fastqc/hpgl0840_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
kmer_content=${kmer_content_tmp:-0}

stat_string=$(printf "hpgl0840,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s" "${total_reads}" "${poor_quality}" "${per_quality}" "${per_base_content}" "${per_sequence_gc}" "${per_base_n}" "${per_seq_length}" "${over_rep}" "${adapter_content}" "${kmer_content}")
echo "$stat_string" >> outputs/fastqc_forward_stats.csv

4.2 Trimomatic

Trimomatic is perhaps the most important step in the process, as it removes the sometimes very abundant sequencing adapters, which if forgotten, would destroy all the following steps.

The following is preprocessing/hpgl0840/scripts/05trim-hpgl0840.sh

This job has a little bit of logic in it in case trimomatic fails.

## This call to trimomatic removes illumina and epicentre adapters from hpgl0840_forward.fastq.gz:hpgl0840_reverse.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.

## Trimomatic_Pairwise: In case a trimming needs to be redone...
if [[ ! -r "hpgl0840_forward.fastq.gz" ]]; then
  if [[ -r "sequences/hpgl0840_forward.fastq.xz" ]]; then
    mv sequences/hpgl0840_forward.fastq.xz . && pxz -d hpgl0840_forward.fastq.xz && pigz hpgl0840_forward.fastq &&\
       mv sequences/hpgl0840_reverse.fastq.xz . && pxz -d hpgl0840_reverse.fastq.xz && pigz hpgl0840_reverse.fastq
  else
    echo "Missing files. Did not find hpgl0840_forward.fastq.gz nor sequences/hpgl0840_forward.fastq.xz"
    exit 1
  fi
fi
## Note that trimomatic prints all output and errors to STDERR, so send both to output
trimomatic PE -threads 1 -phred33 hpgl0840_forward.fastq.gz hpgl0840_reverse.fastq.gz hpgl0840_forward-trimmed_paired.fastq.gz hpgl0840_forward-trimmed_unpaired.fastq.gz hpgl0840_reverse-trimmed_paired.fastq.gz hpgl0840_reverse-trimmed_unpaired.fastq.gz \
    ILLUMINACLIP:/cbcbhomes/abelew/libraries/adapters.fa:2:20:4 SLIDINGWINDOW:4:25 \
    1>outputs/hpgl0840-trimomatic.out 2>&1
excepted=$(grep "Exception" outputs/hpgl0840-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 hpgl0840_forward.fastq.gz hpgl0840_reverse.fastq.gz hpgl0840_forward-trimmed_paired.fastq.gz hpgl0840_forward-trimmed_unpaired.fastq.gz hpgl0840_reverse-trimmed_paired.fastq.gz hpgl0840_reverse-trimmed_unpaired.fastq.gz SLIDINGWINDOW:4:25 \
    1>outputs/hpgl0840-trimomatic.out 2>&1
fi
sleep 10
mv hpgl0840_forward-trimmed_paired.fastq.gz hpgl0840_forward-trimmed.fastq.gz && mv hpgl0840_reverse-trimmed_paired.fastq.gz hpgl0840_reverse-trimmed.fastq.gz

4.2.1 Trimomatic statistics

In the same way, the text output of trimomatic is parsed to generate some statistics.

The following is preprocessing/hpgl0840/scripts/06trimst-hpgl0840.sh

## This is a stupidly simple job to collect trimomatic statistics

if [ ! -r outputs/trimomatic_stats.csv ]; then
  echo "total_reads,surviving_both,surviving_forward,surviving_reverse,dropped_reads" > outputs/trimomatic_stats.csv
fi
total_reads_tmp=$(grep "^Input Read Pairs" outputs/hpgl0840-trimomatic.out | awk '{print $4}')
total_reads=${total_reads_tmp:-0}
surviving_both_tmp=$(grep "^Input Read Pairs" outputs/hpgl0840-trimomatic.out | awk '{print $7}')
surviving_both=${surviving_both_tmp:-0}
surviving_forward_tmp=$(grep "^Input Read Pairs" outputs/hpgl0840-trimomatic.out | awk '{print $12}')
surviving_forward=${surviving_forward_tmp:-0}
surviving_reverse_tmp=$(grep "^Input Read Pairs" outputs/hpgl0840-trimomatic.out | awk '{print $17}')
surviving_reverse=${surviving_reverse_tmp:-0}
dropped_reads_tmp=$(grep "^Input Read Pairs" outputs/hpgl0840-trimomatic.out | awk '{print $20}')
dropped_reads=${dropped_reads_tmp:-0}

stat_string=$(printf "hpgl0840,%s,%s,%s,%s,%s" "${total_reads}" "${surviving_both}" "${surviving_forward}" "${surviving_reverse}" "${dropped_reads}")
echo "$stat_string" >> outputs/trimomatic_stats.csv

4.3 Biopieces

Biopieces is an interesting project to provide a suite of simple tools which may be connected to perform interesting tasks. I use it to plot some graphs of the sequencing quality/attributes.

The following is preprocessing/hpgl0840/scripts/02biop-hpgl0840trimmed.sh

## This script uses biopieces to draw some simple graphs of the sequence.

## Do not forget that _only_ the last command in a biopieces string is allowed to have the -x.
mkdir -p outputs/biopieces
less hpgl0840_reverse-trimmed.fastq | read_fastq -i - -e base_33 |\
 plot_scores -T 'Quality Scores' -t svg -o outputs/biopieces/hpgl0840_reverse-trimmed_quality_scores.svg |\
 plot_nucleotide_distribution -T 'NT. Distribution' -t svg -o outputs/biopieces/hpgl0840_reverse-trimmed_ntdist.svg |\
 plot_lendist -T 'Length Distribution' -k SEQ_LEN -t svg -o outputs/biopieces/hpgl0840_reverse-trimmed_lendist.svg |\
 analyze_gc |\
     bin_vals -b 20 -k 'GC%' |\
     plot_distribution -k 'GC%_BIN' -t svg -o outputs/biopieces/hpgl0840_reverse-trimmed_gcdist.svg |\
 analyze_gc |\
     mean_vals -k 'GC%' -o outputs/biopieces/hpgl0840_reverse-trimmed_gc.txt |\
 count_records -o outputs/biopieces/hpgl0840_reverse-trimmed_count.txt -x

4.4 Kallisto Mouse

Kallisto is a super-fast aligner for estimating relative transcript abundance. I increasingly think it is better than tophat.

I added a little logic to simplify its output to provide a count table-like file.

The following is preprocessing/hpgl0840/scripts/30kall_mmusculus-hpgl0840trimmed.sh

## This is a kallisto pseudoalignment of hpgl0840_forward-trimmed.fastq.gz hpgl0840_reverse-trimmed.fastq.gz against
## /cbcbhomes/abelew/libraries/genome/indexes/mmusculus.idx.
## This jobs depended on:
## Other candidates for making a pretty count table include:
##  perl -F'\t' -a -n -e 'print "$F[0] $F[3]\n"' outputs/kallisto_mmusculus/abundance.tsv > outputs/kallisto_mmusculus/abundance.count
##   awk '{printf("%s %s\n", $1, $4)}' outputs/kallisto_mmusculus/abundance.tsv > outputs/kallisto_mmusculus/abundance.count
## because kallisto is so fast

mkdir -p outputs/kallisto_mmusculus && sleep 10 && \
  kallisto quant  --plaintext -b 100 -o outputs/kallisto_mmusculus -i /cbcbhomes/abelew/libraries/genome/indexes/mmusculus.idx \
    hpgl0840_forward-trimmed.fastq.gz hpgl0840_reverse-trimmed.fastq.gz 2>outputs/kallisto_mmusculus/kallisto_mmusculus.stderr 1>outputs/kallisto_mmusculus/kallisto_mmusculus.txt && \
  cut -d "  " -f 1,4 outputs/kallisto_mmusculus/abundance.tsv > outputs/kallisto_mmusculus/abundance.count && \
  gzip outputs/kallisto_mmusculus/abundance.count && xz -9e *.tsv

I currently do not have a job to create kallisto statistics, which is dumb, fix that.

4.5 Tophat Mouse

This step is where I started to see that something is strange in the data. But before I get into a discussion of what is odd, here is a tophat invocation:

Note that the example is actually from hpgl0850 because 840 is promastigote and unlikely to be interesting for the mouse genome (though I did run it for that sample).

There is a little bit of interesting logic in this:

1. I make sure to sort the accepted_hits.bam because it is annoying to not be able to look at these in IGV.
2. I take an extra step to extract the perfectly paired aligned reads (The samtools view -b -f 2 line). In the best case scenario, this removes nothing.

The following is preprocessing/hpgl0850/scripts/31th_mmusculus-hpgl0850trimmed.sh

## I still have no clue what I am doing when I use tophat...
## However, I know that -g 1 will allow only 1 hit in the case of multihits, but randomly place it
## From the manual:  "If there are more alignments with the same score than this
## number, TopHat will randomly report only this many alignments"
## -N 1 will discard anything with >1 mismatch (default is 2)
## -r adjusts the allowable mean distance between the paired reads
## --mate-std-dev sets the deviation of -r
## --microexon-search will tell it to search short exons for reads >=50

mkdir -p outputs/tophat_mmusculus && tophat  -g 1 --b2-very-sensitive  \
  -G /cbcbhomes/abelew/libraries/genome/mmusculus.gtf \
  -p 4 -o outputs/tophat_mmusculus \
  /cbcbhomes/abelew/libraries/genome/indexes/mmusculus \
  hpgl0850_forward-trimmed.fastq.gz hpgl0850_reverse-trimmed.fastq.gz 2>outputs/tophat.out 1>&2 && \
 samtools sort -l 9 -n outputs/tophat_mmusculus/accepted_hits.bam outputs/tophat_mmusculus/accepted_sorted && \
 mv outputs/tophat_mmusculus/accepted_sorted.bam outputs/tophat_mmusculus/accepted_hits.bam && \
 samtools index outputs/tophat_mmusculus/accepted_hits.bam && \
 samtools sort -l 9 -n outputs/tophat_mmusculus/unmapped.bam outputs/tophat_mmusculus/unmapped_sorted && \
 mv outputs/tophat_mmusculus/unmapped_sorted.bam outputs/tophat_mmusculus/unmapped.bam && \
 samtools index outputs/tophat_mmusculus/unmapped.bam

if [ -r "outputs/tophat_mmusculus/accepted_hits.bam" ]; then
  samtools view -b -f 2 outputs/tophat_mmusculus/accepted_hits.bam > outputs/tophat_mmusculus/accepted_paired.bam && \
    samtools index outputs/tophat_mmusculus/accepted_paired.bam
fi

For my initial run, I left all the defaults in place and aligned against all reads in the sample. My first moment of worry therefore happened when I looked at the alignment statistics, for hpgl0850 they are:

Keep in mind that hpgl0850 is a 12 hour infected PMN sample, I therefore assumed a significant majority of the full set of reads would map to the mouse genome.

Left reads: Input : 5557510 Mapped : 253591 ( 4.6% of input) of these: 15051 ( 5.9%) have multiple alignments (36717 have >1) Right reads: Input : 5557510 Mapped : 386499 ( 7.0% of input) of these: 15051 ( 3.9%) have multiple alignments (56680 have >1) 5.8% overall read mapping rate.

Aligned pairs: 183046 of these: 15051 ( 8.2%) have multiple alignments 30642 (16.7%) are discordant alignments 2.7% concordant pair alignment rate.

Clearly this is not true, I think I will be able to answer the why in a little bit…

4.5.1 Tophat statistics

I generate a similar csv file for the tophat alignments:

The following is preprocessing/hpgl0850/scripts/33tpstats_mmusculus-hpgl0850trimmed.sh

## This is a stupidly simple job to collect tophat alignment statistics.

if [ ! -r outputs/tophat_stats.csv ]; then
  echo "basename,species,original_reads,aligned_reads,failed_reads,rpm,count_table" > outputs/tophat_stats.csv
fi
bamtools stats < outputs/tophat_mmusculus/accepted_hits.bam 2>outputs/tophat_mmusculus/accepted_hits.bam.stats 1>&2 && bamtools stats < outputs/tophat_mmusculus/unmapped.bam 2>outputs/tophat_mmusculus/unmapped.bam.stats 1>&2
original_reads_tmp=$(grep "^Input Reads" outputs/hpgl0850-trimmed-trimomatic.out | awk '{print $3}' | sed 's/ //g')
original_reads=${original_reads_tmp:-0}
reads_tmp=$(grep "^reads_in " outputs/tophat_mmusculus/prep_reads.info | awk -F= '{print $2}' | sed 's/ //g')
reads=${reads_tmp:-0}
aligned_tmp=$(grep "^Total reads" outputs/tophat_mmusculus/accepted_hits.bam.stats | awk '{print $3}' | sed 's/ .*//g')
aligned=${aligned_tmp:-0}
failed_tmp=$(grep "^Total reads" outputs/tophat_mmusculus/unmapped.bam.stats | awk '{print $3}' | sed 's/ .*//g')
failed=${failed_tmp:-0}
rpm_tmp=$(perl -e "printf(1000000 / ${aligned})" 2>/dev/null)
rpm=${rpm_tmp:-0}
stat_string=$(printf "hpgl0850-trimmed,mmusculus,%s,%s,%s,%s,%s,accepted_hits.count.xz" "${original_reads}" "${reads}" "${aligned}" "${failed}" "$rpm")
echo "$stat_string" >> outputs/tophat_stats.csv

4.5.2 HTSeq for tophat in M.musculus

I usually use htseq’s default parameters against the perfectly paired reads when making count tables from tophat mappings, in this data set, that might be a bad choice.

It is worth noting that when possible (as it is for the mouse genome), I do separate countings for the following feature types: 3’ UTRs, 5’ UTRs, exons, rRNA, snRNA, snoRNA, pseudogenes, lincRNA, miRNA, NMD targets, introns, antisense transcripts, retained introns, exons, genes, and all.

But I will show only the all script, as in most cases that is all we care about.

The following is preprocessing/hpgl0850/scripts/32htall_mmusculus-hpgl0840trimmed.sh

## Counting the number of hits in outputs/tophat_mmusculus/accepted_paired.bam for each feature found in /cbcbhomes/abelew/libraries/genome/mmusculus_all.gtf
## Is this stranded? yes.  The defaults of htseq are:
##  --order=name --idattr=gene_id --minaqual=10 --type=exon --stranded=yes --mode=union

htseq-count -q -f bam -s yes -i ID  -t exon \
  outputs/tophat_mmusculus/accepted_paired.bam /cbcbhomes/abelew/libraries/genome/mmusculus_all.gtf \
  1>outputs/htseq/accepted_paired_all.count 2>outputs/tophat_mmusculus/accepted_paired_htseq.err && \
    xz -f -9e outputs/htseq/accepted_paired_all.count

4.6 Tophat L.major

This is the same as above, but using my indexes/genome for Leishmania major from the TriTrypDb version 27.

4.7 Bowtie2 L.major

Here is where things get interesting and I started to figure out what is wrong.

I did a little reading and found the following about how bowtie2 scores multi-hit reads:

Taken from: http://biofinysics.blogspot.com/2014/05/how-does-bowtie2-assign-mapq-scores.html#bt2expt

“In bowtie2, the best a true multiread (AS=XS) can get is MAPQ=1 regardless of how low or high its multiplicity. This occurs when there are 0 or 1 mismatches over perfect base calls in the read, or when AS=XS goes down to -6. When there are 2-5 mismatches over perfect base calls (or the AS=XS <= -12 —- i.e. -12 to -30.6), the MAPQ becomes 0.

If someone wanted to exclude “true multireads” from their data set, using MAPQ >= 2 would work. This would also exclude any uniquely mapping reads with >=4 mismatches over high quality bases. In terms of high quality bases and unireads, MAPQ >= 3 allows up to 3 mismatches, MAPQ >= 23 allows up to 2 mismatches, MAPQ >= 40 allows up to 1 mismatch, and MAPQ >= 42 allows 0 mismatches. There will also be other “maxireads” in most or all of these sets."

This is relevant because of the following gem in the newest version of the htseq-count documentation:

http://www-huber.embl.de/HTSeq/doc/count.html

“skip all reads with alignment quality lower than the given minimum value (default: 10 — Note: the default used to be 0 until version 0.5.4.)”

With that in mind, take a look at the bowtie2 script, The following is preprocessing/scripts/hpgl0850/15bt2_lmajor-hpgl0850trimmed.sh

## This is a bowtie2 alignment of  -1 hpgl0850_forward-trimmed.fastq -2 hpgl0850_reverse-trimmed.fastq  against
## /cbcbhomes/abelew/libraries/genome/indexes/lmajor using arguments:  --very-sensitive .
## This jobs depended on: 1319885.cbcbsrv.umiacs.umd.edu

mkdir -p outputs/bowtie2_lmajor && sleep 10 && bowtie2 -x /cbcbhomes/abelew/libraries/genome/indexes/lmajor  --very-sensitive  \
  -p 1 \
  -q   -1 hpgl0850_forward-trimmed.fastq -2 hpgl0850_reverse-trimmed.fastq  \
  --un outputs/bowtie2_lmajor/hpgl0850-trimmed_unaligned_lmajor.fastq \
  --al outputs/bowtie2_lmajor/hpgl0850-trimmed_aligned_lmajor.fastq \
  -S outputs/bowtie2_lmajor/hpgl0850-trimmed.sam \
  2>outputs/bowtie2_lmajor/hpgl0850-trimmed.err \
  1>outputs/bowtie2_lmajor/hpgl0850-trimmed.out

4.7.1 Converting to an indexed sorted bam alignment

I hate sam files, so I immediately convert them before counting:

Just like with tophat, I have an extra step to extract the properly paired aligned reads.

The following is preprocessing/scripts/hpgl0850/19s2b_bt2_lmajor-hpgl0850trimmed.shell

## Converting the text sam to a compressed, sorted, indexed bamfile.
## Also printing alignment statistics to outputs/bowtie2_lmajor/hpgl0840-trimmed.bam.stats
## This job depended on: 1323883.cbcbsrv.umiacs.umd.edu
samtools view -u -t /cbcbhomes/abelew/libraries/genome/lmajor.fasta \
  -S outputs/bowtie2_lmajor/hpgl0840-trimmed.sam -o outputs/bowtie2_lmajor/hpgl0840-trimmed.bam 1>outputs/bowtie2_lmajor/hpgl0840-trimmed.bam.out 2>&1
samtools sort -l 9 outputs/bowtie2_lmajor/hpgl0840-trimmed.bam -o outputs/bowtie2_lmajor/hpgl0840-trimmed-sorted.bam 2>outputs/bowtie2_lmajor/hpgl0840-trimmed-sorted.bam.out 1>&2 && \
  rm outputs/bowtie2_lmajor/hpgl0840-trimmed.bam && \
  rm outputs/bowtie2_lmajor/hpgl0840-trimmed.sam && \
  mv outputs/bowtie2_lmajor/hpgl0840-trimmed-sorted.bam outputs/bowtie2_lmajor/hpgl0840-trimmed.bam && \
  samtools index outputs/bowtie2_lmajor/hpgl0840-trimmed.bam && \
  samtools view -b -f 2 -o outputs/bowtie2_lmajor/hpgl0840-trimmed.bam_paired.bam outputs/bowtie2_lmajor/hpgl0840-trimmed.bam && \
  samtools index outputs/bowtie2_lmajor/hpgl0840-trimmed.bam_paired.bam
bamtools stats -in outputs/bowtie2_lmajor/hpgl0840-trimmed.bam 2>outputs/bowtie2_lmajor/hpgl0840-trimmed.bam.stats 1>&2 && \
  bamtools stats -in outputs/bowtie2_lmajor/hpgl0840-trimmed.bam_paired.bam 2>outputs/bowtie2_lmajor/hpgl0840-trimmed.bam_paired.bam.stats 1>&2

4.7.2 Counting the bowtie2 results

In the following block, I am actually going to print my test data set’s script. In the test data, I took just the first 25,000 reads. I will follow that up with a small section of the count table generated from it.

I think this will show where the wheels fell off the car.

There is an important caveat related to the note I pasted above. In the following script I added the ‘-a 1’ argument, which tells htseq to count all pairs with alignment qualities >= 1 rather than the default of 10.

## Counting the number of hits in outputs/bowtie2_lmajor/test-trimmed.bam for each feature found in /cbcbhomes/abelew/libraries/genome/lmajor.gff
## Is this stranded? yes.  The defaults of htseq are:
##  --order=name --idattr=gene_id --minaqual=10 --type=exon --stranded=yes --mode=union

htseq-count -q -f bam -s yes  -i ID  -t exon -a 1 \
  outputs/bowtie2_lmajor/test-trimmed.bam /cbcbhomes/abelew/libraries/genome/lmajor.gff \
  1>outputs/bowtie2_lmajor/test-trimmed.count 2>outputs/bowtie2_lmajor/test-trimmed_htseq.err && \
    xz -f -9e outputs/bowtie2_lmajor/test-trimmed.count

Here are a few representative lines from the resulting count table, the main thing to recall is that in most of our data sets, the L. major rRNA genes are mostly excluded, and the genes surrounding the rRNA loci are usually 3-10x higher in counts:

exon_LmjF.27.2590-1     1
exon_LmjF.27.2600-1     0
exon_LmjF.27.2610-1     0
exon_LmjF.27.2620-1     0
exon_LmjF.27.2630-1     0
exon_LmjF.27.2640-1     0
exon_LmjF.27.2650-1     0
exon_LmjF.27.2660-1     0
exon_LmjF.27.ncRNA1-1   0
exon_LmjF.27.rRNA.01-1  1298
exon_LmjF.27.rRNA.02-1  1388
exon_LmjF.27.rRNA.03-1  1365
exon_LmjF.27.rRNA.04-1  1437
exon_LmjF.27.rRNA.05-1  1343
exon_LmjF.27.rRNA.06-1  1410

Oh crap. I did not see these in my first run because htseq is dropping them all as (I think) many of them are multi-hitting to a significant subset of the 52 tandem rRNA loci in the L. major genome.