# Finding LOD peaks

## Overview

Teaching:30 min

Exercises:30 minQuestions

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().