Finding LOD peaks

Overview

Teaching: 30 min
Exercises: 30 min

Questions

• How do I locate LOD peaks above a certain threshold value?

Objectives

• Locate LOD peaks above a threshold value throughout the genome.

• Identify the Bayesian credible interval for a QTL peak.

Once we have LOD scores from a genome scan and a significance threshold, we can look for QTL peaks associated with the phenotype. High LOD scores indicate the neighborhood of a QTL but don’t give its precise position. To find the exact position of a QTL, we define an interval that is likely to hold the QTL.

We’ll use the Bayesian credible interval, which is a method for defining QTL intervals. It describes the probability that the interval contains the true value. Credible intervals make a probabilistic statement about the true value, for example, a 95% credible interval states that there is a 95% chance that the true value lies within the interval.

To find peaks above a given threshold LOD value, use the function find_peaks(). It can also provide Bayesian credible intervals by using the argument prob (the nominal coverage for the Bayes credible intervals). Set the argument expand2markers = FALSE to keep from expanding the interval out to typed markers, or exclude this argument if you’d like to include flanking markers.

You need to provide both the scan1() output, the marker/pseudomarker map and a threshold. We will use the 95% threshold from the permutations in the previous lesson.

operm <- scan1perm(pr, iron$pheno, Xcovar=Xcovar, n_perm=1000)
thr <- summary(operm)
find_peaks(scan1_output = out, map = map, threshold = thr, prob = 0.95, expand2markers = FALSE)

  lodindex lodcolumn chr  pos       lod ci_lo ci_hi
1        1     liver   2 56.8  4.957564  54.3  70.3
2        1     liver   7 50.1  4.050766  17.1  53.6
3        1     liver   8 40.0  3.802511  34.0  55.0
4        1     liver  16 28.6  7.681569  21.6  32.6
5        2    spleen   8 13.6  4.302919   5.0  23.0
6        2    spleen   9 56.6 12.063226  54.6  58.6

The find_peaks() function can also pick out multiple peaks on a chromosome: each peak must exceed the chosen threshold, and the argument peakdrop indicates the amount that the LOD curve must drop between the lowest of two adjacent peaks. Use this feature with caution.

find_peaks(scan1_output = out, map = map, threshold = thr, peakdrop = 1.8, prob = 0.95, expand2markers = FALSE)

  lodindex lodcolumn chr  pos       lod ci_lo ci_hi
1        1     liver   2 56.8  4.957564  54.3  70.3
2        1     liver   7 25.1  4.040021  15.1  27.1
3        1     liver   7 50.1  4.050766  41.1  53.6
4        1     liver   8 40.0  3.802511  34.0  55.0
5        1     liver  16 28.6  7.681569  21.6  32.6
6        2    spleen   8 13.6  4.302919   5.0  23.0
7        2    spleen   9 56.6 12.063226  54.6  58.6

Each row shows a different peak; the columns show the peak location, LOD score and the lower and upper interval endpoints.

Challenge 1

Find peaks in the genome scan object called out that meet a threshold of 3 and are in the interval described by a 2 point LOD drop from the peak. How many peaks meet the LOD threshold of 3 and lie within the interval defined by a 2 point LOD drop from the maximum peaks on each chromosome?

Solution to Challenge 1

find_peaks(out, map, threshold=3, drop=2) produces 7 peaks on 6 different chromosomes that meet a LOD threshold of 3 and are within the interval defined by a 2-LOD drop from the maximum peak on each chromosome.

Challenge 2

1). Calculate the 90% Bayes credible interval. For chromosome 16 in liver, what is the range of the interval that has a 90% chance of containing the true QTL position?
2). Compare with the 95% Bayes credible interval calculated earlier. How does the interval change as you increase the probability? Why?

Solution to Challenge 2

1). find_peaks(out, map, prob=0.90, expand2markers = FALSE) produces a range from 25.1 to 40.4.
2). find_peaks(out, map, prob=0.95, expand2markers = FALSE) produces a range from 6.6 to 40.4, which is much broader than that of a 90% probability. The interval widens because the probability that the interval contains the true QTL position has increased.

Key Points

• LOD peaks and support intervals can be identified with find_peaks().