5. Read mapping

5.1. Preface

In this section we will use our skill on the command-line interface to map our reads from the evolved line to our ancestral reference genome.

Note

You will encounter some To-do sections at times. Write the solutions and answers into a text-file.

5.2. Overview

The part of the workflow we will work on in this section can be viewed in Fig. 5.1.

../_images/workflow3.png

Fig. 5.1 The part of the workflow we will work on in this section marked in red.

5.3. Learning outcomes

After studying this section of the tutorial you should be able to:

  1. Explain the process of sequence read mapping.
  2. Use bioinformatics tools to map sequencing reads to a reference genome.
  3. Filter mapped reads based on quality.

5.4. Before we start

Lets see how our directory structure looks so far:

cd ~/analysis
# create a mapping result directory
mkdir mappings
ls -1F
assembly/
data/
mappings/
SolexaQA/
SolexaQA++
trimmed/
trimmed-fastqc/
trimmed-solexaqa/

5.5. Mapping sequence reads to a reference genome

We want to map the sequencing reads to the ancestral reference genome we created in the section Genome assembly. We are going to use the quality trimmed forward and backward DNA sequences of the evolved line and use a program called BWA to map the reads.

Todo

  1. Discuss briefly why we are using the ancestral genome as a reference genome as opposed to a genome for the evolved line.

5.5.1. Installing the software

We are going to use a program called BWA fo map our reads to a genome.

It is simple to install and use.

source activate ngs
conda install samtools
conda install bamtools
conda install bedtools
conda install bowtie2
conda install bwa
conda install picard

5.6. Bowtie2

5.6.1. Overview

Bowtie2 is a short read aligner, that can take a reference genome and map single- or paired-end data to it [TRAPNELL2009]. It requires an indexing step in which one supplies the reference genome and Bowtie2 will create an index that in the subsequent steps will be used for aligning the reads to the reference genome. The general command structure of the Bowtie2 tools we are going to use are shown below:

# bowtie2 help
bowtie2-build

# indexing
bowtie2-build genome.fasta PATH_TO_INDEX_PREFIX

# paired-end mapping
bowtie2 -X 1000 -x PATH_TO_INDEX_PREFIX -1 read1.fq.gz -2 read2.fq.gz -S aln-pe.sam
  • -X: Adjust the maximum fragment size (length of paired-end alignments + insert size) to 1000bp. This might be useful if you do not know the exact insert size of your data. The Bowtie2 default is set to 500 which is often considered too short.

5.6.2. Creating a reference index for mapping

Todo

Create an Bowtie2 index for our reference genome assembly. Attention! Remember which file you need to submit to Bowtie2.

Hint

Should you not get it right, try the commands in Code: Bowtie2 indexing.

Note

Should you be unable to run Bowtie2 indexing on the data, you can download the index from Downloads. Unarchive and uncompress the files with tar -xvzf bowtie2-index.tar.gz.

5.6.3. Mapping reads in a paired-end manner

Now that we have created our index, it is time to map the filtered and trimmed sequencing reads of our evolved line to the reference genome.

Todo

Use the correct bowtie2 command structure from above and map the reads of the evolved line to the reference genome.

Hint

Should you not get it right, try the commands in Code: Bowtie2 mapping.

Note

Bowtie2 does give very cryptic error messages without telling much why it did not want to run. The most likely reason is that you specified the paths to the files and result file wrongly. Check this first. Use tab completion a lot!

5.7. BWA

Attention

If the mapping did not succeed with Bowtie2. We can use the aligner BWA explained in this section. If the mapping with Bowtie2 did work, you can jump this section.

5.7.1. Overview

BWA is a short read aligner, that can take a reference genome and map single- or paired-end data to it. It requires an indexing step in which one supplies the reference genome and BWA will create an index that in the subsequent steps will be used for aligning the reads to the reference genome. The general command structure of the BWA tools we are going to use are shown below:

# bwa index help
bwa index

# indexing
bwa index reference-genome.fa

# bwa mem help
bwa mem

# single-end mapping
bwa mem reference-genome.fa reads.fq > aln-se.sam

# paired-end mapping
bwa mem reference-genome.fa read1.fq read2.fq > aln-pe.sam

5.7.2. Creating a reference index for mapping

Todo

Create an BWA index for our reference genome assembly. Attention! Remember which file you need to submit to BWA.

Hint

Should you not get it right, try the commands in Code: BWA indexing.

Note

Should you be unable to run BWA indexing on the data, you can download the index from Downloads. Unarchive and uncompress the files with tar -xvzf bwa-index.tar.gz.

5.7.3. Mapping reads in a paired-end manner

Now that we have created our index, it is time to map the filtered and trimmed sequencing reads of our evolved line to the reference genome.

Todo

Use the correct bwa mem command structure from above and map the reads of the evolved line to the reference genome.

Hint

Should you not get it right, try the commands in Code: BWA mapping.

5.8. The sam mapping file-format

Bowtie2 and BWA will produce a mapping file in sam-format. Have a look into the sam-file that was created by either program. A quick overview of the sam-format can be found here and even more information can be found here. Briefly, first there are a lot of header lines. Then, for each read, that mapped to the reference, there is one line.

The columns of such a line in the mapping file are described in Table 5.1.

Table 5.1 The sam-file format fields.
Col Field Description
1 QNAME Query (pair) NAME
2 FLAG bitwise FLAG
3 RNAME Reference sequence NAME
4 POS 1-based leftmost POSition/coordinate of clipped sequence
5 MAPQ MAPping Quality (Phred-scaled)
6 CIAGR extended CIGAR string
7 MRNM Mate Reference sequence NaMe (‘=’ if same as RNAME)
8 MPOS 1-based Mate POSition
9 ISIZE Inferred insert SIZE
10 SEQ query SEQuence on the same strand as the reference
11 QUAL query QUALity (ASCII-33 gives the Phred base quality)
12 OPT variable OPTional fields in the format TAG:VTYPE:VALUE

One line of a mapped read can be seen here:

M02810:197:000000000-AV55U:1:1101:10000:11540   83      NODE_1_length_1419525_cov_15.3898       607378  60      151M    =       607100  -429    TATGGTATCACTTATGGTATCACTTATGGCTATCACTAATGGCTATCACTTATGGTATCACTTATGACTATCAGACGTTATTACTATCAGACGATAACTATCAGACTTTATTACTATCACTTTCATATTACCCACTATCATCCCTTCTTTA FHGHHHHHGGGHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHGHHHHHGHHHHHHHHGDHHHHHHHHGHHHHGHHHGHHHHHHFHHHHGHHHHIHHHHHHHHHHHHHHHHHHHGHHHHHGHGHHHHHHHHEGGGGGGGGGFBCFFFFCCCCC NM:i:0  MD:Z:151        AS:i:151        XS:i:0

It basically defines, the read and the position in the reference genome where the read mapped and a quality of the map.

5.9. Mapping post-processing

5.9.1. Fix mates and compress

Because aligners can sometimes leave unusual SAM flag information on SAM records, it is helpful when working with many tools to first clean up read pairing information and flags with SAMtools. We are going to produce also compressed bam output for efficient storing of and access to the mapped reads.

samtools fixmate -O bam mappings/evolved-6.sam mappings/evolved-6.fixmate.bam
  • -O bam: specifies that we want compressed bam output

Attention

The step of sam to bam-file conversion might take a few minutes to finish, depending on how big your mapping file is.

We will be using the SAM flag information later below to extract specific alignments.

Hint

A very useful tools to explain flags can be found here.

Once we have bam-file, we can also delete the original sam-file as it requires too much space.

rm mappings/evolved-6.sam

5.9.2. Sorting

We are going to use SAMtools again to sort the bam-file into coordinate order:

# convert to bam file and sort
samtools sort -O bam -o mappings/evolved-6.sorted.bam mappings/evolved-6.fixmate.bam
  • -o: specifies the name of the output file.
  • -O bam: specifies that the output will be bam-format

5.9.3. Remove duplicates

In this step we remove duplicate reads. The main purpose of removing duplicates is to mitigate the effects of PCR amplification bias introduced during library construction. It should be noted that this step is not always recommended. It depends on the research question. In SNP calling it is a good idea to remove duplicates, as the statistics used in the tools that call SNPs sub-sequently expect this (most tools anyways). However, for other research questions that use mapping, you might not want to remove duplicates, e.g. RNA-seq.

picard MarkDuplicates REMOVE_DUPLICATES=true METRICS_FILE=mappings/evolved-6.marked_dup_metrics.txt INPUT=mappings/evolved-6.sorted.bam OUTPUT=mappings/evolved-6.sorted.dedup.bam

Todo

Figure out what "PCR amplification bias" means.

Note

Should you be unable to do the post-processing steps, you can download the mapped data from Downloads.

5.10. Mapping statistics

5.10.1. Stats with SAMtools

Lets get an mapping overview:

samtools flagstat mappings/evolved-6.sorted.dedup.bam

Todo

Look at the mapping statistics and understand their meaning. Discuss your results. Explain why we may find mapped reads that have their mate mapped to a different chromosome/contig? Can they be used for something?

For the sorted bam-file we can get read depth for at all positions of the reference genome, e.g. how many reads are overlapping the genomic position.

samtools depth mappings/evolved-6.sorted.dedup.bam | gzip > mappings/evolved-6.depth.txt.gz

Todo

Extract the depth values for contig 20 and load the data into R, calculate some statistics of our scaffold.

zcat mappings/evolved-6.depth.txt.gz | egrep '^NODE_20_' | gzip >  mappings/NODE_20.depth.txt.gz

Now we quickly use some R to make a coverage plot for contig NODE20. Open a R shell by typing R on the command-line of the shell.

x <- read.table('mappings/NODE_20.depth.txt.gz', sep='\t', header=FALSE,  strip.white=TRUE)

# Look at the beginning of x
head(x)

# calculate average depth
mean(x[,3])
# std dev
sqrt(var(x[,3]))

# mark areas that have a coverage below 20 in red
plot(x[,2], x[,3], col = ifelse(x[,3] < 20,'red','black'), pch=19, xlab='postion', ylab='coverage')

# to save a plot
png('mappings/covNODE20.png', width = 1200, height = 500)
plot(x[,2], x[,3], col = ifelse(x[,3] < 20,'red','black'), pch=19, xlab='postion', ylab='coverage')
dev.off()

The result plot will be looking similar to the one in Fig. 5.2

../_images/covNODE20.png

Fig. 5.2 A example coverage plot for a contig with highlighted in red regions with a coverage below 20 reads.

Todo

Look at the created plot. Explain why it makes sense that you find relatively bad coverage at the beginning and the end of the contig.

5.10.2. Stats with QualiMap

For a more in depth analysis of the mappings, one can use QualiMap.

QualiMap examines sequencing alignment data in SAM/BAM files according to the features of the mapped reads and provides an overall view of the data that helps to the detect biases in the sequencing and/or mapping of the data and eases decision-making for further analysis.

Installation:

conda install qualimap

Run QualiMap with:

qualimap bamqc -bam mappings/evolved-6.sorted.dedup.bam

This will create a report in the mapping folder. See this webpage to get help on the sections in the report.

Todo

Install QualiMap and investigate the mapping of the evolved sample. Write down your observations.

5.11. Sub-selecting reads

It is important to remember that the mapping commands we used above, without additional parameters to sub-select specific alignments (e.g. for Bowtie2 there are options like --no-mixed, which suppresses unpaired alignments for paired reads or --no-discordant, which suppresses discordant alignments for paired reads, etc.), are going to output all reads, including unmapped reads, multi-mapping reads, unpaired reads, discordant read pairs, etc. in one file. We can sub-select from the output reads we want to analyse further using SAMtools.

Todo

Explain what concordant and discordant read pairs are? Look at the Bowtie2 manual.

5.11.1. Concordant reads

We select read-pair that have been mapped in a correct manner (same chromosome/contig, correct orientation to each other).

samtools view -h -b -f 2 mappings/evolved-6.sorted.dedup.bam > mappings/evolved-6.sorted.concordant.bam
  • -h: Include the sam header
  • -b: Output will be bam-format
  • -f 2: Only extract correctly paired reads. -f extracts alignments with the specified SAM flag set.

5.11.2. Quality-based sub-selection

In this section we want to sub-select reads based on the quality of the mapping. It seems a reasonable idea to only keep good mapping reads. As the SAM-format contains at column 5 the \(MAPQ\) value, which we established earlier is the "MAPping Quality" in Phred-scaled, this seems easily achieved. The formula to calculate the \(MAPQ\) value is: \(MAPQ=-10*log10(p)\), where \(p\) is the probability that the read is mapped wrongly. However, there is a problem! While the MAPQ information would be very helpful indeed, the way that various tools implement this value differs. A good overview can be found here. Bottom-line is that we need to be aware that different tools use this value in different ways and the it is good to know the information that is encoded in the value. Once you dig deeper into the mechanics of the \(MAPQ\) implementation it becomes clear that this is not an easy topic. If you want to know more about the \(MAPQ\) topic, please follow the link above.

For the sake of going forward, we will sub-select reads with at least medium quality as defined by Bowtie2:

samtools view -h -b -q 20 mappings/evolved-6.sorted.dedup.concordant.bam > mappings/evolved-6.sorted.concordant.q20.bam
  • -h: Include the sam header
  • -q 20: Only extract reads with mapping quality >= 20

5.11.3. Unmapped reads

We could decide to use Kraken like in section Taxonomic investigation to classify all unmapped sequence reads and identify the species they are coming from and test for contamination.

Lets see how we can get the unmapped portion of the reads from the bam-file:

samtools view -b -f 4 mappings/evolved-6.sorted.dedup.bam > mappings/evolved-6.sorted.unmapped.bam

# count them
samtools view -c mappings/evolved-6.sorted.unmapped.bam
  • -b: indicates that the output is BAM.
  • -f INT: only include reads with this SAM flag set. You can also use the command samtools flags to get an overview of the flags.
  • -c: count the reads

Lets extract the fastq sequence of the unmapped reads for read1 and read2.

bamToFastq -i mappings/evolved-6.sorted.unmapped.bam -fq mappings/evolved-6.sorted.unmapped.R1.fastq -fq2  mappings/evolved-6.sorted.unmapped.R2.fastq

References

[TRAPNELL2009]Trapnell C, Salzberg SL. How to map billions of short reads onto genomes. Nat Biotechnol. (2009) 27(5):455-7. doi: 10.1038/nbt0509-455.