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, /work2/projects/BioITeam/common/script/align_bwa_illumina.sh. Type in the script name to see its usage.
align_bwa_illumina.sh 2021_06_05 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:
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.stampede2.tacc.utexas.edu", exit the idev session # 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 -s -f ../fastq # Copy the alignment commands file and submit the batch job cp $CORENGS/tacc/aln_script.cmds . launcher_creator.py -j aln_script.cmds -n aln_script -t 02:00:00 -w 4 -a UT-2015-05-18 -q normal sbatch --reservation=BIO_DATA_week_1 aln_script.slurm showq -u |
While we're waiting for the job to complete, lets look at the aln_script.cmds file.
/work2/projects/BioITeam/common/script/align_bwa_illumina.sh global ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bwa_global sacCer3 1 50 /work2/projects/BioITeam/common/script/align_bwa_illumina.sh local ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bwa_local sacCer3 1 /work2/projects/BioITeam/common/script/align_bowtie2_illumina.sh global ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bt2_global sacCer3 1 50 /work2/projects/BioITeam/common/script/align_bowtie2_illumina.sh local ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bt2_local sacCer3 1 |
Notes:
This alignment pipeline script performs the following steps:
There are a number of output files, with the most important being those desribed below.
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 - 2020-06-14 23:19:38 .. samstats file 'bwa_global.samstats.txt' exists and is not empty - 2020-06-14 23:19:38 =============================================================================== ## Cleaning up files (keep 0) - 2020-06-14 23:19:38 =============================================================================== ckRes 0 cleanup =============================================================================== ## All bwa alignment tasks completed successfully! - 2020-06-14 23:19:38 =============================================================================== |
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 |
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.
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?
That information is in the *.iszinfo.txt files, on the line labeled Mode. The median value is th 3rd value in the 5 fivnum values; it is the 7th whitespace-separated field on the Mode line. |
Use grep to isolate the Mode line, and awk to isolate the median value field:
|
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 stampede2, with its 68 physical cores per node, they are designed to run best with no more than 4 tasks per node.
These alignment scripts should always be run with a wayness of 4 (-w 4) in the stampede2 batch system, meaning at most 4 commands per node. |
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.
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.
|
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.
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!
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:
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.
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_global.log |
Notes:
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: Repeat the alignment performing a local alignment, and write the output in BAM format. How do the alignment statistics compare?
--local |
Reports these alignment statistics:
Interestingly, the local alignment rate here is lower than we saw with the gloabl alignment. |