/
bwa_mapping.txt
159 lines (117 loc) · 6.78 KB
/
bwa_mapping.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
Mapping with bwa
================
One of the most common operations when dealing with NGS data is mapping sequences
to other sequences, and most commonly, mapping reads to a reference. Because
short read alignment is so common, there is a deluge of programs available for doing
it, and it is well-understood problem. Most of the efficient programs work by first
building an *index* of the sequence that will be mapped to, which allows for
extremely fast lookups, and then running a separate module which uses that index
to align short sequences against the reference. Alignment is used for visualization,
coverage estimation, SNP calling, expression analysis, and a slew of other problems.
For a high-level overview, try `this NCBI review <http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2836519/>`__.
Getting the Dependencies
````````````````````````
bwa is one of the many available read mappers. bwa is a good "reference" version
of a common alignment algorithm based on the Burrows-Wheeler transform; those
of you looking to dive into a relatively thick paper on it can read about it
`here <http://bioinformatics.oxfordjournals.org/content/25/14/1754.short>`__.
Firstly, let's download bwa::
cd /root
curl -L "http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.5a.tar.bz2?r=http%3A%2F%2Fsourceforge.net%2Fprojects%2Fbio-bwa%2Ffiles%2F&ts=1379347638&use_mirror=softlayer-dal" > bwa-0.7.5a.tar.bz2
And then decompress it::
tar xjf bwa-0.7.5a.tar.bz2
Now, change into the directory with the source code and compile it::
cd bwa-0.7.5a/
make
Finally, move the compiled executable to a location where the operating system
can find it::
cp bwa /usr/local/bin
So what did we actually just do? bwa is offered as an open-source program -- that
is, the actual code which defines its functionality is freely available. However,
bwa is written in C, which, unlike Python, is not "interpreted." Instead, the
source code needs to be parsed and converted into machine language which can be
run on the processor. This allows compiled programs to run much more efficiently
than their interpreted counterparts, at the cost of needing to be recompiled for
different systems. Although pre-compiled executables can be found for many programs,
Unix programs are often distributed as source. 'make' is a program for tracking
dependencies amongst files, and is used to manage the compilation of larger
projects with many files.
'make' is also very useful for setting up pipelines, as it is not limited to
source files -- it can be used for any situation where one program needs to output
of another to run. You can read about it `here <http://www.gnu.org/software/make/>`__ if you're interested, though
it's outside the scope of this tutorial.
Getting the Data
````````````````
We'll be using the data we downloaded during the reads and quality control session;
if you missed that, you'll eventually want to run through it in :doc:`reads_and_qc`,
but for now, you can just grab the quality-controlled data with::
cd /mnt
curl -O https://s3.amazonaws.com/public.ged.msu.edu/SRR390202.pe.qc.fq.gz
curl -O https://s3.amazonaws.com/public.ged.msu.edu/SRR390202.se.qc.fq.gz
You'll also need a reference genome, which can be acquired with::
curl -O http://ftp.ncbi.nlm.nih.gov/genomes/Bacteria/Escherichia_coli_O104_H4_2011C_3493_uid176127/NC_018658.fna
Mapping the Reads
`````````````````
To speed up the demonstration, we will just map a subset of the reads rather than
the entire file, which is somewhat large (though small compared to many datasets).
The head command outputs the first n lines of a file, by default 4::
cd /mnt/ecoli
gunzip -c ../SRR390202.pe.qc.fq.gz | head -400000 > ecoli_pe.fq
gunzip -c ../SRR390202.se.qc.fq.gz | head -400000 > ecoli_se.fq
We've got our reads and a reference, so we're ready to get started. First, we
build an index of the reference genome using bwa::
mv ../NC_018658.fna .
bwa index -a bwtsw NC_018658.fna
The -a flag tells bwa which indexing algorithm to use. The program will automatically
output some files with set extensions, which the main alignment program knows the
format of. Thus, we run the alignment like so::
bwa mem -p NC_018658.fna ecoli_pe.fq > aln.x.ecoli_NC_018658.sam
which aligns the left and right reads against the reference, and outputs them
to the given SAM file. SAM is a common format for alignments which is understood
by many programs, along with BAM. It's often useful to have both, so we'll use
a utility called samtools to produce a sorted BAM file as well. First,
install samtools::
cd /mnt
curl -O -L http://sourceforge.net/projects/samtools/files/samtools/0.1.18/samtools-0.1.18.tar.bz2
tar xvfj samtools-0.1.18.tar.bz2
cd samtools-0.1.18
make
cp samtools /usr/local/bin
cd misc/
cp *.pl maq2sam-long maq2sam-short md5fa md5sum-lite wgsim /usr/local/bin/
cd ..
cd bcftools
cp *.pl bcftools /usr/local/bin/
Then, run samtools to do the conversion::
cd /mnt/ecoli
samtools view -uS aln.x.ecoli_NC_018658.sam > aln.x.ecoli_NC_018658.bam
samtools sort aln.x.ecoli_NC_018658.bam aln.x.ecoli_NC_018658.bam.sorted
samtools index aln.x.ecoli_NC_018658.bam.sorted.bam
For additional resources on these tools, check out:
* the bwa `manual <http://bio-bwa.sourceforge.net/bwa.shtml>`__
* info on `samtools <http://samtools.sourceforge.net/>`__
* the SAM `format spec <http://samtools.sourceforge.net/SAM1.pdf>`__
Visualizing your Data with Tablet
`````````````````````````````````
First, install dropbox using :doc:`amazon/installing-dropbox`.
Now, copy your mapping files and the reference to your dropbox folder::
cp aln.x.ecoli_NC_018658.bam.sorted.bam* /root/Dropbox
cp NC_018658.fna /root/Dropbox
Although you can do many things with your alignments, one useful thing is to simply
view them through a graphical interface. To demonstrate this, we'll use a program
called Tablet, which can be downloaded `here <http://bioinf.scri.ac.uk/tablet/download.shtml>`__.
Tablet claims to run on Windows, Linux, and OSX, though I have only tested it out
on Linux. Because this is a GUI-driven program, we'll be running it on our
local machines instead of our EC2 instances. So, go ahead and grab the appropriate
version for your system, and install it.
Once you have it installed, open it up. You'll want to load your mapping file
and reference genome:
.. image:: images/tablet_open.png
:width: 75%
You'll then get some loading bars, and potentially an error about the indexing file
which can be ignored. You need to select a contig on the left to view; ecoli has
a very good reference, and is only one contig:
.. image:: images/tablet_view.png
:width: 75%
You can move left and right along the contig, as well as zoom. Tablet can view
other information like gene structure, but we won't get into that.