Plotting Linkage Disequilibrium Decay Across The Genome
2
5
Entering edit mode
8.1 years ago
Javier2013 ▴ 130

Hi everyone:

I am trying to plot the linkage disequilibrium decay across the genome using genome-wide snps and get a plot that looks something like this:

http://postimg.org/image/5bkdju6jx/

In the methods sections the authors specified that they used the --r2 command in plink. I read the details of this command in the plink website, but I am still not sure how should I use the output file to plot a similar thing.

Any help would be really appreciated.

Thank you in advance.

Javier

Ps: Heres the article's title: Population Genetic Structure of the People of Qatar (2010)

plink ld snps plot • 17k views
ADD COMMENT
0
Entering edit mode

Can you post the first 4 lines of the plink.ld file containing the pai wise correlations?

ADD REPLY
0
Entering edit mode

Hi Irsan:

Heres the output file I get from plink using default values

http://postimg.org/image/87hq8nodp/

I got an answer from another post, and what I need to do is calculate the distance between SNP A and SNP B and calculate the r2 averages for every distance binned for example every 1kb. Do you know how to do a script that could do that?

Thank you!

ADD REPLY
0
Entering edit mode

Hi

i have two file:ped and map

which command of plink is requir to build input file for R?(for linkage disequilibrium decay across the genome)

ADD REPLY
0
Entering edit mode

Your question is unrelated to this old thread and should be asked in a separate question.

ADD REPLY
3
Entering edit mode
8.1 years ago
Irsan ★ 7.4k

In the plink.ld file you have colums CHR_A, BP_A, SNP_A, CHR_B, BP_B, SNP_B and R2.

in R do

# install and load ggplot2
install.packages("ggplot2")
library(ggplot2)

# import the data
ld <- read.table("plink.ld", sep="\t",header=T)

# plot the average correlation for each snp distance
ggplot(ld) +
geom_line(aes(x=BP_B - BP_A, y = R2))
ADD COMMENT
0
Entering edit mode

If you don't like scripting you can also do this with Excel pivot tables

ADD REPLY
0
Entering edit mode

Hi Irsan,

Thank you for the reply.

I just tried plotting my .ld file but I couldn't get any clear result. I just get a plot where everything is black http://s24.postimg.org/9mhyytpqt/Screen_Shot_2013_10_27_at_3_43_34_PM.png

So I did some scripting on my own and found out that the r2 average is too high even at 70kb apart. http://s12.postimg.org/812gxfza5/Screen_Shot_2013_10_27_at_3_49_53_PM.png

I usually do not use the ggplot package so I am not very sure why am I getting that plot...

But I am also thinking that theres something wrong with the my output file. Maybe something with the commands? Have you done a LD decay plot before??

Here is my log file:

Options in effect: --noweb --bfile testLWK --r2 --ld-window-kb 70 --ld-window-r2 0 --ld-window 10 --out LWK_70kb_10snps

Thanks again.

ADD REPLY
0
Entering edit mode

No i have never used plink before, for me this is just a visualization question. And it should be really straightforward. Check if your data frame is alright by dim(ld), str(ld), summary(ld)

ADD REPLY
0
Entering edit mode

In R version 3.3.2 (2016-10-31) the line

ggplot(ld) + geom_line(aes(x=BP_B - BP_A, y = R2))

gives

Error in eval(expr, envir, enclos) : object 'BP_B' not found

ADD REPLY
1
Entering edit mode
3.8 years ago
hfang4 ▴ 10

I had the same problem as described by veronicaschroeder78:

Error in eval(expr, envir, enclos) : object 'BP_B' not found.

I solved the problem by modifying a little be in Iran's codes: change sep="\t" to sep="", others keep the same, and also change "geom_line" to "geom_point" will give you scatter plot of the LD.

ADD COMMENT
0
Entering edit mode

Please suggest how to add non-linear regression tendency line to the LD plot.

ADD REPLY

Login before adding your answer.

Traffic: 1479 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6