More Alignment exercises

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.

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.

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

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 prebuilt 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 similar arguments to the BWA alignment script.

    • again, we specify that reads should first be trimmed to 50 bp.

  • all output files associated with this command will be named with the prefix bt2_global.

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

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

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:

This will show something like:

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

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:

If this command returns 4 (the number of alignment tasks we performed), all went well, and we're done.  Another way is to list all the filenames that contain that pattern using grep -l:

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 the special grep -L option to do this work.

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 files that don't contain that message:

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:

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%).

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

will produce output like this:

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

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

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

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.

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

1. Obtain reference FASTA for vibCho from NCBI GenBank records

2. Prepare the vibCho reference index for bowtie2 from the GenBank FASTA

3. Align reads using bowtie2, producing a SAM file

4. Convert the SAM file to a BAM file (samtools view) 

5. Sort the BAM file by genomic location (samtools sort)

6. Index the BAM file (samtools index)

7. Gather simple alignment statistics (samtools flagstat and samtools idxstats)

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

To look for a Genome at NCBI, navigate to: https://www.ncbi.nlm.nih.gov/datasets/genome/

1. Enter Vibrio cholerae as the Search term

   1. The first assembly listed is ASM836960v1
      
      

      1. It has a checkmark beside it, meaning this is the "preferred" Genome

2. Click on the ASM836960v1 link

   1. Takes you to the Genome assembly ASM836960v1 page

3. Click on the FTP link

   1. Right click on the GCF_008369605.1_ASM836960v1_genomic.fna.gz link and Copy link address 

4. In your $SCRATCH/core_ngs/references/fasta directory, type wget then Paste

   1. you should see:
      wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/008/369/605/GCF_008369605.1_ASM836960v1/GCF_008369605.1_ASM836960v1_genomic.fna.gz 

   2. Press Enter

   3. You should see the file GCF_008369605.1_ASM836960v1_genomic.fna.gz.

5. Go back to the FTP page

   1. Right click on the GCF_008369605.1_ASM836960v1_genomic.gtf.gz link and Copy link address

6. In your $SCRATCH/core_ngs/references/fasta directory, type wget then Paste

   1. you should see:
      wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/008/369/605/GCF_008369605.1_ASM836960v1/GCF_008369605.1_ASM836960v1_genomic.gtf.gz
      
      

   2. Press Enter

   3. You should see the file GCF_008369605.1_ASM836960v1_genomic.gtf.gz.

7. Unzip and rename the two GZIP files:
   cd $SCRATCH/core/ngs/references/fasta
gunzip GCF_008369605.1_ASM836960v1_genomic.fna.gz > vibCho_ASM836960v1.fna
gunzip GCF_008369605.1_ASM836960v1_genomic.gtf.gz> vibCho_ASM836960v1.gtf
   
   
   

Exercise: How many contigs are in the Vibrio cholerae FASTQ reference genome?

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

First make sure you're in an idev session:

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

 Now that it's loaded, check out the options (bowtie2 --help | more). 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!

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_ASM836960v1.fna FASTA from GenBank records

• bt2_index_base is the prefix of all the bowtie2-build output file

To build the reference index for alignment we only need the FASTA file, and a prefix that we want all the output files to have.

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

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

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

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.

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:

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

Exercise: How many sequences does the ./fq/cholera_rnaseq.fastq.gz file 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 preprocessing 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:

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

Notes: