Examples of Read Mapping with Karma and BWA

#  Some instructions for read mapping and variant calling using the 
#  University of Michigan tools and procedures.  Please view this as source.

	#  -- Paul Anderson and Tom Blackwell, December 11, 2009  --

#  These instructions will cover three components:

#  (1)  read mapping using karma
#  (2)  read mapping using bwa
#  (3)  samtools processing from .sam to .glf

#  These instructions are specifically with reference to the CEU trio 
#  chromosome 20 test data set that has been distributed.  Please 
#  see some general discussion at the beginning of item 2.  Everything 
#  that's not commented out should be runnable code, although you 
#  will need to alter the directory and file names.

#  (1)  Command line procedure for karma read mapping of the CEU trio 
#  chromosome 20 test data set.  Paul Anderson writes:  

#  Start with four subdirectories:  

#  indiv  containing single and paired end Illumina .fastq sequence files, 
#  ab     containing AB SOLiD color space reads as in the test data set, 
#  k.ref  containing the gzipped human genome reference sequence,  and 
#  k.out  an empty directory for the resulting karma .sam and .stats files.

#  First, build karma's binary word index files:

cd k.ref zcat human_g1k_v37.fasta.gz > human_g1k_v37.fa karma --reference human_g1k_v37.fa --createIndex --occurrenceCutoff 5000

#  then map Illumina paired end or single end reads as, for example:

cd ../k.out karma --reference ../k.ref/human_g1k_v37.fa --pairedReads --maxInsert 2000 \ ../indiv/SRR010941_1.recal.fastq.gz ../indiv/SRR010941_2.recal.fastq.gz

karma --reference ../k.ref/human_g1k_v37.fa --maxInsert 2000 \ ../indiv/SRR010936.recal.fastq.gz

#  The read lengths in the Illumina data I saw are long enough so that the 
#  only parameter worth tweaking might be  --occurrenceCutoff  in the index 
#  structures.  I will need to study the mapping performance to see if this 
#  makes any appreciable difference in speed.

#  LS 454 fastq files can also be mapped using the same index, but the 
#  mapping success rate will be much lower.

#  For AB Solid data, a color space reference and set of index files must be 
#  created.  For the lengths of reads I see, we will create two -- one with a 
#  12-mer index, and the other with a 15-mer.

cd ../k.ref ln -s human_g1k_v37.fa human_g1k_v37_12CS.fa ln -s human_g1k_v37.fa human_g1k_v37_15CS.fa karma --reference human_g1k_v37_12CS.fa --colorSpace --createIndex --wordSize 12 karma --reference human_g1k_v37_15CS.fa --colorSpace --createIndex --wordSize 15

#  Then, short color space reads (25/26-mers) should be mapped using the 12-mer index:

cd ../k.out karma --reference ../k.ref/human_g1k_v37.fa \ --csreference ../k.ref/human_g1k_v37_12CS.fa \ --colorSpace --pairedReads --maxInsert 2000 \ ../ab/TG150_1.color.space.fastq.gz \ ../ab/TG150_2.color.space.fastq.gz

#  To map longer color space reads (>30-mer), use the bigger index:

karma --reference ../k.ref/human_g1k_v37.fa \ --csreference ../k.ref/human_g1k_v37_15CS.fa \ --colorSpace --pairedReads --maxInsert 2000 \ ../ab/TG152_1.color.space.fastq.gz \ ../ab/TG152_2.color.space.fastq.gz

#  You can use the shorter reference to map any length reads, but it is slower.  
#  The 15-mer color space index (human_g1k_v37_15CS.fa) is a good compromise 
#  for mapping reads 30 or more bases, and for machines with 24GB of RAM.

#  The maxInsert of 2000 is somewhat arbitrary - you want a number that includes 
#  the bulk of the paired end read insert distances you are interested in.  More 
#  is slower, so you don't want an arbitrarily large value.

#  (2)  Command line procedure for bwa read mapping of the CEU trio 
#  chromosome 20 test data set.

#  As a general principle, one wants to keep all sequence data separated 
#  by sequencing technology and by individual throughout the entire mapping 
#  and variant calling process.  Different technologies may have different 
#  requirements for read mapping and different characteristics for variant 
#  calling.  This can be done either using the directory structure or with a 
#  file tracking database.  

#  As a further complication, the Broad Institute Illumina sequencing runs 
#  in the current test data set benefit from removing leading and trailing Ns 
#  before read mapping with bwa.  All Illumina data for the CEU trio that 
#  have SRR identifiers are Broad Institute runs.  (There are no Broad 
#  runs for the YRI trio.  There, SRR identifiers will indicate either Beijing 
#  or Wash U runs.)  I have used exactly the same bwa command lines with 
#  or without N-trimming.

#  In what follows, I will use somewhat simplified directory and file names.
#  At this stage, it's easiest to have the command line prompt remain at top 
#  level in the directory structure, and refer to all files using relative 
#  pathnames, relative to the location of the prompt.  At the start, suppose 
#  that there are two subdirectories:  'bwa.ref'  containing the gzipped human 
#  genome reference sequence, and  'indiv'  containing many .fastq files for 
#  the target individual generated by the appropriate sequencing technology.  
#  The bwa program has an inconvenient habit of writing to std.err.  I routinely 
#  redirect that to a log file so that I can keep on working while bwa runs in 
#  the background.  

#  In general, bwa read mapping requires initial indexing of the genome 
#  reference sequence, followed by two passes for each .fastq file.  The 
#  software components of this process are:  

#  bwa  index  --  Do this only when a new version of the genome reference 
#                  sequence is released.  Takes just under two hours.
#  bwa   aln   --  Run this on every .fastq file individually.  Highly 
#                  variable timings -- takes between 10 minutes and many 
#                  hours per .fastq file.  Better data runs quicker.
#  bwa  samse  --  Converts a single .fastq / .sai pair to .sam alignment 
#                  format.  Usually under 1 minute per file.
#  bwa  sampe  --  Converts paired end .fastq / .sai pairs (four files total) 
#                  to a single .sam alignment file.  1 - 2 minutes per run.  
#                  (In general, the CEU trio chromosome 20 test data set 
#                  contains VERY small .fastq files.  Chromosome 20 is just 
#                  over 2% of the entire genome.)

#  Sample command lines.  This is csh syntax.  

#  The parameter string  " -n 0.002  -M 7  -R 25 "  used in bwa aln below 
#  is just my initial guess.  I have made an entire run using this string.  
#  I might make another run using Heng Li's default values for all three 
#  parameters and compare the results, but I have no answers from this yet.  

bwa index -a bwtsw -p bwa.ref/ncbi.v37.ref \ bwa.ref/human_g1k_v37.fasta.gz >>& logfile

mkdir bwa.sai bwa.sam

( nice +20 bwa aln -n 0.002 -M 7 -R 25 bwa.ref/ncbi.v37.ref \ indiv/SRRx_1.fastq.gz > bwa.sai/SRRx_1.sai ) >>& logfile

#  ... and the same for  SRRx_2.fastq.gz  or  SRRx.fastq.gz.
#  Then, depending on whether these are paired end reads or not, either:

( nice +20 bwa sampe bwa.ref/ncbi.v37.ref \ bwa.sai/SRRx_1.sai bwa.sai/SRRx_2.sai \ indiv/SRRx_1.fastq.gz indiv/SRRx_2.fastq.gz \ > bwa.sam/SRRx.pair.sam ) >>& logfile

#  or, for single end reads:

( nice +20 bwa samse bwa.ref/ncbi.v37.ref \ bwa.sai/SRRx.sai indiv/SRRx.fastq.gz \ > bwa.sam/SRRx.single.sam ) >>& logfile

#  For LS 454 sequence data, bwa read mapping is a one-step process.  
#  This uses the same genome reference sequence index as for Illumina 
#  data above, and does not use paired end information in the mapping.  
#  (The directory  bwa.sam  shown here should be separate from that 
#  created above for Illumina data.  Similarly for AB SOLiD below.)  

( nice +20 bwa dbwtsw bwa.ref/ncbi.v37.ref indiv/SRRx.fastq.gz \ > bwa.sam/SRRx.sam ) >>& logfile

#  For AB SOLiD data, one must build a separate index structure for 
#  the genome reference sequence, and (completely undocumented) one 
#  must rewrite all of the .fastq files replacing  0,1,2,3,"."  with 
#  A,C,G,T,N  in every line of color space sequence and omitting the 
#  first two characters from each line of converted sequence and from 
#  the base call quality strings.  (An awk script does this conversion 
#  really quickly.)

bwa index -a bwtsw -p bwa.ref/ncbi.color.ref -c \ bwa.ref/human_g1k_v37.fasta.gz >>& logfile

mkdir bwa.sai bwa.sam

( nice +20 bwa aln -n 0.002 -M 7 -R 25 -o 0 -c bwa.ref/ncbi.color.ref \ indiv/TG145_1.cvt.gz > bwa.sai/TG145_1.sai ) >>& logfile

#  ... and the same for  TG145_2.cvt.gz  or  TG145.cvt.gz.
#  Then, depending on whether these are paired end reads or not, either:

( nice +20 bwa sampe bwa.ref/ncbi.color.ref \ bwa.sai/TG145_1.sai bwa.sai/TG145_2.sai \ indiv/TG145_1.cvt.gz indiv/TG145_2.cvt.gz \ > bwa.sam/SRRx.pair.sam ) >>& logfile

#  or:

( nice +20 bwa samse bwa.ref/ncbi.color.ref \ bwa.sai/TG145.sai indiv/TG145.cvt.gz \ > bwa.sam/SRRx.single.sam ) >>& logfile

#  (3)  Here is the process which builds the .glf format files used for 
#  SNP calling out of .sam format files for individual sequencing runs 
#  produced using either read mapping algorithm.  This will use only 
#  samtools utilities and contains nothing specific to either read mapper.  

#  For illustration, I will assume two of the directories left over from 
#  bwa read mapping in the preceding section.  These are:  

#  bwa.ref  containing files  human_g1k_v37.fasta.gz, human_g1k_v37.fasta.fai, and 
#  bwa.sam  containing all of the single-end and paired-end .sam files generated 
#  	 	 	for one individual and sequencing technology.

#  The full set of command lines is:  (Note csh syntax again -- and this time, 
#     I will change directories for convenience in file naming.)

mkdir bwa.bam indiv.bam cd bwa.bam set minq=17

foreach file ( ../bwa.sam/*.sam )

  set unsorted=`basename $file .sam`.nosort.bam
  samtools view -bhuS  -o $unsorted  -q $minq  $file	>>& logfile
  samtools sort  $unsorted  `basename $file .sam`	>>& logfile
  rm $unsorted
  end

unset minq

#  This is the point where one would insert steps to recalibrate 
#  base call quality values, check genotype identity or remove 
#  duplicate sequence reads.  Each of these steps is specific to 
#  an individual sequencing run.

cd ../indiv.bam samtools merge person.bam ../bwa.bam/*.bam samtools index person.bam samtools view person.bam 20 | samtools pileup \ -f ../bwa.ref/human_g1k_v37.fasta.gz \ -t ../bwa.ref/human_g1k_v37.fasta.fai \ -g - > person.glf