Question: R: Manhattan Plot Of Fst Values Instead Of -Log(P)
3
gravatar for Abdel
7.3 years ago by
Abdel340
Amsterdam
Abdel340 wrote:

Hi all,

I have a tab delimited file with the following columns:

CHR  SNP  bp  Fst

where Fst is the fixation index. Does anyone know how to make a Manhattan plot of these Fst values in R? Many thanks!

R visualization fst • 10k views
ADD COMMENTlink modified 4.5 years ago by CB10 • written 7.3 years ago by Abdel340
6
gravatar for John
7.3 years ago by
John1.5k
John1.5k wrote:

I found the following blog useful, you need tweak the codes to fit F-value not -log10(p) value -

http://gettinggeneticsdone.blogspot.com/2011/04/annotated-manhattan-plots-and-qq-plots.html

Edits:

If you would like to use lattice or ggplot, we can create similar plot:

# just an example 
Fst <- rnorm(10000, 10, 5)
chr <- c(rep(1, 2500), rep(2, 2500), rep(3, 2500), rep(4, 2500))
BP <- c(1:2500, 1:2500, 1:2500, 1:2500)
mydf <- data.frame (Fst, chr, BP)
 require(ggplot2)
   qplot(BP, Fst, facets = . ~ chr, col = factor(chr), data = mydf)

I like bw theme better

qplot(BP, Fst, facets = . ~ chr, col = factor(chr), data = mydf) + theme_bw()

You need to reshape plot area or change arrangement of graph to make all facets in single line as Manhattan plot. The graph is pretty close but not exactly as in above blog code.

ADD COMMENTlink modified 7.3 years ago • written 7.3 years ago by John1.5k

You just have to replace bp by BP and Fst by P and it should work.

ADD REPLYlink written 7.3 years ago by Maxime Lamontagne2.1k
3
gravatar for Zev.Kronenberg
7.3 years ago by
United States
Zev.Kronenberg11k wrote:

I have done the exact plot you want for FST.

here is where [R] becomes awesome:

dat<-read.csv("your.data.txt", header=FALSE, sep="\t")
plot(dat$FST, col=as.factor(dat$CHR))

If you need the true coordinates ie BP I will send you a script I wrote to take that into account. Otherwise the code above will suffice.

Pretty much R gives the X values a number 1:length(dat$FST) and then colors the points by the dat$CHR

ADD COMMENTlink written 7.3 years ago by Zev.Kronenberg11k

Hi, I tried this code for data in the same format as above, am I missing something or should these two R commands be sufficient for plotting? Error in plot.window(...) : need finite 'xlim' values. Thanks

ADD REPLYlink modified 6.0 years ago • written 6.0 years ago by Rubal7760

please post "head(dat)"

ADD REPLYlink written 6.0 years ago by Zev.Kronenberg11k

Hi all,

I need to do the same plot as described above but I am facing the same problem.

 plot(dat$FST, col=as.factor(dat$CHR))
Error in plot.window(...) : need finite 'xlim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

 

what I have to do? what do you mean by "head(dat)"?

Thanks in advance

ADD REPLYlink written 4.8 years ago by Tohamy70

Would you mind posting the script that you wrote to take  BP into account? Cheers! :)

ADD REPLYlink written 5.1 years ago by learning2speakpirate.R!!! 0

https://github.com/jewmanchue/zk_tools/blob/master/CHR_POS_GENOMIC_POS

ADD REPLYlink written 4.8 years ago by Zev.Kronenberg11k

First of all thanks for your help.

I just need to use this code only

plot(dat$FST, col=as.factor(dat$CHR))

but it give me the prevoius error.

I dont need the true coordinates.

 

ADD REPLYlink written 4.8 years ago by Tohamy70
0
gravatar for CB
4.5 years ago by
CB10
US
CB10 wrote:

I want to create a Manhattan plot  to plot Fst values of each SNP. I used John's first code to plot Fst values for 3 different chromosomes and it worked well. But the size of chromosomes are fixed, and I want to set the maximum size for each chromosome. Is it possible? The code that I used was:

Fst <- rnorm(48167,0.13, 0.2)

chr <- c(rep(1, 14552), rep(2, 33582), rep(10, 33))

BP <- c(1:14552, 1:33582, 1:33)

data <- data.frame (Fst, chr, BP) qplot(BP, Fst, facets = . ~ chr, col = factor(chr), data = data)

 

I tried to use the package qqman to plot Fst values instead of P values and did not work.

Or I need to use the perl script that Zev.Kronenberg mentioned? Convert chrmosomes positions into genomic coordiantes?

ADD COMMENTlink modified 4.5 years ago • written 4.5 years ago by CB10

With ggplot2 you can facet on the chromosome :  something like +facet_grid(.~CHR).  You will have to play with the margins on the facets, but it does work.  If you use my script it make the process easier as you don't have to facet or fiddle with the coordinates. 

 

 

ADD REPLYlink modified 4.5 years ago • written 4.5 years ago by Zev.Kronenberg11k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1502 users visited in the last hour