Page History

...

Based on the usage description, we only need to specify two things:

• the The name of the FASTA file  
• whether Whether to use the bwtsw or is algorithm for indexing

...

mkdir -p $SCRATCH/core_ngs/alignment/fastq
mkdir -p $SCRATCH/core_ngs/references/fasta
cp $CORENGS/alignment/*fastq.gz $SCRATCH/core_ngs/alignment/fastq/
cp $CORENGS/references/*.*fa $SCRATCH/core_ngs/references/fasta/

...

First, the single quotes around the pattern – this tells the bash shell to pass the exact string contents to grep.

As part of its friendly command line parsing and evaluation, the shell will often look for special characters on the command line that mean something to it (for example, the $ in front of an environment variable name, like in $SCRATCH). Well, regular expressions treat the $ specially too – but in a completely different way! Those single quotes tell the shell "don't look inside here for special characters – treat this as a literal string and pass it to the program". The shell will obey, will strip the single quotes off the string, and will pass the actual pattern, ^>, to the grep program. (Note that the shell does look inside double quotes ( " ) for certain special signals, such as looking for environment variable names to evaluate. Read more about Quoting in the shell.)

...

Now, we're ready to execute the actual alignment, with the goal of initially producing a SAM file from the input FASTQ files and reference. First prepare a directory for this exercise and link the sacCer3 reference directories there (this will make our commands more readable).

# Copy the pre-built references
mkdir -p $SCRATCH/core_ngs/references
rsynccp -avrP $CORENGS/references/*.fa $SCRATCH/core_ngs/references/fasta/

# Get the FASTQ to align
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz $SCRATCH/core_ngs/alignment/fastq/

...

For a basic alignment like this, we can just go with the default alignment parameters.

Also note Note that bwa writes its (binary) output to standard output by default, so we need to redirect that to a .sai file.

...

idev -p normal -m 180 -N 1 -n 24 -A UT-2015-05-18 --reservation=CCBBintro_NGS

For simplicity, we will just execute these commands directly, one at a time. Each command should only take few minutes and you will see bwa's progress messages in your terminal.

...

bwa sampe sacCer3/sacCer3.fa yeast_R1.sai yeast_R2.sai \
  fastq/Sample_Yeast_L005_R1.cat.fastq.gz \
  fastq/Sample_Yeast_L005_R2.cat.fastq.gz > yeast_pairedend.sam

...

grep -P -c '^HWI' yeast_pairedend.sam

grep -P -v -c '^@' yeast_pairedend.sam

...

• Do both R1 and R2 reads have separate alignment records?
• Does the SAM file contain both mapped and un-mapped reads?
• What is the order of the alignment records in this SAM file?
  
  

...

You may have noticed that some alignment records contain contig names (e.g. chrV) in field 3 while others contain an asterisk ( * ). Usually the The * means the record didn't map. (This isn't always true – later you'll see how to properly distinguish between mapped and unmapped reads using samtools.) We're going to use this heuristic along with cut to see about how many records represent aligned sequences. (Note this is not the strictly correct method of finding unmapped reads because not all unmapped reads have an asterisk in field 3. Later you'll see how to properly distinguish between mapped and unmapped reads using samtools.)

First we need to make sure that we don't look at fields in the SAM header lines. We're going to end up with a series of pipe operations, and the best way to make sure you're on track is to enter them one at a time piping to head:

# the ^HWI^@ pattern matches lines starting with HWI@ (the start of all read names in column 1)only header lines), 
# and -v says output lines that don't match
grep -v -P '^HWI^@' yeast_pairedend.sam | head

...

grep -v -P '^HWI^@' yeast_pairedend.sam | cut -f 3 | head

...

grep -v -P '^HWI^@' yeast_pairedend.sam | cut -f 3 | grep -v '*' | head

...

grep -v -P '^HWI^@' yeast_pairedend.sam | cut -f 3 | grep -v '*' | wc -l

...