S.cerevisiae 2016: Preprocessing RNAseq samples.

1 Preprocessing some S. cerevisiae samples.

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

2.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
  cyoa --task rnaseq --method fastq --input ${i}_forward.fastq.gz
  echo "Finishing ${i}"
  cd ${start}
done

2.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 a few (pseudo)aligners: bowtie2, tophat, and kallisto.

#!/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"
  cyoa --task rnaseq --method kallisto --input "${i}_forward-trimmed.fastq.gz" --species scerevisiae
  cyoa --task rnaseq --method tophat --input "${i}_forward-trimmed.fastq.gz" --species scerevisiae
  cyoa --task rnaseq --method bt2 --input "${i}_forward-trimmed.fastq.gz" --species scerevisiae
  echo "Finishing ${i}"
  cd "${start}" || exit
done

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

3.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/hpgl0564/scripts/00fqc-hpgl0564.sh

#!/usr/bin/env bash
#PBS -V -S /usr/bin/bash -q workstation
#PBS -d /cbcb/nelsayed-scratch/atb/rnaseq/scerevisiae_cbf5v1/preprocessing/hpgl0564
#PBS -N fqc-hpgl0564 -l mem=6gb -l walltime=10:00:00 -l ncpus=8
#PBS -o /cbcb/nelsayed-scratch/atb/rnaseq/scerevisiae_cbf5v1/preprocessing/hpgl0564/outputs/qsub/fqc-hpgl0564.qsubout -j oe -V -m n
echo "####Started /cbcb/nelsayed-scratch/atb/rnaseq/scerevisiae_cbf5v1/preprocessing/hpgl0564/scripts/00fqc-hpgl0564.sh at $(date)" >> outputs/log.txt
cd /cbcb/nelsayed-scratch/atb/rnaseq/scerevisiae_cbf5v1/preprocessing/hpgl0564 || 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 hpgl0564_forward.fastq.gz \
  2>outputs/fastqc.out 1>&2

## The following lines give status codes and some logging
echo $? > outputs/status/fqc-hpgl0564.status
echo "###Finished ${PBS_JOBID} 00fqc-hpgl0564.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

3.1.1 Fastqc statistics

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

The following is preprocessing/hpgl0480/scripts/01fqcst_forward-hpgl0564.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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_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/hpgl0564_forward.gz_fastqc/fastqc_data.txt | awk -F '\\t' '{print $2}')
kmer_content=${kmer_content_tmp:-0}

stat_string=$(printf "hpgl0564,%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

3.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/hpgl0564/scripts/05trim-hpgl0564.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 hpgl0564_forward.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 "hpgl0564_forward.fastq.gz" ]]; then
  if [[ -r "sequences/hpgl0564_forward.fastq.xz" ]]; then
    mv sequences/hpgl0564_forward.fastq.xz . && pxz -d hpgl0564_forward.fastq.xz && pigz hpgl0564_forward.fastq
  else
    echo "Missing files. Did not find hpgl0564_forward.fastq.gz nor sequences/hpgl0564_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 hpgl0564_forward.fastq.gz hpgl0564_forward-trimmed_paired.fastq.gz hpgl0564_forward-trimmed_unpaired.fastq.gz \
    ILLUMINACLIP:/cbcbhomes/abelew/libraries/adapters.fa:2:20:4 SLIDINGWINDOW:4:25 \
    1>outputs/hpgl0564-trimomatic.out 2>&1
excepted=$(grep "Exception" outputs/hpgl0564-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 hpgl0564_forward.fastq.gz hpgl0564_forward-trimmed_paired.fastq.gz hpgl0564_forward-trimmed_unpaired.fastq.gz SLIDINGWINDOW:4:25 \
    1>outputs/hpgl0564-trimomatic.out 2>&1
fi
sleep 10
mv hpgl0564_forward-trimmed_paired.fastq.gz hpgl0564_forward-trimmed.fastq.gz

3.2.1 Trimomatic statistics

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

The following is preprocessing/hpgl0564/scripts/06trimst-hpgl0564.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/hpgl0564-trimomatic.out | awk '{print $4}')
total_reads=${total_reads_tmp:-0}
surviving_both_tmp=$(grep "^Input Read Pairs" outputs/hpgl0564-trimomatic.out | awk '{print $7}')
surviving_both=${surviving_both_tmp:-0}
surviving_forward_tmp=$(grep "^Input Read Pairs" outputs/hpgl0564-trimomatic.out | awk '{print $12}')
surviving_forward=${surviving_forward_tmp:-0}
surviving_reverse_tmp=$(grep "^Input Read Pairs" outputs/hpgl0564-trimomatic.out | awk '{print $17}')
surviving_reverse=${surviving_reverse_tmp:-0}
dropped_reads_tmp=$(grep "^Input Read Pairs" outputs/hpgl0564-trimomatic.out | awk '{print $20}')
dropped_reads=${dropped_reads_tmp:-0}

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

3.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/hpgl0564/scripts/02biop-hpgl0564trimmed.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 hpgl0564_forward-trimmed.fastq | read_fastq -i - -e base_33 |\
 plot_scores -T 'Quality Scores' -t svg -o outputs/biopieces/hpgl0564_forward-trimmed_quality_scores.svg |\
 plot_nucleotide_distribution -T 'NT. Distribution' -t svg -o outputs/biopieces/hpgl0564_forward-trimmed_ntdist.svg |\
 plot_lendist -T 'Length Distribution' -k SEQ_LEN -t svg -o outputs/biopieces/hpgl0564_forward-trimmed_lendist.svg |\
 analyze_gc |\
     bin_vals -b 20 -k 'GC%' |\
     plot_distribution -k 'GC%_BIN' -t svg -o outputs/biopieces/hpgl0564_forward-trimmed_gcdist.svg |\
 analyze_gc |\
     mean_vals -k 'GC%' -o outputs/biopieces/hpgl0564_forward-trimmed_gc.txt |\
 count_records -o outputs/biopieces/hpgl0564_forward-trimmed_count.txt -x

3.4 Kallisto Yeast

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/hpgl0564/scripts/30kall_scerevisiae-hpgl0564trimmed.sh

## This is a kallisto pseudoalignment of hpgl0564_forward-trimmed.fastq.gz against
## /cbcbhomes/abelew/libraries/genome/indexes/scerevisiae.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_scerevisiae/abundance.tsv > outputs/kallisto_scerevisiae/abundance.count
##   awk '{printf("%s %s\n", $1, $4)}' outputs/kallisto_scerevisiae/abundance.tsv > outputs/kallisto_scerevisiae/abundance.count
## because kallisto is so fast

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

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

3.5 Tophat Yeast

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/hpgl0564/scripts/31th_scerevisiae-hpgl0564trimmed.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_scerevisiae && tophat  -g 1 --b2-very-sensitive  \
  -G /cbcbhomes/abelew/libraries/genome/scerevisiae.gtf \
  -p 4 -o outputs/tophat_scerevisiae \
  /cbcbhomes/abelew/libraries/genome/indexes/scerevisiae \
  hpgl0564_forward-trimmed.fastq.gz 2>outputs/tophat.out 1>&2 && \
 samtools sort -l 9 -n outputs/tophat_scerevisiae/accepted_hits.bam outputs/tophat_scerevisiae/accepted_sorted && \
 mv outputs/tophat_scerevisiae/accepted_sorted.bam outputs/tophat_scerevisiae/accepted_hits.bam && \
 samtools index outputs/tophat_scerevisiae/accepted_hits.bam && \
 samtools sort -l 9 -n outputs/tophat_scerevisiae/unmapped.bam outputs/tophat_scerevisiae/unmapped_sorted && \
 mv outputs/tophat_scerevisiae/unmapped_sorted.bam outputs/tophat_scerevisiae/unmapped.bam && \
 samtools index outputs/tophat_scerevisiae/unmapped.bam

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

3.5.1 Tophat statistics

I generate a similar csv file for the tophat alignments:

The following is preprocessing/hpgl0564/scripts/33tpstats_scerevisiae-hpgl0564trimmed.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_scerevisiae/accepted_hits.bam 2>outputs/tophat_scerevisiae/accepted_hits.bam.stats 1>&2 && bamtools stats < outputs/tophat_scerevisiae/unmapped.bam 2>outputs/tophat_scerevisiae/unmapped.bam.stats 1>&2
original_reads_tmp=$(grep "^Input Reads" outputs/hpgl0564-trimmed-trimomatic.out | awk '{print $3}' | sed 's/ //g')
original_reads=${original_reads_tmp:-0}
reads_tmp=$(grep "^reads_in " outputs/tophat_scerevisiae/prep_reads.info | awk -F= '{print $2}' | sed 's/ //g')
reads=${reads_tmp:-0}
aligned_tmp=$(grep "^Total reads" outputs/tophat_scerevisiae/accepted_hits.bam.stats | awk '{print $3}' | sed 's/ .*//g')
aligned=${aligned_tmp:-0}
failed_tmp=$(grep "^Total reads" outputs/tophat_scerevisiae/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 "hpgl0564-trimmed,scerevisiae,%s,%s,%s,%s,%s,accepted_hits.count.xz" "${original_reads}" "${reads}" "${aligned}" "${failed}" "$rpm")
echo "$stat_string" >> outputs/tophat_stats.csv

3.5.2 HTSeq for tophat in S.cerevisiae

I usually use htseq’s default parameters against the perfectly paired reads when making count tables from tophat mappings. Sometimes this is a bad idea, but in yeast I think it is ok.

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 for yeast I just do the set of all annotations, as in most cases that is all we care about.

The following is preprocessing/hpgl0564/scripts/32htall_scerevisiae-hpgl0564trimmed.sh

## Counting the number of hits in outputs/tophat_scerevisiae/accepted_paired.bam for each feature found in /cbcbhomes/abelew/libraries/genome/scerevisiae_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_scerevisiae/accepted_paired.bam /cbcbhomes/abelew/libraries/genome/scerevisiae_all.gtf \
  1>outputs/htseq/accepted_paired_all.count 2>outputs/tophat_scerevisiae/accepted_paired_htseq.err && \
    xz -f -9e outputs/htseq/accepted_paired_all.count

3.6 Bowtie2

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/hpgl0564/15bt2_scerevisiae-hpgl0564trimmed.sh

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

mkdir -p outputs/bowtie2_scerevisiae && sleep 10 && bowtie2 -x /cbcbhomes/abelew/libraries/genome/indexes/scerevisiae  --very-sensitive  \
  -p 1 \
  -q   hpgl0564_forward-trimmed.fastq  \
  --un outputs/bowtie2_scerevisiae/hpgl0564-trimmed_unaligned_scerevisiae.fastq \
  --al outputs/bowtie2_scerevisiae/hpgl0564-trimmed_aligned_scerevisiae.fastq \
  -S outputs/bowtie2_scerevisiae/hpgl0564-trimmed.sam \
  2>outputs/bowtie2_scerevisiae/hpgl0564-trimmed.err \
  1>outputs/bowtie2_scerevisiae/hpgl0564-trimmed.out

3.6.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/hpgl0564/19s2b_bt2_scerevisiae-hpgl0564trimmed.shell

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

3.6.2 Counting the bowtie2 results

## Counting the number of hits in outputs/bowtie2_scerevisiae/hpgl0564-trimmed.bam for each feature found in /cbcbhomes/abelew/libraries/genome/scerevisiae.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 \
  outputs/bowtie2_scerevisiae/hpgl0564-trimmed.bam /cbcbhomes/abelew/libraries/genome/scerevisiae.gff \
  1>outputs/bowtie2_scerevisiae/hpgl0564-trimmed.count 2>outputs/bowtie2_scerevisiae/hpgl0564-trimmed_htseq.err && \
    xz -f -9e outputs/bowtie2_scerevisiae/hpgl0564-trimmed.count