Page Comparison

Exercise #4: Bowtie2 global alignment - Vibrio cholerae RNA-seq

While we have focused on aligning eukaryotic data, the same tools can be used with prokaryotic data. The major differences are less about the underlying data and much more about the external/public databases that store and distribute reference data. If we want to study a prokaryote, the reference data is usually downloaded from a resource like GenBank.  

In this exercise, we will use some RNA-seq data from Vibrio cholerae, published on GEO here, and align it to a reference genome.

Overview of Vibrio cholerae alignment workflow with Bowtie2

Alignment of this prokaryotic data follows the workflow below. Here we will concentrate on steps 1 and 2.

1. Prepare the vibCho reference index for bowtie2 from GenBank records
2. Align reads using bowtie2, producing a SAM file
3. Convert the SAM file to a BAM file (samtools view) 
4. Sort the BAM file by genomic location (samtools sort)
5. Index the BAM file (samtools index)
6. Gather simple alignment statistics (samtools flagstat and samtools idxstat)

Obtaining the GenBank records

First prepare a directory to work in, and change to it:

mkdir -p $SCRATCH/core_ngs/references/vibCho
cd $SCRATCH/core_ngs/references/vibCho

V. cholerae has two chromosomes. We download each separately.

1. Navigate to http://www.ncbi.nlm.nih.gov/nuccore/NC_012582
   • click on the Send to down arrow (top right of page)
     • select Complete Record
     • select File as Destination, and Format FASTA
     • click Create File
       
   • in the Opening File dialog, select Save File then OK
     • Save the file on your local computer as NC_012582.fa
2. Back on the main http://www.ncbi.nlm.nih.gov/nuccore/NC_012582 page
   • click on the Send to down arrow (top right of page)
     • select Complete Record
     • select File as Destination, and Format GFF3
     • click Create File
       
   • in the Opening File dialog, select Save File then OK
     • Save the file on your local computer as NC_012582.gff3
3. Repeat steps 1 and 2 for the 2nd chromosome
   
   • NCBI URL is http://www.ncbi.nlm.nih.gov/nuccore/NC_012583
   • use NC_012583 as the filename prefix for the files you save
   • you should now have 4 files:
     • NC_012582.fa, NC_012582.gff3
     • NC_012583.fa, NC_012583.gff3
       
4. Transfer the files from your local computer to TACC
   • to the ~/scratch/core_ngs/references/vibCho directory created above
     
     • On a Mac or Windows with a WSL shell, use scp from your laptop
       
     • On Windows, use the pscp.exe PuTTy tool
     • See Copying files between TACC and your laptop

mkdir -p $SCRATCH/core_ngs/references/vibCho
cd $SCRATCH/core_ngs/references/vibCho
cp $CORENGS/ncbi/vibCho/NC_* .

Once you have the 4 files locally in your $SCRATCH/core_ngs/references/vibCho directory, combine them using cat:

cd $SCRATCH/core_ngs/references/vibCho
cat NC_01258[23].fa   > vibCho.O395.fa
cat NC_01258[23].gff3 > vibCho.O395.gff3

Now we have a reference sequence file that we can use with the bowtie2 reference builder, and ultimately align sequence data against.

Introducing bowtie2

First make sure you're in an idev session:

idev -m 180 -p normal -A UT-2015-05-18 -N 1 -n 68 

Go ahead and load the bowtie2 module so we can examine some help pages and options.

module biocontainers
module load bowtie

 Now that it's loaded, check out the options. There are a lot of them! In fact for the full range of options and their meaning, Google "Bowtie2 manual" and bring up that page (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml). The Table of Contents is several pages long! Ouch!

This is the key to using bowtie2 - it allows you to control almost everything about its behavior, which make it the go-to aligner for specialized alignment tasks (e.g. aligning miRNA or other small reads). But it also makes it is much more challenging to use than bwa – and it's easier to screw things up too!

Building the bowtie2 vibCho index

Before the alignment, of course, we've got to build a bowtie2- specific index using bowtie2-build. Go ahead and check out its options. Unlike for the aligner itself, we only need to worry about a few things here:

• reference_in file is just the vibCho.O395.fa FASTA we built from GenBank records
• bt2_index_base is the prefix of all the bowtie2-build output file

Here, to build the reference index for alignment, we only need the FASTA file. (This is not always true - extensively spliced transcriptomes requires splice junction annotations to align RNA-seq data properly.)

First create a directory specifically for the bowtie2 index, then build the index using bowtie-build.

mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cd $SCRATCH/core_ngs/references/bt2/vibCho

# Symlink to the fasta file you created
ln -sf $SCRATCH/core_ngs/references/vibCho.O395.fa
# or, to catch up:
ln -sf $CORENGS/references/vibCho.O395.fa

bowtie2-build vibCho.O395.fa vibCho.O395

This should also go pretty fast. You can see the resulting files using ls like before.

Performing the bowtie2 alignment

We'll set up a new directory to perform the V. cholerae data alignment. But first make sure you have the FASTQ file to align and the vibCho bowtie2 index:

# Get the FASTQ to align
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz $SCRATCH/core_ngs/alignment/fastq/

# Set up the bowtie2 index
mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/idx/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho

Make sure you're in an idev session with the bowtie2 BioContainers module loaded:

idev -p normal -m 120 -A UT-2015-05-18 -N 1 -n 68
module load biocontainers
module load bowtie

Now set up a directlry to do this alignment, with symbolic links to the bowtie2 index directory and the directory containing the FASTQ to align:

mkdir -p $SCRATCH/core_ngs/alignment/vibCho
cd $SCRATCH/core_ngs/alignment/vibCho
ln -sf ../../references/bt2/vibCho
ln -sf ../../alignment/fastq fq

We'll be aligning the V. cholerae reads now in ./fq/cholera_rnaseq.fastq.gz (how many sequences does it contain?)

Note that here the data is from standard mRNA sequencing, meaning that the DNA fragments are typically longer than the reads. There is likely to be very little contamination that would require using a local rather than global alignment, or many other pre-processing steps (e.g. adapter trimming). Thus, we will run bowtie2 with default parameters, omitting options other than the input, output, and reference index. This performs a global alignment.

As you can tell from looking at the bowtie2 help message, the general syntax looks like this:

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

So execute this bowtie2 global, single-end alignment command:

cd $SCRATCH/core_ngs/alignment/vibCho
bowtie2 -x vibCho/vibCho.O395 -U fq/cholera_rnaseq.fastq.gz -S cholera_rnaseq.sam 2>&1 | tee aln.log

Notes:

• -x  vibCho/vibCho.O395.fa – prefix path of index files
• -U fq/cholera_rnaseq.fastq.gz – FASTQ file for single-end (Unpaired) alignment
• -S cholera_rnaseq.sam – tells bowtie2 to report alignments in SAM format to the specified file
• 2>&1 redirects standard error to standard output
  • while the alignment data is being written to the cholera_rnaseq.sam file, bowtie2 will report its progress to standard error.
• | tee aln.log takes the bowtie2 progress output and pipes it to the tee program
  • tee takes its standard input and writes it to the specified file and also to standard output
  • that way, you can see the progress output now, but also save it to review later (or supply to MultiQC)

Since the FASTQ file is not large, this should not take too long, and you will see progress output like this:

89006 reads; of these:
  89006 (100.00%) were unpaired; of these:
    20675 (23.23%) aligned 0 times
    38226 (42.95%) aligned exactly 1 time
    30105 (33.82%) aligned >1 times
76.77% overall alignment rate

When the job is complete you should have a cholera_rnaseq.sam file that you can examine using whatever commands you like.  Remember, to further process it downstream, you should create a sorted, indexed BAM file from this SAM output.

Exercise #5: Bowtie2 local alignment - Human microRNA-seq

Now we're going to switch over to a different aligner that was originally designed for very short reads and is frequently used for RNA-seq data. Accordingly, we have prepared another test microRNA-seq dataset for you to experiment with (not the same one you used cutadapt on). This data is derived from a human H1 embryonic stem cell (H1-hESC) small RNA dataset generated by the ENCODE Consortium – its about a half million reads.

However, there is a problem!  We don't know (or, well, you don't know) what the adapter structure or sequences were. So, you have a bunch of 36 base pair reads, but many of those reads will include extra sequence that can impede alignment – and we don't know where! We need an aligner that can find subsections of the read that do align, and discard (or "soft-clip") the rest away – an aligner with a local alignment mode. Bowtie2 is just such an aligner.

Overview miRNA alignment workflow with bowtie2

If the adapter structure were known, the normal workflow would be to first remove the adapter sequences with cutadapt. Since we can't do that, we will instead perform a local alignment of the single-end miRNA sequences using bowtie2. This workflow has the following steps:

1. Prepare the mirbase v20 reference index for bowtie2 (one time) using bowtie2-build 
2. Perform local alignment of the R1 reads with bowtie2, producing a SAM file directly
3. Convert the SAM file to a BAM file (samtools view)
4. Sort the BAM file by genomic location (samtools sort)
5. Index the BAM file (samtools index)
6. Gather simple alignment statistics (samtools flagstat and samtools idxstat)

This looks so much simpler than bwa – only one alignment step instead of three! We'll see the price for this "simplicity" in a moment...

As before, we will just do the alignment steps leave samtools for the next section.

Mirbase is a collection of all known microRNAs in all species (and many speculative miRNAs). We will use the human subset of that database as our alignment reference.  This has the advantage of being significantly smaller than the human genome, while likely containing almost all sequences likely to be detected in a miRNA sequencing run.

These are the four reference genomes we will be using today, with some information about them (and here is information about many more genomes):

Building the bowtie2 mirbase index

Before the alignment, of course, we've got to build a mirbase index using bowtie2-build (go ahead and check out its options). Unlike for the aligner itself, we only need to worry about a few things here:

bowtie2-build <reference_in> <bt2_index_base>

• reference_in file is just the FASTA file containing mirbase v20 sequences
• bt2_index_base is the prefix of where we want the files to go

 Following what we did earlier for BWA indexing, namely move our FASTA into place, create the index directory, and establish our symbolic links.

mkdir -p $WORK/core_ngs/references/bt2/mirbase.v20
mv $WORK/core_ngs/references/hairpin_cDNA_hsa.fa $WORK/core_ngs/references/fasta
cd $WORK/core_ngs/references/bt2/mirbase.v20
ln -s -f ../../fasta/hairpin_cDNA_hsa.fa
ls -la 
module load perl
module load bowtie/2.2.0

Now build the mirbase index with bowtie2-build like we did for the V. cholerae index:

bowtie2-build hairpin_cDNA_hsa.fa hairpin_cDNA_hsa.fa

That was very fast!  It's because the mirbase reference genome is so small compared to what programs like this are used to dealing with, which is the human genome (or bigger).  You should see the following files:

hairpin_cDNA_hsa.fa
hairpin_cDNA_hsa.fa.1.bt2
hairpin_cDNA_hsa.fa.2.bt2
hairpin_cDNA_hsa.fa.3.bt2
hairpin_cDNA_hsa.fa.4.bt2
hairpin_cDNA_hsa.fa.rev.1.bt2
hairpin_cDNA_hsa.fa.rev.2.bt2

Performing the bowtie2 local alignment

Now, we're ready to actually try to do the alignment.  Remember, unlike BWA, we actually need to set some options depending on what we're after. Some of the important options for bowtie2 are:

To decide how we want to go about doing our alignment, check out the file we're aligning with less:

cd $SCRATCH/core_ngs/alignment
less fastq/human_mirnaseq.fastq.gz

Lots of reads have long strings of A's, which must be an adapter or protocol artifact.  Even though we see how we might be able to fix it using some tools we've talked about, what if we had no idea what the adapter sequence was, or couldn't use cutadapt or other programs to prepare the reads?

 In that case, we need a local alignment where the seed length smaller than the expected insert size. Here, we are interested in finding any sections of any reads that align well to a microRNA, which are between 16 and 24 bases long, with most 20-22. So an acceptable alignment should have at least 16 matching bases, but could have more.

If we're also interested in detecting miRNA SNPs, we might want to allow a mismatch in the seed. So, a good set of options might look something like this:

-N 1 -L 16 --local

Let's make a link to the mirbase index directory to make our command line simpler:

cd $SCRATCH/core_ngs/alignment
ln -s -f $WORK/core_ngs/references/bt2/mirbase.v20 mb20

Putting this all together we have a command line that looks like this.

bowtie2 --local -N 1 -L 16 -x mb20/hairpin_cDNA_hsa.fa -U fastq/human_mirnaseq.fastq.gz -S human_mirnaseq.sam

Create a commands file called bt2.cmds with this task definition then generate and submit a batch job for it (time 1 hour, development queue).

Use nano to create the bt2.cmds file. Then:

launcher_creator.py -n bt2 -j bt2.cmds -t 01:00:00 -A UT-2015-05-18
sbatch bt2.slurm; showq -u

#copy from Amelia's scratch:
cp /scratch/01786/awh394/core_ngs.test/alignment/human_mirnaseq.sam .

When the job is complete you should have a human_mirnaseq.sam file that you can examine using whatever commands you like. An example alignment looks like this.

TUPAC_0037_FC62EE7AAXX:2:1:2607:1430#0/1  0  hsa-mir-302b  50  22 3S20M13S * 0 0
    TACGTGCTTCCATGTTTTANTAGAAAAAAAAAAAAG  ZZFQV]Z[\IacaWc]RZIBVGSHL_b[XQQcXQcc
    AS:i:37 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:16G3       YT:Z:UU

Notes:

• This is one alignment record, although it has been broken up below for readability.
• This read mapped to the mature microRNA sequence hsa-mir-302b, starting at base 50 in that contig.
• Notice the CIGAR string is 3S20M13S, meaning that 3 bases were soft clipped from one end (3S), and 13 from the other (13S).
  
  • If we did the same alignment using either bowtie2 --end-to-end mode, or using bwa aln as in Exercise #1, very little of this file would have aligned.
• The 20M part of the CIGAR string says there was a block of 20 read bases that mapped to the reference.
  • If we had not lowered the seed parameter of bowtie2 from its default of 22, we would not have found many of the alignments like this one that only matched for 20 bases.

Such is the nature of bowtie2 – it it can be a powerful tool to sift out the alignments you want from a messy dataset with limited information, but doing so requires careful tuning of the parameters, which can take quite a few trials to figure out.

Exercise: About how many records in human_mirnaseq.sam represent aligned reads?

grep -P -v '^@' human_mirnaseq.sam | cut -f 3 | grep -v '*' | wc -l

Use sort and uniq to create a histogram of mapping qualities

The mapping quality score is in field 5 of the human_mirnaseq.sam file. We can do this to pull out only that field:

grep -P -v '^@' human_mirnaseq.sam | cut -f 5 | head

We will use the uniq create a histogram of these values. The first part of the --help for uniq says:

Usage: uniq [OPTION]... [INPUT [OUTPUT]]
Filter adjacent matching lines from INPUT (or standard input),
writing to OUTPUT (or standard output).

With no options, matching lines are merged to the first occurrence.
Mandatory arguments to long options are mandatory for short options too.
  -c, --count           prefix lines by the number of occurrences

To create a histogram, we want to organize all equal mapping quality score lines into an adjacent block, then use uniq -c option to count them. The sort -n command does the sorting into blocks (-n means numerical sort). So putting it all together, and piping the output to a pager just in case, we get:

grep -P -v '^@' human_mirnaseq.sam | cut -f 5 | sort -n | uniq -c | more

Exercise: What is the flaw in this "program"?

Exercise #6: BWA-MEM - Human mRNA-seq

After bowtie2 came out with a local alignment option, it wasn't long before bwa developed its own local alignment algorithm called BWA-MEM (for Maximal Exact Matches), implemented by the bwa mem command. bwa mem has the following advantages:

• It incorporates a lot of the simplicity of using bwa with the complexities of local alignment, enabling straightforward alignment of datasets like the mirbase data we just examined
• It can align different portions of a read to different locations on the genome
  • In a long RNA-seq experiment, reads will (at some frequency) span a splice junction themselves, or a pair of reads in a paired-end library will fall on either side of a splice junction. We want to be able to align reads that do this for many reasons, from accurate transcript quantification to novel fusion transcript discovery.

Thus, our last exercise will be the alignment of a human long RNA-seq dataset composed (by design) almost exclusively of reads that cross splice junctions.

bwa mem was made available when we loaded the bwa module, so take a look at its usage information. The most important parameters, similar to those we've manipulated in the past two sections, are the following:

...

There are many more parameters to control the scoring scheme and other details, but these are the most essential ones to use to get anything of value at all.

The test file we will be working with is just the R1 file from a paired-end total RNA-seq experiment, meaning it is (for our purposes) single-end. Go ahead and take a look at it, and find out how many reads are in the file.

cd $SCRATCH/core_ngs/alignment
ls fastq
gunzip -c fastq/human_rnaseq.fastq.gz | echo $((`wc -l`/4))

A word about real splice-aware aligners

...

# Request a 180 minute interactive node on the normal queue using our reservation
idev -m 120 -N 1 -A OTH21164 -r CoreNGS
idev -m 120 -N 1 -A TRA23004 -r CoreNGS

# Request a 120 minute idev node on the development queue 
idev -m 120 -N 1 -A OTH21164 -p development
idev -m 120 -N 1 -A TRA23004 -p development

# Using our reservation
sbatch --reseservation=CoreNGS <batch_file>.slurm

Exercise #3: PE alignment with BioITeam scripts

Now that you've done everything the hard way, let's see how to do run an alignment pipeline using a BWA alignment script maintained by the BioITeam,  /work/projects/BioITeam/common/script/align_bwa_illumina.sh. Type in the script name to see its usage.

align_bwa_illumina.sh 2022_06_10
Align Illumina SE or PE data with bwa. Produces a sorted, indexed,
duplicate-marked BAM file and various statistics files. Usage:

align_bwa_illumina.sh <aln_mode> <in_file> <out_pfx> <assembly> [ paired trim_sz trim_sz2 seq_fmt qual_fmt ]

Required arguments:
  aln_mode  Alignment mode, either global (bwa aln) or local (bwa mem).
  in_file   For single-end alignments, path to input sequence file.
            For paired-end alignments using fastq, path to the the R1
            fastq file which must contain the string 'R1' in its name.
            The corresponding 'R2' must have the same path except for 'R1'.
  out_pfx   Desired prefix of output files in the current directory.
  assembly  One of hg38, hg19, hg38, mm10, mm9, sacCer3, sacCer1, ce11, ce10,
            danRer7, hs_mirbase, mm_mirbase, or reference index prefix.
Optional arguments:
  paired    0 = single end alignment (default); 1 = paired end.
  trim_sz   Size to trim reads to. Default 0 (no trimming)
  trim_sz2  Size to trim R2 reads to for paired end alignments.
            Defaults to trim_sz
  seq_fmt   Format of sequence file (fastq, bam or scarf). Default is
            fastq if the input file has a '.fastq' extension; scarf
            if it has a '.sequence.txt' extension.
  qual_type Type of read quality scores (sanger, illumina or solexa).
            Default is sanger for fastq, illumina for scarf.
Environment variables:
  show_only  1 = only show what would be done (default not set)
  aln_args   other bowtie2 options (e.g. '-T 20' for mem, '-l 20' for aln)
  no_markdup 1 = don't mark duplicates (default 0, mark duplicates)
  run_fastqc 1 = run fastqc (default 0, don't run). Note that output
             will be in the directory containing the fastq files.
  keep       1 = keep unsorted BAM (default 0, don't keep)
  bwa_bin    BWA binary to use. Default bwa 0.7.x. Note that bwa 0.6.2
             or earlier should be used for scarf and other short reads.
  also: NUM_THREADS, BAM_SORT_MEM, SORT_THREADS, JAVA_MEM_ARG

Examples:
  align_bwa_illumina.sh local  ABC_L001_R1.fastq.gz my_abc hg38 1
  align_bwa_illumina.sh global ABC_L001_R1.fastq.gz my_abc hg38 1 50
  align_bwa_illumina.sh global sequence.txt old sacCer3 0 '' '' scarf solexa

There are lots of bells and whistles in the arguments, but the most important are the first few:

1. aln_mode – whether to perform a global or local alignment (the 1st argument must be one of those words)
   
   • global mode uses the bwa aln workflow as we did above
   • local mode uses the bwa mem command
2. in_file – full or relative path to the FASTQ file (just the R1 fastq if paired end). Can be compressed (.gz)
3. out_pfx – prefix for all the output files produced by the script. Should relate back to what the data is.
4. assembly – genome assembly to use.
   • there are pre-built indexes for some common eukaryotes (hg38, hg19, mm10, mm9, danRer7, sacCer3) that you can use
   • or provide a full path for a bwa reference index you have built somewhere
5. paired flag – 0 means single end (the default); 1 means paired end
6. trim_sz – if you want the FASTQ hard trimmed down to a specific length before alignment, supply that number here

We're going to run this script and a similar Bowtie2 alignment script, on the yeast data using the TACC batch system. In a new directory, copy over the commands and submit the batch job. We ask for 2 hours (-t 02:00:00) with 4 tasks/node (-w 4); since we have 4 commands, this will run on 1 compute node.

# Make sure you're not in an idev session by looking at the hostname
hostname
# If the hostname looks like "c455-004.ls6.tacc.utexas.edu", exit the idev session

# Copy over the Yeast data if needed
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/Sample_Yeast*.gz $SCRATCH/core_ngs/alignment/fastq/

# Make a new alignment directory for running these scripts
mkdir -p $SCRATCH/core_ngs/alignment/bwa_script
cd $SCRATCH/core_ngs/alignment/bwa_script
ln -sf ../fastq

# Copy the alignment commands file and submit the batch job
cd $SCRATCH/core_ngs/alignment/bwa_script 
cp $CORENGS/tacc/aln_script.cmds .

# Use -a TRA23004 below if OTH21164 hasn't been working for you...
launcher_creator.py -j aln_script.cmds -n aln_script -t 01:00:00 -w 4 -a OTH21164 -q normal
sbatch --reservation=CoreNGS-Thu aln_script.slurm 

# or
launcher_creator.py -j aln_script.cmds -n aln_script -t 01:00:00 -w 4 -a OTH21164 -q development
sbatch aln_script.slurm  

showq -u

While we're waiting for the job to complete, lets look at the aln_script.cmds file.

/work/projects/BioITeam/common/script/align_bwa_illumina.sh     global ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bwa_global sacCer3 1 50
/work/projects/BioITeam/common/script/align_bwa_illumina.sh     local  ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bwa_local  sacCer3 1
/work/projects/BioITeam/common/script/align_bowtie2_illumina.sh global ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bt2_global sacCer3 1 50
/work/projects/BioITeam/common/script/align_bowtie2_illumina.sh local  ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bt2_local  sacCer3 1

Notes:

• The 1st command performs a paired-end BWA global alignment (similar to above), but asks that the 100 bp reads be trimmed to 50 first.
  
  • we refer to the pre-built index for yeast by name: sacCer3
    • this index is located in the /work/projects/BioITeam/ref_genome/bwa/bwtsw/sacCer3/ directory
  • we provide the name of the R1 FASTQ file
    • because we request a PE alignment (the 1 argument) the script will look for a similarly-named R2 file.
  • all output files associated with this command will be named with the prefix bwa_global.
• The 2nd command performs a paired-end BWA local alignment.
  
  • all output files associated with this command will be named with the prefix bwa_local.
  • no trimming is requested because the local alignment should ignore 5' and 3' bases that don't match the reference genome
• The 3rd command performs a paired-end Bowtie2 global alignment.
  
  • the Bowtie2 alignment script has the same first arguments as the BWA alignment script.
  • all output files associated with this command will be named with the prefix bt2_global.
  • again, we specify that reads should first be trimmed to 50 bp.
• The 4th command performs a paired-end Bowtie2 local alignment.
  
  • all output files associated with this command will be named with the prefix bt2_local.
  • again, no trimming is requested for the local alignment.

Output files

This alignment pipeline script performs the following steps:

• Hard trims FASTQ, if optionally specified (fastx_trimmer)
  
• Performs the global or local alignment (here, a PE alignment)
  
  • BWA global: bwa aln the R1 and R2 separately, then bwa sampe to produce a SAM file
    
  • BWA local: call bwa mem with both R1 and R2 to produce a SAM file
    
  • Bowtie2 global: call bowtie2 in its default global (end-to-end) mode on both R1 and R2 to produce a SAM file
    
  • Bowtie2 local: call bowtie2 --local with both R1 and R2 to produce a SAM file
• Converts SAM to BAM (samtools view)
• Sorts the BAM (samtools sort)
• Marks duplicates (Picard MarkDuplicates)
• Indexes the sorted, duplicate-marked BAM (samtools index)
• Gathers statistics (samtools idxstats, samtools flagstat, plus a custom statistics script of Anna's)
• Removes intermediate files
  

There are a number of output files, with the most important being those desribed below.

1. <prefix>.align.log – Log file of the entire alignment process.
   • check the tail of this file to make sure the alignment was successful
2. <prefix>.sort.dup.bam – Sorted, duplicate-marked alignment file.
3. <prefix>.sort.dup.bam.bai – Index for the sorted, duplicate-marked alignment file
4. <prefix>.flagstat.txt – samtools flagstat output
   
5. <prefix>.idxstats.txt – samtools idxstats output
6. <prefix>.samstats.txt – Summary alignment statistics from Anna's stats script
7. <prefix>.iszinfo.txt – Insert size statistics (for paired-end alignments) from Anna's stats script

Verifying alignment success

The alignment log will have a  "I ran successfully" message at the end if all went well, and if there was an error, the important information should also be at the end of the log file. So you can use tail to check the status of an alignment. For example:

tail bwa_global.align.log

This will show something like:

------------------------------------------------------------------
..Done alignmentUtils.pl bamstats - 2022-06-10 12:59:05
.. samstats file 'bwa_global.samstats.txt' exists and is not empty - 2022-06-10 12:59:05
===============================================================================
## Cleaning up files (keep 0) - 2022-06-10 12:59:05
===============================================================================
ckRes 0 cleanup
===============================================================================
## All bwa alignment tasks completed successfully! - 2022-06-10 12:59:06
===============================================================================

Notice that success message:  "All bwa alignment tasks completed successfully!". It should only appear once in any successful alignment log.

When multiple alignment commands are run in parallel it is important to check them all, and you can use grep looking for part of the unique success message to do this. For example:

grep 'completed successfully!' *align.log | wc -l

If this command returns 4 (the number of alignment tasks we performed), all went well, and we're done.

But what if something went wrong? How can we tell which alignment task was not successful? You could tail the log files one by one to see which one(s) don't have the message, but you can also use a special grep option to do this work.

grep -L 'completed successfully' *.align.log

The -L option tells grep to only print the filenames that don't contain the pattern. Perfect! To see happens in the case of failure, try it on a file that doesn't contain that message:

grep -L 'completed successfully' aln_script.cmds

Checking alignment statistics

The <prefix>.samstats.txt statistics files produced by the alignment pipeline has a lot of good information in one place. If you look at bwa_global.samstats.txt you'll see something like this:

-----------------------------------------------
             Aligner:       bwa
     Total sequences:   1184360
        Total mapped:    539079 (45.5 %)
      Total unmapped:    645281 (54.5 %)
             Primary:    539079 (100.0 %)
           Secondary:
          Duplicates:    249655 (46.3 %)
          Fwd strand:    267978 (49.7 %)
          Rev strand:    271101 (50.3 %)
          Unique hit:    503629 (93.4 %)
           Multi hit:     35450 (6.6 %)
           Soft clip:
           All match:    531746 (98.6 %)
              Indels:      7333 (1.4 %)
             Spliced:
-----------------------------------------------
       Total PE seqs:   1184360
      PE seqs mapped:    539079 (45.5 %)
        Num PE pairs:    592180
   F5 1st end mapped:    372121 (62.8 %)
   F3 2nd end mapped:    166958 (28.2 %)
     PE pairs mapped:     80975 (13.7 %)
     PE proper pairs:     16817 (2.8 %)
-----------------------------------------------

Since this was a paired end alignment there is paired-end specific information reported. Note that this Yeast dataset was of poor quality, especially the R2 reads. You can see this by the relatively low R1 alignment rate (~63%) and the really low R2 alignment rate (~28%).

You can also view statistics on insert sizes for properly paired reads in the bwa_global.iszinfo.txt file. This tells you the average (mean) insert size, standard deviation, mode (most common value), and fivenum values (minimum, 1st quartile, median, 3rd quartile, maximum).

  Insert size stats for: bwa_global
        Number of pairs: 16807 (proper)
 Number of insert sizes: 406
        Mean [-/+ 1 SD]: 296 [176 416]  (sd 120)
         Mode [Fivenum]: 228  [51 224 232 241 500]

A quick way to check alignment stats if you have run multiple alignments is again to use grep. For example:

grep 'Total mapped' *samstats.txt

will produce output like this:

bt2_global.samstats.txt:        Total mapped:    602893 (50.9 %)
bt2_local.samstats.txt:        Total mapped:    788069 (66.5 %)
bwa_global.samstats.txt:        Total mapped:    539079 (45.5 %)
bwa_local.samstats.txt:        Total mapped:   1008000 (76.5 %

Exercise: How would you list the median insert size for all the alignments?

grep 'Mode' *.iszinfo.txt | awk '{print $1,"Median insert size:",$7}'

TACC batch system considerations

The great thing about pipeline scripts like this is that you can perform alignments on many datasets in parallel at TACC, and they are written to take advantage of having multiple cores on TACC nodes where possible.

On the ls6 the pipeline scripts are designed to run best with no more than 4 tasks per node. Although each ls6 node has 128 physical cores per node, the alignment workflow is heavily I/O bound overall, and we don't want to overload the file system.

Exercise #4: Bowtie2 alignment - Vibrio cholerae RNA-seq

While we have focused on aligning eukaryotic data, the same tools can be used with prokaryotic data. The major differences are less about the underlying data and much more about the external/public databases that store and distribute reference data. If we want to study a prokaryote, the reference data is usually downloaded from a resource like GenBank.  

In this exercise, we will use some RNA-seq data from Vibrio cholerae, published on GEO here, and align it to a reference genome.

Overview of Vibrio cholerae alignment workflow with Bowtie2

Alignment of this prokaryotic data follows the workflow below. Here we will concentrate on steps 1 and 2.

1. Prepare the vibCho reference index for bowtie2 from GenBank records
2. Align reads using bowtie2, producing a SAM file
3. Convert the SAM file to a BAM file (samtools view) 
4. Sort the BAM file by genomic location (samtools sort)
5. Index the BAM file (samtools index)
6. Gather simple alignment statistics (samtools flagstat and samtools idxstats)

Obtaining the GenBank records

First prepare a directory for the vibCho fasta, and change to it:

mkdir -p $SCRATCH/core_ngs/references/fasta
cd $SCRATCH/core_ngs/references/fasta

V. cholerae has two chromosomes. We download each separately.

1. Navigate to http://www.ncbi.nlm.nih.gov/nuccore/NC_012582
   • click on the Send to down arrow (top right of page)
     • select Complete Record
     • select File as Destination, and Format FASTA
     • click Create File
       
   • in the Opening File dialog, select Save File then OK
     • Save the file on your local computer as NC_012582.fa
2. Back on the main http://www.ncbi.nlm.nih.gov/nuccore/NC_012582 page
   • click on the Send to down arrow (top right of page)
     • select Complete Record
     • select File as Destination, and Format GFF3
     • click Create File
       
   • in the Opening File dialog, select Save File then OK
     • Save the file on your local computer as NC_012582.gff3
3. Repeat steps 1 and 2 for the 2nd chromosome
   
   • NCBI URL is http://www.ncbi.nlm.nih.gov/nuccore/NC_012583
   • use NC_012583 as the filename prefix for the files you save
   • you should now have 4 files:
     • NC_012582.fa, NC_012582.gff3
     • NC_012583.fa, NC_012583.gff3
       
4. Transfer the files from your local computer to TACC
   • to the ~/scratch/core_ngs/references/vibCho directory created above
     
     • On a Mac or Windows 10 or later, use scp from your laptop
       
     • Otherwise on Windows, use the pscp.exe PuTTy tool
     • See Copying files between TACC and your laptop

mkdir -p $SCRATCH/core_ngs/references/fasta
cd $SCRATCH/core_ngs/references/fasta
cp $CORENGS/references/fasta/NC_* .

Once you have the 4 files locally in your $SCRATCH/core_ngs/references/vibCho directory, combine them using cat:

cd $SCRATCH/core_ngs/references/fasta
cat NC_01258[23].fa   > vibCho.O395.fa
cat NC_01258[23].gff3 > vibCho.O395.gff3

# verify there are 2 contigs in vibCho.O395.fa
grep -P '^>' vibCho.O395.fa

Now we have a reference sequence file that we can use with the bowtie2 reference builder, and ultimately align sequence data against.

Introducing bowtie2

First make sure you're in an idev session:

idev -m 120 -A OTH21164 -N 1 -r CoreNGS      # or -A TRA23004 
# or
idev -m 90 -A OTH21164 -N 1 -p development   # or -A TRA23004 

Go ahead and load the bowtie2 module so we can examine some help pages and options.

module load biocontainers
module load bowtie2

 Now that it's loaded, check out the options. There are a lot of them! In fact for the full range of options and their meaning, Google "Bowtie2 manual" and bring up that page (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml). The Table of Contents is several pages long! Ouch!

This is the key to using bowtie2 - it allows you to control almost everything about its behavior, which make it the go-to aligner for specialized alignment tasks (e.g. aligning miRNA or other small reads). But it also makes it is much more challenging to use than bwa – and it's easier to screw things up too!

Building the bowtie2 vibCho index

Before the alignment, of course, we've got to build a bowtie2- specific index using bowtie2-build. Go ahead and check out its options. Unlike for the aligner itself, we only need to worry about a few things here:

• reference_in file is just the vibCho.O395.fa FASTA we built from GenBank records
• bt2_index_base is the prefix of all the bowtie2-build output file

Here, to build the reference index for alignment, we only need the FASTA file. (This is not always true - extensively spliced transcriptomes requires splice junction annotations to align RNA-seq data properly.)

First create a directory specifically for the bowtie2 index, then build the index using bowtie-build.

mkdir -p $SCRATCH/core_ngs/references/fasta
cd $SCRATCH/core_ngs/references/fasta
cp $CORENGS/references/fasta/vibCho* .

mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cd $SCRATCH/core_ngs/references/bt2/vibCho

# Symlink to the fasta file you created using relative path syntax
ln -sf ../../fasta/vibCho.O395.fa

bowtie2-build vibCho.O395.fa vibCho.O395

This should also go pretty fast. You can see the resulting files using ls like before.

Performing the bowtie2 alignment

Make sure you're in an idev session with the bowtie2 BioContainers module loaded:

idev -m 120 -A OTH21164  -N 1 -r CoreNGS    # or -A TRA23004
# or
idev -m 90 -A OTH21164 -N 1 -p development  # or -A TRA23004

module load biocontainers
module load bowtie2

We'll set up a new directory to perform the V. cholerae data alignment. But first make sure you have the FASTQ file to align and the vibCho bowtie2 index.

# Get a pre-built vibCho index if you didn't already build one
mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/references/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho/

# Get the FASTQ to align
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz  $SCRATCH/core_ngs/alignment/fastq/

Now set up a directory to do this alignment, with symbolic links to the bowtie2 index directory and the directory containing the FASTQ to align:

mkdir -p $SCRATCH/core_ngs/alignment/vibCho
cd $SCRATCH/core_ngs/alignment/vibCho
ln -sf ../../references/bt2/vibCho
ln -sf ../../alignment/fastq fq

We'll be aligning the V. cholerae reads now in ./fq/cholera_rnaseq.fastq.gz (how many sequences does it contain?)

Note that here the data is from standard mRNA sequencing, meaning that the DNA fragments are typically longer than the reads. There is likely to be very little contamination that would require using a local rather than global alignment, or many other pre-processing steps (e.g. adapter trimming). Thus, we will run bowtie2 with default parameters, omitting options other than the input, output, and reference index. This performs a global alignment.

As you can tell from looking at the bowtie2 help message, the general syntax looks like this:

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

So execute this bowtie2 global, single-end alignment command:

mkdir -p $SCRATCH/core_ngs/references/fasta
cp $CORENGS/references/fasta/vibCho*  $SCRATCH/core_ngs/references/fasta/

mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/references/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho/

mkdir -p $SCRATCH/core_ngs/alignment/vibCho
cd $SCRATCH/core_ngs/alignment/vibCho
ln -sf ../../references/bt2/vibCho
ln -sf ../../alignment/fastq fqmkdir -p $SCRATCH/core_ngs/alignment/vibCho

cd $SCRATCH/core_ngs/alignment/vibCho
bowtie2 -x vibCho/vibCho.O395 -U fq/cholera_rnaseq.fastq.gz \
  -S cholera_rnaseq.sam 2>&1 | tee aln_global.log

Notes:

• -x  vibCho/vibCho.O395.fa – prefix path of index files
• -U fq/cholera_rnaseq.fastq.gz – FASTQ file for single-end (Unpaired) alignment
• -S cholera_rnaseq.sam – tells bowtie2 to report alignments in SAM format to the specified file
• 2>&1 redirects standard error to standard output
  • while the alignment data is being written to the cholera_rnaseq.sam file, bowtie2 will report its progress to standard error.
• | tee aln.log takes the bowtie2 progress output and pipes it to the tee program
  • tee takes its standard input and writes it to the specified file and also to standard output
  • that way, you can see the progress output now, but also save it to review later (or supply to MultiQC)

Since the FASTQ file is not large, this should not take too long, and you will see progress output like this:

89006 reads; of these:
  89006 (100.00%) were unpaired; of these:
    5902 (6.63%) aligned 0 times
    51483 (57.84%) aligned exactly 1 time
    31621 (35.53%) aligned >1 times
93.37% overall alignment rate

When the job is complete you should have a cholera_rnaseq.sam file that you can examine using whatever commands you like.  Remember, to further process it downstream, you should create a sorted, indexed BAM file from this SAM output.

Exercise: Repeat the alignment performing a local alignment, and write the output in BAM format. How do the alignment statistics compare?

module load samtools
cd $SCRATCH/core_ngs/alignment/vibCho
bowtie2 --local -x vibCho/vibCho.O395 -U fq/cholera_rnaseq.fastq.gz 2>aln_local.log \
 | samtools view -b > cholera_rnaseq.local.bam

89006 reads; of these:
  89006 (100.00%) were unpaired; of these:
    13359 (15.01%) aligned 0 times
    46173 (51.88%) aligned exactly 1 time
    29474 (33.11%) aligned >1 times
84.99% overall alignment rate

Exercise #5: BWA-MEM - Human mRNA-seq

After bowtie2 came out with a local alignment option, it wasn't long before bwa developed its own local alignment algorithm called BWA-MEM (for Maximal Exact Matches), implemented by the bwa mem command. 

bwa mem has the following advantages:

• It provides the simplicity of using bwa without the complexities of local alignment
• It can align different portions of a read to different locations on the genome
  • In a total RNA-seq experiment, reads will (at some frequency) span a splice junction themselves
    • or a pair of reads in a paired-end library will fall on either side of a splice junction.
  • We want to be able to align these splice-adjacent reads for many reasons, from accurate transcript quantification to novel fusion transcript discovery.

This exercise will align a human total RNA-seq dataset that includes numerous reads that cross splice junctions.

A word about real splice-aware aligners

Using bwa mem for RNA-seq alignment is sort of a "poor man's" RNA-seq alignment method. Real splice-aware aligners like tophat2 STAR, hisat2 or star tophat have more complex algorithms (as shown below) – and take a lot more time!

Image RemovedImage Added

...

In the transcriptome-

...

Now, try aligning it with bwa aln like we did in Example #1, but first link to the hg19 bwa index directory.  In this case, due to the size of the hg19 index, we are linking to Anna's scratch area INSTEAD of our own work area containing indexes that we built ourselves.

cd $SCRATCH/core_ngs/alignment
ln -s -f /scratch/01063/abattenh/ref_genome/bwa/bwtsw/hg19
ls hg19

You should see a set of files analogous to the yeast files we created earlier, except that their universal prefix is hg19.fa.  

Go ahead and try to do a single-end alignment of the file to the human genome using bwa aln like we did in Exercise #1, saving intermediate files with the prefix human_rnaseq_bwa. Go ahead and just execute on the command lineaware alignment above, reads that span splice junctions are reported in the SAM file with genomic coordinates that start in the first exon and end in the second exon (the CIGAR string uses the N operator, e.g. 30M1000N60M).

BWA MEM does not know about the exon structure of the genome. But it can align different sub-sections of a read to two different locations, producing two alignment records from one input read. One of the two will be marked as secondary (0x100 flag).

Setup for BWA mem

First set up our working directory for this alignment. Since it takes a long time to build a bwa index for a large genome (here human hg38/GRCh38), we'll use one that the BioITeam maintains in its /work/projects/BioITeam/ref_genome area.

bwa aln hg19/hg19.fa fastq/human_rnaseq.fastq.gz > human_rnaseq_bwa.sai
bwa samse hg19/hg19.fa human_rnaseq_bwa.sai fastq/human_rnaseq.fastq.gz > human_rnaseq_bwa.sam

Once this is complete use less to take a look at the contents of the SAM file, using the space bar to leaf through them. You'll notice a lot of alignments look basically like this:

HWI-ST1097:228:C21WMACXX:8:1316:10989:88190     4       *       0       0       *       *       0       0
  AAATTGCTTCCTGTCCTCATCCTTCCTGTCAGCCATCTTCCTTCGTTTGATCTCAGGGAAGTTCAGGTCTTCCAGCCGCTCTTTGCCACTGATCTCCAGCT
  CCCFFFFFHHHHHIJJJJIJJJJIJJJJHJJJJJJJJJJJJJJIIIJJJIGHHIJIJIJIJHBHIJJIIHIEGHIIHGFFDDEEEDDCDDD@CDEDDDCDD

Notice that the contig name (field 3) is just an asterisk ( * ) and the alignment flags value is a 4 (field 2), meaning the read did not align (decimal 4 = hex 0x4 = read did not map).

Essentially, nothing (with a few exceptions) aligned. Why?

RNA-seq alignment with bwa mem

Exercise: use bwa mem to align the same data

Based on the following syntax and the above reference path, use bwa mem to align the same file, saving output files with the prefix human_rnaseq_mem. Go ahead and just execute on the command line.

# Make sure you're in an idev session
idev -m 120 -N 1 -A OTH21164 -r CoreNGS      # or -A TRA23004
# or
idev -m 90 -N 1 -A OTH21164 -p development   # or -A TRA23004

# Load the modules we'll need
module load biocontainers
module load bwa
module load samtools

# Copy over the FASTQ data if needed
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*.gz $SCRATCH/core_ngs/alignment/fastq/

# Make a new alignment directory for running these scripts
cds
mkdir -p core_ngs/alignment/bwamem
cd core_ngs/alignment/bwamem
ln -sf ../fastq
ln -sf /work/projects/BioITeam/ref_genome/bwa/bwtsw/hg38

Now take a look at bwa mem usage (type bwa mem with no arguments, or bwa mem 2>&1 | more). The most important parameters are the following:

RNA-seq alignment with bwa mem

Based on its help info, this is the structure of the bwa mem command we will use:

bwa mem -M <ref.fa> <reads.fq> > outfile.sam

After performing the setup above, execute the following command in your idev session:

cd $SCRATCH/core_ngs/alignment/bwamem
bwa mem -M hg38/hg38.fa fastq/human_rnaseq.fastq.gz 

Check the length of the SAM file you generated with wc -l. Since there is one alignment per line, there must be 586266 alignments (minus no more than 100 header lines), which is more than the number of sequences in the FASTQ file. This is bwa mem can report multiple alignment records for the same read, hopefully on either side of a splice junction. These alignments can still be tied together because they have the same read ID.

# This gives you a view where each read is listed next to the number of entries it has in the SAM file
cut -f 1 human_rnaseq_mem.sam | sort | uniq -c | less

#This gives essentially a histogram of the number of times each read aligned - a plurality of reads aligned twice, which seems reasonable since these are all reads crossing a junction, but plenty aligned more or less
cut -f 1 human_rnaseq_mem.sam | sort | uniq -c | awk '{print $1}' | sort | uniq -c | less	

#This gives a better idea of the alignment rate, which is how many reads aligned at least once to the genome.  Divided by the number of reads in the original file, the real alignment rate is much higher.
cut -f 1 human_rnaseq_mem.sam | sort | uniq | wc -l	

# NOTE: Some of these one-liners are only reasonably fast if the files are relatively small (around a million reads or less). For bigger files, there are better ways to get this information, mostly using samtools.

BWA-MEM vs Tophat

Another approach to aligning long RNA-seq data is to use an aligner that is more explicitly concerned with sensitivity to splice sites, namely a program like Tophat. Tophat uses either bowtie (tophat) or bowtie2 (tophat2) as the actual aligner, but performs the following steps:

• aligns reads to the genome
• reads that do not align to the genome are aligned against a transcriptome, if provided
  • if they align, the transcriptome coordinates are converted back to genomic coordinates, with gaps represented in the CIGAR string, for example as 196N
• reads that do not align to the transcriptome are split into smaller pieces, each of which Tophat attempts to map to the genome

Note that Tophat also reports secondary alignments, but they have a different meaning. Tophat always reports spliced alignments as one alignment records with the N CIGAR string operator indicating the gaps. Secondary alignments for Tophat (marked with the 0x100 BAM flag) represent alternate places in the genome where a  read (spliced or not) may have mapped.

...

2>hs_rna.bwamem.log | \
  samtools view -b | \
  samtools sort -O BAM -T human_rnaseq.tmp > human_rnaseq.sort.bam

This multi-pipe command performs three steps:

• The bwa mem alignment
  • the program's progress output (on standard error) is redirected to a log file (2>hs_rna.bwamem.log)
  • its alignment records (on standard output) is piped to the next step (conversion to BAM)
• Conversion of bwa mem's SAM output to BAM format
  • recall that the -b option to samtools view says to output in BAM format
• Sorting the BAM file
  • samtools sort takes the binary output from samtools view and writes a sorted BAM file.

Because the progress output is being redirected to a log file, we won't see anything until the command completes. Then you should have a human_rnaseq.sort.bam file and an hs_rna.bwamem.log logfile.

Exercise: Compare the number of original FASTQ reads to the number of alignment records.

cd $SCRATCH/core_ngs/alignment/bwamem
zcat ./fastq/human_rnaseq.fastq.gz | wc -l | awk '{print $1/4}'

cd $SCRATCH/core_ngs/alignment/bwamem
samtools flagstat human_rnaseq.sort.bam | tee hs_rnaseq.flagstat.txt

133570 + 0 in total (QC-passed reads + QC-failed reads)
33570 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
133450 + 0 mapped (99.91% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)