View source for 6.1 Genome indexing for bowtie

To  map  reads  to  a  reference  genome,  all  mapping  tools  require  the  genome  to  be  converted  into  particular  structures  for  quick  access  and  to  keep  the  memory  footprint  small. In this tutorial we will use '''tophat2''' for the transcriptome reconstruction, which  in turn requires '''bowtie2''' to run for the read mapping to the genome. '''Bowtie2''' needs  the  genome  to  be  indexed  using  the  BurrowsAWheeler  transform,  and  provides  a  tool  ('''bowtie2-build''')  to  obtain  this  transformation  starting  from  the  genome  sequence  stored  in  a  text  file  in  fasta  format.  You  will  also  try  another  mapper,  '''STAR''',  which  performs quite well and can be easily integrated with '''cufflinks/cuffdif''', as you will  see later.

Let’s first look at the bowtie2-build options: 
 
 $ module load bowtie/2.2.3
 $ bowtie2-build --help
 
The general usage is: 
 
Usage: 
 bowtie2-build [options]* <reference_in> <bt2_index_base>
 
You  need  to  specify  the  genomic  sequence  file  ('''reference_in''')  and  a  label  to  identify 
the index ('''bt2_index_base'''), which will be the prefix of all files written by '''bowtie2-build'''.
Copy a folder containing the genomic sequence with the following command: 
 
 $ cp –r /ibers/repository/public/courses/Rna-Seq/genome .

Check the content of the genome folder that you just copied: 
 
 $ ls –l genome

Have a look at the file with the commands more or head. This is a file in fasta format,  widely  used  for  storing  nucleotide  or  amino  acid  sequences.
The  first  row,  and  any  row starting  with  >,  is  the  name  of  the  sequence.  For  this  tutorial,  only  the  zebrafish chromosome 12 is in this file, and the reads subset that we are using is already filtered for mapping on it.  
Try to count the number of nucleotides in the chromosome 12: knowing that each row in the fasta file contains 60 nucleotides, how would you know the chromosome size?  
To  keep  things  well  organized,  you  can  create  a  folder  to  store  the  index  (name  it  for 
example '''bt2Index''' or any other name you like):
 $ mkdir bt2Index
Once  you  created  the  folder,  let’s  open  a  text  file  with  a  text  editor.  You  can  call  the  file  any name you want, for example '''bowtie2index.sh'''. Write at the beginning of the file  the header for the scheduler, and then load the required modules: 

 #$ -S /bin/sh
 #$ -cwd
 #$ -q amd.q,large.q,intel.q
 #$ -l h_vmem=10G
 #$ -e run_bowtie_index.e
 #$ -N run_bowtie_index
 #$ -o run_bowtie_index.o
 module load bowtie/2.2.3

Now write the Unix command to instruct the script to go in  the folder in which to execute the commands. For my account, I would write in the script:   
 cd /ibers/ernie/home/vpl/zebrafish/bt2Index 
 
If  you  are  not  sure  the  full  path,  i.e.  where  your  folder  is,  check  the  path  of  the  folder  using '''pwd'''. 
Now you can write the command to index the genome: 
 bowtie2-build -f ../genome/zebrafish.fa Z12

'''-f''': specify that the genome is in fast format
 
'''../genome/zebrafish.fa''': path  to  the  genome  fasta  file  (remember that the command .. indicates the parent folder, i.e. the one right above that  in which you currently are)

'''Z12''': the prefix that I chose for all files written by bowtie2-build (but you can chose a  different one)

Save and close the file. Now you can submit it to the scheduler using the '''qsub''' command,  like this: 
 
 $ qsub bowtie2index.sh
The job is submitted to the scheduler, and the job identifier is written on screen. You can  check the status of the job using the '''qstat''' command. This command returns the status  of all jobs from all users, hence to check only your job you can write something like this: 
 
 $ qstat | grep <your user name>

Remember that the operator '''|''' pipes the output of a command ('''qstat''' in this case) to  another command ('''grep''' in this case). The command '''grep''' will filter the '''qstat''' output  to select only those lines containing your user name.  Once the job is finished, look at the content of the folder (with '''ls –l'''). The bunch of files  having extension '''.bt2''' and prefix '''Z12''' contain the '''index''' of the genome (or in your case  only of the chromosome 12) in a format that '''bowtie2''' can use.   Now we are ready to map the reads onto the chromosome 12.