Bwa

From wiki
Revision as of 13:17, 13 March 2017 by Rf (talk | contribs)
Jump to: navigation, search

Introduction

Samtools originator Heng Li coded this aligner.

Usage

As with samtools, bwa also went through some re-structuring, so that it has an old-style tw-step (aln and sam{s,p}e) usage, characterised by the following typical sequence of commands:

bwa index reference.fa
bwa aln -I -t 8 reference.fa s_1.txt > out.sai
bwa samse reference.fa out.sai s_1.txt > out.sam
samtools view -bSu out.sam | samtools sort -  out.sorted

There is a more modern usage which consists of just one step: bwa mem. Again the ouput must be directed to a sam file name of your choosing. However the aln and samse or sampe methods are still useful for certain conditions. Please consult documentation.

Note how there are no option switches. The prefix for the reference must come first, and the reads must come second.

Also note that if there are several readsets, it is best to align and obtain a sam/bam file for each individually (of course, not pairs, the pairs count as one readset). Canonical usage

bwa mem index_prefix [input_reads.fastq|input_reads_pair_1.fastq input_reads_pair_2.fastq] [options]

Indexing

Before alignment, indexing the reference is necessary. It is a relatively simple task:

bwa <name_of_reference_file>


In the case of of reference filename being "unplaced.scaf.fa", the output files will be:

  • unplaced.scaf.fa.amb, a text file
  • unplaced.scaf.fa.ann, a text file
  • unplaced.scaf.fa.bwt, a binary file
  • unplaced.scaf.fa.pac, a binary file
  • unplaced.scaf.fa.sa, a binary file

There is only one argument in this case, so when bwa indexes the reference, it will use the whole filename and generate output index files with extensions added onto this name. A second argument is allowed which allows a custom prefix to be used for these output index files.

bwa index input_reference.fasta index_prefix

This index_prefix name can then used for the actual alignment step. Otherwise the filename of the reference file will be used not as a full names but in a name prefix context.

Indexing may take a long time with large reference files. Here is some example output

[bwa_index] Pack FASTA... 150.69 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=4745806978, availableWord=345932252
[BWTIncConstructFromPacked] 10 iterations done. 99999986 characters processed.
[BWTIncConstructFromPacked] 20 iterations done. 199999986 characters processed.
.
.
.
[BWTIncConstructFromPacked] 530 iterations done. 4719199682 characters processed.
[BWTIncConstructFromPacked] 540 iterations done. 4741303618 characters processed.
[bwt_gen] Finished constructing BWT in 543 iterations.
[bwa_index] 2470.36 seconds elapse.
[bwa_index] Update BWT... 129.37 sec
[bwa_index] Pack forward-only FASTA... 302.45 sec
[bwa_index] Construct SA from BWT and Occ... 1220.18 sec
[main] Version: 0.7.12-r1039
[main] CMD: bwa index unplaced.scaf.fa
[main] Real time: 5398.987 sec; CPU: 4273.064 sec


Actual alignment step

bwa mem index_prefix input_reads_pair_1.fastq input_reads_pair_2.fastq

In this case we have paired read-files. With single reads of course only one name would be required.