Page Comparison

...

As you can imagine from this series of steps, Tophat is very computationally intensive and takes much longer than bwa mem – very large alignments (hundreds of millions of reads) may not complete in stampede's 48 hour maximum job time!

Exercise #5: Simple SAMtools Utilities

We have used several alignment methods that all generate results in the form of the near-universal SAM/BAM file format.  The SAMtools program is a ubiquitously used set of tools that allow a user to manipulate SAM/BAM files in many different ways, ranging from simple tasks (like SAM/BAM interconversion) to more complex functions (like removal of PCR duplicates).  It is available in the TACC module system in the typical fashion.

In this exercise, we will use five very simple utilities provided by samtools: view, sort, index, flagstat, and idxstats. Each of these is executed in one line for a given SAM/BAM file. In the SAMtools/BEDtools section tomorrow we will explore samtools in more in depth.

For the sake of time and simplicity, here we are only going to run these commands on the yeast paired-end alignment file. The same commands can be run on the other files by changing the names, so feel free to try them on other SAM files. Indeed, it is very common in practice to use bash loops to generate many commands for a large set of alignments and deposit those commands into a batch job cmds file for submission.

To start, we will move to the directory containing our SAM files, among other things, and load up samtools using the module system. After loading it, just run the samtools command to see what the available tools are (and to see what the syntax of an actual command is).

cd $SCRATCH/core_ngs/alignment
ls -la
module load samtools
samtools

You will see the following screen after running samtools with no other options:

Program: samtools (Tools for alignments in the SAM format)
Version: 1.2 (using htslib 1.2.1)
Usage:   samtools <command> [options]
Commands:
  -- indexing
         faidx       index/extract FASTA
         index       index alignment
  -- editing
         calmd       recalculate MD/NM tags and '=' bases
         fixmate     fix mate information
         reheader    replace BAM header
         rmdup       remove PCR duplicates
         targetcut   cut fosmid regions (for fosmid pool only)
  -- file operations
         bamshuf     shuffle and group alignments by name
         cat         concatenate BAMs
         merge       merge sorted alignments
         mpileup     multi-way pileup
         sort        sort alignment file
         split       splits a file by read group
         bam2fq      converts a BAM to a FASTQ
  -- stats
         bedcov      read depth per BED region
         depth       compute the depth
         flagstat    simple stats
         idxstats    BAM index stats
         phase       phase heterozygotes
         stats       generate stats (former bamcheck)
  -- viewing
         flags       explain BAM flags
         tview       text alignment viewer
         view        SAM<->BAM<->CRAM conversion

Samtools view

The utility samtools view provides a way of converting SAM (text format) files to BAM (binary, compressed) files directly. It also provides many, many other functions which we will discuss lster. To get a preview, execute samtools view without any other arguments. You should see:

Usage:   samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]
Options: -b       output BAM
         -C       output CRAM (requires -T)
         -1       use fast BAM compression (implies -b)
         -u       uncompressed BAM output (implies -b)
         -h       include header in SAM output
         -H       print SAM header only (no alignments)
         -c       print only the count of matching records
         -o FILE  output file name [stdout]
         -U FILE  output reads not selected by filters to FILE [null]
         -t FILE  FILE listing reference names and lengths (see long help) [null]
         -T FILE  reference sequence FASTA FILE [null]
         -L FILE  only include reads overlapping this BED FILE [null]
         -r STR   only include reads in read group STR [null]
         -R FILE  only include reads with read group listed in FILE [null]
         -q INT   only include reads with mapping quality >= INT [0]
         -l STR   only include reads in library STR [null]
         -m INT   only include reads with number of CIGAR operations
                  consuming query sequence >= INT [0]
         -f INT   only include reads with all bits set in INT set in FLAG [0]
         -F INT   only include reads with none of the bits set in INT
                  set in FLAG [0]
         -x STR   read tag to strip (repeatable) [null]
         -B       collapse the backward CIGAR operation
         -s FLOAT integer part sets seed of random number generator [0];
                  rest sets fraction of templates to subsample [no subsampling]
         -@ INT   number of BAM compression threads [0]
         -?       print long help, including note about region specification
         -S       ignored (input format is auto-detected)

That is a lot to process! For now, we just want to read in a SAM file and output a BAM file. The input format is auto-detected, so we don't need to say that we're inputing a SAM instead of a BAM. We just need to tell the tool to output the file in BAM format, and provide the name of the destination BAM file. This command is as follows:

samtools view -b yeast_pairedend.sam -o yeast_pairedend.bam

• the -b option tells the tool to output BAM format
• the -o option specifies the name of the output BAM file that will be created

How do you look at the BAM file contents now? That's simple. Just use samtools view. without the -b option. Remember to pipe output to a pager!

samtools view yeast_pairedend.bam | more

Samtools sort

Look at the SAM file briefly using less. You will notice, if you scroll down, that the alignments are in no particular order, with chromosomes and start positions all mixed up. This makes searching through the file very inefficient. samtools sort is a piece of samtools that provides the ability to re-order entries in the SAM file either by coordinate position or by read name.

If you execute samtools sort without any options, you see its help page:

Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read name
  -o FILE    Write final output to FILE rather than standard output
  -O FORMAT  Write output as FORMAT ('sam'/'bam'/'cram')   (either -O or
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam       -T is required)
  -@ INT     Set number of sorting and compression threads [1]
Legacy usage: samtools sort [options...] <in.bam> <out.prefix>
Options:
  -f         Use <out.prefix> as full final filename rather than prefix
  -o         Write final output to stdout rather than <out.prefix>.bam
  -l,m,n,@   Similar to corresponding options above

In most cases you will be sorting a BAM file by position rather than by name. You can use either -o or reidrection with > to control the output.

To sort the paired-end yeast BAM file by coordinate, and get a BAM file named yeast_pairedend.sort.bam as output, execute the following command:

samtools sort -O bam -T yeast_pairedend.sort yeast_pairedend.bam > yeast_pairedend.sort.bam

• The -O options says the output format should be BAM
• The -T options gives a prefix for temporary files produced during sorting
• By default sort writes its output to standard output, so we use > to redirect to a file named yeast_pairedend.sort.bam

Samtools index

Many tools (like the UCSC Genome Browser) only need to use sub-sections of the BAM file at a given point in time. For example, if you are viewing alignments that are within a particular gene alignments on other chromosomes generally do not need to be loaded. In order to speed up access, BAM files are indexed, producing BAI files which allow other programs to navigate directly to the alignments of interest. This is especially important when you have many alignments.

The utility samtools index creates an index that has the exact name as the input BAM file, with suffix .bai appended. The help page, if you execute samtools index with no arguments, is as follows:

Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]
Options:
  -b       Generate BAI-format index for BAM files [default]
  -c       Generate CSI-format index for BAM files
  -m INT   Set minimum interval size for CSI indices to 2^INT [14]

Here, the syntax is way, way easier. We want a BAI-format index which is the default. (CSI-format is used with extremely long contigs, which we aren't considering here - the most common use case are highly polyploid plant genomes).

So all we have to type is:

samtools index yeast_pairedend.sort.bam

This will produce a file named yeast_pairedend.bam.bai.

Most of the time when an index is required, it will be automatically located as long as it is in the same directory as its BAM file and shares the same name up until the .bai extension.

Samtools idxstats

Now that we have a sorted, indexed BAM file, we might like to get some simple statistics about the alignment run. For example, we might like to know how many reads aligned to each chromosome/contig. The samtools idxstats is a very simple tool that provides this information. If you type the command without any arguments, you will see that it could not be simpler - just type the following command:

samtools idxstats yeast_pairedend.sort.bam

The output is a text file with four tab-delimited columns with the following meanings:

1. chromosome name
2. chromosome length
3. number of mapped reads
4. number of unmapped reads

The reason that the "unmapped reads" field for the named chromosomes is not zero is that, if one half of a pair of reads aligns while the other half does not, the unmapped read is still assigned to the chromosome its mate mapped to, but is flagged as unmapped.

Samtools flagstat

Finally, we might like to obtain some other statistics, such as the percent of all reads that aligned to the genome. The samtools flagstat tool provides very simple analysis of the SAM flag fields, which includes information like whether reads are properly paired, aligned or not, and a few other things. Its syntax is identical to that of samtools idxstats:

samtools flagstat yeast_pairedend.sort.bam

You should see something like this:

1184360 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
547664 + 0 mapped (46.24%:-nan%)
1184360 + 0 paired in sequencing
592180 + 0 read1
592180 + 0 read2
473114 + 0 properly paired (39.95%:-nan%)
482360 + 0 with itself and mate mapped
65304 + 0 singletons (5.51%:-nan%)
534 + 0 with mate mapped to a different chr
227 + 0 with mate mapped to a different chr (mapQ>=5)

Ignore the "+ 0" addition to each line - that is a carry-over convention for counting QA-failed reads that is no longer necessary.

The most important statistic is the mapping rate, but this readout allows you to verify that some common expectations (e.g. that about the same number of R1 and R2 reads aligned, and that most mapped reads are proper pairs) are met.

Exercise #6: Yeast BWA PE alignment with Anna's script

Now that you've done everything the hard way, let's see how to do run an alignment pipeline using Anna's script.

First the setup:

mkdir -p $SCRATCH/core_ngs/align2/fastq
cd $SCRATCH/core_ngs/align2/fastq
cp /corral-repl/utexas/BioITeam/core_ngs_tools/alignment/*fastq.gz .

Before executing the script you need to have one environment variable set. We'll do it at the command line here, but you could put it in your .bashrc file.

export path_code=/work/01063/abattenh/code

Now change into the directory and call the script with no arguments to see usage

cd $SCRATCH/core_ngs/align2
$path_code/script/align/align_bwa_illumina.sh

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

1. FASTQ file – full or relative path to the FASTQ file (just the R1 fastq if paired end). Can be compressed (.gz)
2. output prefix – prefix for all the output files produced by the script. Should relate back to what the data is.
3. assembly – genome assembly to use.
   • there are pre-built indexes for some common eukaryotes (hg19, hg18, mm10, mm9, danRer7, sacCer3) that you can use
   • or provide a full path for a bwa reference index you have built somewhere
4. paired flag – 0 means single end (the default); 1 means paired end
5. hard trim length – if you want the FASTQ hard trimmed down to a specific length, supply that number here
   

Now run the pipeline. By piping the output to tee <filename> we can see the script's progress at the terminal, and it also is written to <filename>.

$path_code/script/align/align_bwa_illumina.sh ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz yeast_chip sacCer3 1 2>&1 | tee aln.yeast_chip.log

Output files

This alignment pipeline script performs the following steps:

• Hard trims FASTQ, if optionally specified (fastx_trimmer)
  
• Aligns the R1 FASTQ (bwa aln)
  
• Aligns the R2 FASTQ, if paired end alignment specified (bwa aln)
• Reports the alignments as SAM (bwa samse for single end, or bwa sampe for paired end)
  
• 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. aln.<prefix>.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>.samstats.txt – Summary alignment statistics 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 aln.yeast_chip.log

This will show something like:

..samstats file 'yeast_chip.samstats.txt' exists Thu May 28 16:36:01 CDT 2015
..samstats file file 'yeast_chip.samstats.txt' size ok Thu May 28 16:36:01 CDT 2015
---------------------------------------------------------
Cleaning up files...
---------------------------------------------------------
ckRes 0 cleanup
---------------------------------------------------------
All bwa alignment tasks completed successfully!
Thu May 28 16:36:01 CDT 2015
---------------------------------------------------------

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, suppose I have run 6 alignments and have these 6 log files:

aln.delswr1_htz1_tap1t0.log   aln.delswr1_htz1_tap1t30.log  aln.wt_htz1_tap1t15.log
aln.delswr1_htz1_tap1t15.log  aln.wt_htz1_tap1t0.log        aln.wt_htz1_tap1t30.log

I can check that all 6 completed with this command:

grep 'completed successfully' aln.*.log | wc -l

If this command returns 6, I'm done. But what if it doesn't? If you grep -v (lines that don't contain the pattern), you'll get every line in every log file except the success message line, which is not what you want at all.

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' aln.*.log

The -L option tells grep to only print the filenames that don't contain the pattern. Perfect!

Checking alignment statistics

The <prefix>.samstats.txt statistics file produced by the alignment pipeline has a lot of good information in one place. If you use cat or more to view it you'll see this:

-----------------------------------------------
             Aligner:       bwa
     Total sequences:   1184360
        Total mapped:    547664 (46.2 %)
      Total unmapped:    636696 (53.8 %)
             Primary:    547664 (100.0 %)
           Secondary:
          Duplicates:    324280 (59.2 %)
          Fwd strand:    272898 (49.8 %)
          Rev strand:    274766 (50.2 %)
           Multi hit:     18688 (3.4 %)
           Soft clip:    222451 (40.6 %)
           All match:    319429 (58.3 %)
              Indels:      6697 (1.2 %)
             Spliced:
-----------------------------------------------
       Total PE seqs:   1184360
      PE seqs mapped:    547664 (46.2 %)
        Num PE pairs:    592180
   F5 1st end mapped:    300477 (50.7 %)
   F3 2nd end mapped:    247187 (41.7 %)
     PE pairs mapped:    241180 (40.7 %)
     PE proper pairs:    236557 (39.9 %)
-----------------------------------------------
  Insert size stats for: yeast_chip
        Number of pairs: 236557 (proper)
 Number of insert sizes: 212
        Mean [-/+ 1 SD]: 215 [153 277]  (sd 62)
         Mode [Fivenum]: 223  [105 210 220 229 321]
-----------------------------------------------

Since this was a paired end alignment there is paired-end specific information reported, including insert size statistics: mean/standard deviation, mode (most common insert size value) and fivenum (min, q1, median, q3 max insert sizes).

A quick way to check alignment stats if you have run multiple alignments is again to use grep. For example, for the 6 alignment files shown earlier, running this:

grep 'Total map' *samstats.txt

will produce output like this:

delswr1_htz1_tap1t0.samstats.txt:        Total mapped:  32761761 (86.8 %)
delswr1_htz1_tap1t15.samstats.txt:        Total mapped:  33699464 (89.2 %)
delswr1_htz1_tap1t30.samstats.txt:        Total mapped:  28441655 (87.6 %)
wt_htz1_tap1t0.samstats.txt:        Total mapped:  28454847 (89.5 %)
wt_htz1_tap1t15.samstats.txt:        Total mapped:  33245627 (90.9 %)
wt_htz1_tap1t30.samstats.txt:        Total mapped:  32567026 (90.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.

Anna's alignment pipeline scripts are written to take advantage of having multiple cores on TACC nodes, and are thus designed to run with at most two pipeline commands per TACC node.

Assuming you have your alignment commands in a file called aln.cmds, here's how to create and submit a batch job for the commands.

launcher_creator.py -n aln -j aln.cmds -t 12:00:00 -q normal -w 2
sbatch aln.slurm
showq -u

Note the maximum run time specified here is 12 hours (-t 12:00:00). This is a reasonable value for a higher eukaryote with 20-40 M reads, and is way more than a yeast alignment would need (~ 4 hours). For very deeply sequenced eukaryotes (e.g. human genome re-sequencing with hundresd of millions of reads), you may want to specify the maximum job time of 48 hours.

Exercise: What would alignment commands look like if you were putting it in a batch system .cmds file?

$path_code/script/align/align_bwa_illumina.sh ./fastq/Sample_Yeast_L005_R1.cat.fastq.gz yeast1 sacCer3 1 2>&1 > aln.yeast1.log
$path_code/script/align/align_bwa_illumina.sh ./fastq/Sample_ABCDE_L005_R1.cat.fastq.gz yeast2 sacCer3 1 2>&1 > aln.yeast2.log