Tip | ||
---|---|---|
| ||
Use our summer school reservation (CoreNGS-Thu) when submitting batch jobs to get higher priority on the ls6 normal queue today:
|
Table of Contents |
---|
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, /work2work/projects/BioITeam/common/script/align_bwa_illumina.sh. Type in the script name to see its usage.
Code Block | ||
---|---|---|
| ||
align_bwa_illumina.sh 20212022_0605_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 |
...
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.
Code Block | ||||
---|---|---|---|---|
| ||||
# Make sure you're not in an idev session by looking at the hostname hostname # If the hostname looks like "c455-004.stampede2ls6.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 -s -f ../fastq # Copy the alignment commands file and submit the batch job cd $SCRATCH/core_ngs/alignment/bwa_script cp $CORENGS/tacc/aln_script.cmds . launcher_creator.py -j aln_script.cmds -n aln_script -t 0201:00:00 -w 4 -a UT-2015-05-18OTH21164 -q normal sbatch --reservation=BIO_DATA_week_1CoreNGS-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.
Code Block | ||||
---|---|---|---|---|
| ||||
/work2work/projects/BioITeam/common/script/align_bwa_illumina.sh global ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bwa_global sacCer3 1 50 /work2work/projects/BioITeam/common/script/align_bwa_illumina.sh local ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bwa_local sacCer3 1 /work2work/projects/BioITeam/common/script/align_bowtie2_illumina.sh global ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bt2_global sacCer3 1 50 /work2work/projects/BioITeam/common/script/align_bowtie2_illumina.sh local ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz bt2_local sacCer3 1 |
...
- 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--global with 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
- BWA global: bwa aln the R1 and R2 separately, then bwa sampe 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
...
This will show something like:
Code Block |
---|
..Done ------------------------------------------------------------------ ..Done alignmentUtils.pl bamstats - 20202022-06-1410 2312:1959:3805 .. samstats file 'bwa_global.samstats.txt' exists and is not empty - 20202022-06-1410 2312:1959:3805 =============================================================================== ## Cleaning up files (keep 0) - 20202022-06-1410 2312:1959:3805 =============================================================================== ckRes 0 cleanup =============================================================================== ## All bwa alignment tasks completed successfully! - 20202022-06-1410 2312:1959:3806 =============================================================================== |
...
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 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.
Tip | ||
---|---|---|
| ||
These alignment scripts should always be run with a wayness of 4 (-w 4) in the stampede2 ls6 batch system, meaning at most 4 commands per node. |
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.
...
- Prepare the vibCho reference index for bowtie2 from GenBank records
- Align reads using bowtie2, producing a SAM file
- Convert the SAM file to a BAM file (samtools view)
- Sort the BAM file by genomic location (samtools sort)
- Index the BAM file (samtools index)
- Gather simple alignment statistics (samtools flagstat and samtools idxstatidxstats)
Obtaining the GenBank records
First prepare a directory to work infor the vibCho fasta, and change to it:
Code Block | ||
---|---|---|
| ||
mkdir -p $SCRATCH/core_ngs/references/vibChofasta cd $SCRATCH/core_ngs/references/vibChofasta |
V. cholerae has two chromosomes. We download each separately.
- 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
- click on the Send to down arrow (top right of page)
- 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
- click on the Send to down arrow (top right of page)
- 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
- 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 shell10 or later, use scp from your laptop
- On Otherwise on Windows, use the pscp.exe PuTTy tool
- See Copying files between TACC and your laptop
- On a Mac or Windows with a WSL shell10 or later, use scp from your laptop
- to the ~/scratch/core_ngs/references/vibCho directory created above
...
Expand | |||||
---|---|---|---|---|---|
| |||||
|
Once you have the 4 files locally in your $SCRATCH/core_ngs/references/vibCho directory, combine them using cat:
Code Block | ||
---|---|---|
| ||
cd $SCRATCH/core_ngs/references/vibChofasta 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.
...
Code Block | ||||
---|---|---|---|---|
| ||||
idev -m 180120 -A OTH21164 -pN normal1 -Ar UT-2015-05-18CoreNGS-Thu # or idev -m 90 -A OTH21164 -N 1 -np 68development |
Go ahead and load the bowtie2 module so we can examine some help pages and options.
Code Block | ||
---|---|---|
| ||
module load biocontainers module load bowtiebowtie2 |
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!
...
First create a directory specifically for the bowtie2 index, then build the index using bowtie-build.
Code Blockexpand | |||||
---|---|---|---|---|---|
| |||||
|
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
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 $SCRATCH/core_ngs/references../../fasta/vibCho.O395.fa # or, to catch up: ln -sf $CORENGS/references/vibCho bowtie2-build 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 indexMake sure you're in an idev session with the bowtie2 BioContainers module loaded:
Code Block | ||
---|---|---|
| ||
#idev Get-m the120 FASTQ-A toOTH21164 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/vibChoN 1 -r CoreNGS-Thu # or idev -m 90 -A OTH21164 -N 1 -p development 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.Make sure you're in an idev session with the bowtie2 BioContainers module loaded:
Code Block | ||
---|---|---|
| ||
idev# -pGet normala pre-mbuilt 120vibCho -A UT-2015-05-18 -N 1 -n 68 module load biocontainers module load bowtieindex 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 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:
...
We'll be aligning the V. cholerae reads now in ./fq/cholera_rnaseq.fastq.gz (how many sequences does it contain?)
...
So execute this bowtie2 global, single-end alignment command:
Expand | |||||
---|---|---|---|---|---|
| |||||
|
Notes:
|
Code Block | ||
---|---|---|
| ||
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
- -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:
Code Block |
---|
89006 reads; of these: 89006 (100.00%) were unpaired; of these: 206755902 (236.23%63%) aligned 0 times 3822651483 (4257.95%84%) aligned exactly 1 time 3010531621 (3335.82%53%) aligned >1 times 7693.77%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.
...
Expand | |||||||
---|---|---|---|---|---|---|---|
| |||||||
Reports these alignment statistics:
Interestingly, the local alignment rate here is lower than we saw with the gloabl alignmentglobal alignment. Usually local alignments have higher alignment rates than corresponding global ones. |
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.
- In a total RNA-seq experiment, reads will (at some frequency) span a splice junction themselves
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, hisat2 or STAR have more complex algorithms (as shown below) – and take a lot more time!
In the transcriptome-aware 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).
BWA MEM splits junction-spanning reads into two alignment records |
---|
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.
Code Block | ||||
---|---|---|---|---|
| ||||
# Make sure you're in an idev session
idev -m 120 -N 1 -A OTH21164 -r CoreNGS-Thu
# or
idev -m 90 -N 1 -A OTH21164 -p development
# 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). The most important parameters are the following:
Option | Effect |
---|---|
-k | Controls the minimum seed length (default = 19) |
-w | Controls the "gap bandwidth", or the length of a maximum gap. This is particularly relevant for MEM, since it can determine whether a read is split into two separate alignments or is reported as one long alignment with a long gap in the middle (default = 100) |
-M | For split reads, mark the shorter read as secondary |
-r | Controls how long an alignment must be relative to its seed before it is re-seeded to try to find a best-fit local match (default = 1.5, e.g. the value of -k multiplied by 1.5) |
-c | Controls how many matches a MEM must have in the genome before it is discarded (default = 10000) |
-t | Controls the number of threads to use |
RNA-seq alignment with bwa mem
Based on its help info, this is the structure of the bwa mem command we will use:
Code Block |
---|
bwa mem -M <ref.fa> <reads.fq> > outfile.sam |
After performing the setup above, execute the following command in your idev session:
Code Block | ||
---|---|---|
| ||
cd $SCRATCH/core_ngs/alignment/bwamem
bwa mem -M hg38/hg38.fa fastq/human_rnaseq.fastq.gz 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.
Expand | ||
---|---|---|
| ||
Use the zcat | wc -l | awk idiom to count FASTQ reads. Use samtools flagstat to report alignment statistics. |
Expand | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
| ||||||||||||
Count the FASTQ file reads:
The file has 100,000 reads. Generate alignment statistics from the sorted BAM file:
Output will look like this:
There were 133,570 alignment records reported for the 100,000 input reads. Because bwa mem can split reads and report two alignment records for the same read, there are 33,570 secondary reads reported here. |
Tip |
---|
Be aware that some downstream tools (for example the Picard suite, often used before SNP calling) do not like it when a read name appears more than once in the SAM file. Such reads can be filtered, but only if they can be identified as secondary by specifying the bwa mem -M option as we did above. This option reports the longest alignment normally but marks additional alignments for the read as secondary (the 0x100 BAM flag). This designation also allows you to easily filter out the secondary reads with samtools view -F 0x104 if desired. |