Versions Compared

Old Version 27

changes.mady.by.user Anna Battenhouse

Saved on

compared with

New Version Current

changes.mady.by.user Daryl Barth

Saved on

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.
Comment: Added end of line carriate to bwa mem call


Tip
titleReservations

Use our summer school reservation (

CoreNGSday4

CoreNGS) when submitting batch jobs to get higher priority on the ls6 normal queue

today:sbatch --reservation=CoreNGSday4

.

Code Block
languagebash
titleRequest an interactive (idev) node
# 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


Code Block
languagebash
titleSubmit a batch job
# Using our reservation
sbatch --reseservation=CoreNGS <batch_file>.slurm

Note that the reservation name (CoreNGS) is different from the TACC allocation/project for this class, which is OTH21164.

Table of Contents

Exercise #3: PE alignment with BioITeam scripts

...

Code Block
languagebash
align_bwa_illumina.sh 2022_0506_0510
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

...

Code Block
languagebash
titleRun multiple alignments using the TACC batch system
# 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 -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 .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 normaldevelopment
sbatch --reservation=CoreNGSday4 aln_script.slurm  

showq -u

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

...

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

...

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

...

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

Obtaining the GenBank records

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

Code Block
languagebash
mkdir -p $SCRATCH/core_ngs/references/fasta
cd $SCRATCH/core_ngs/references/fasta

...

First make sure you're in an idevsession:

...

languagebash
titleStart an idev session

session:

Code Block
languagebash
titleStart an idev session
idev -m 120 -A OTH21164 -N 1 -r CoreNGS      # or -A TRA23004 
# or
idev -m 12090 -A OTH21164 -N 1 -p development   # or -rA CoreNGSday4TRA23004

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

...

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

Code Block
languagebash
idev -m 120 -A OTH21164  -N 1 -r CoreNGS    # or -A TRA23004
# or
idev -m 120 90 -A OTH21164 -N 1 -Ap OTH21164development  -N# 1or -rA TRA23004
CoreNGSday4
module load biocontainers
module load bowtie2

...

Code Block
languagebash
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/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:

...

Expand
titleAnswer


Code Block
languagebash
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

Reports these alignment statistics:

Code Block
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

Interestingly, the local alignment rate here is lower than we saw with the global alignment. Usually local alignments have higher alignment rates than corresponding global ones.

...

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.

...

Using bwa mem for RNA-seq alignment is sort of a "poor man's" RNA-seq alignment method. Real splice-aware aligners like tophat2STAR, hisat2 or STAR tophat 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 . One of the two will be marked as secondary (0x100 flag).

...

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
languagebash
titleRun multiple alignments using the TACC batch system
# Make sure you're in an idev session
idev -m 120 -N 1 -A OTH21164 -r CoreNGSday4  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:

...

Code Block
languagebash
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

...

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

...

Expand
titleAnswer:

Count the FASTQ file reads:

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

The file has 100,000 reads.

Generate alignment statistics from the sorted BAM file:

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

Output will look like this:

Code Block
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)

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.

...