R: Plotting MSMC/PSMC results using ggplot2
1
0
Entering edit mode
5.1 years ago
dthorbur ★ 2.5k

I have finished my coalescent simulation analysis for all my populations, and now want to plot the results using ggplot, but I don't know how to parameterise it so that it matches the expectations for such a plot type. The above plot is kind of the expectation, but when I use the same parameters in ggplot, I produce the bottom one.

It's the more "jagged" lines that trail off the plot that I don't know how to create. The labels and axis are easy to modify.

A sample of the code I used for the top is here:

plot(CA_L_fixed_rho2_8cores$left_time_boundary/mu*gen, (1/CA_L_fixed_rho2_8cores$lambda)/(2*mu), log="x",ylim=c(0,100000), type="n", xlab="Years Ago", ylab="Effective Population Size")

lines(CA_L_fixed_rho2_8cores$left_time_boundary/mu*gen, (1/CA_L_fixed_rho2_8cores$lambda)/(2*mu), type="s", col=as.character(Pop_Col_Matrix[grep("CA_L", Pop_Col_Matrix$Population),2]))

lines(CA_R_fixed_rho2_8cores$left_time_boundary/mu*gen, (1/CA_R_fixed_rho2_8cores$lambda)/(2*mu), type="s", col=as.character(Pop_Col_Matrix[grep("CA_R", Pop_Col_Matrix$Population),2]))

And the ggplot is here:

ggplot(output, aes(x = left_time_boundary/mu*gen, y = (1/lambda_00)/(2*mu), linetype = Ecotype, colour = Population)) +
geom_line() +
ylim(0,100000)+
scale_x_log10(limits = c(100,10000)) +
theme_classic()

EDIT** Added 1 population of example data.

time_index left_time_boundary right_time_boundary   lambda_00 Population
        0        0.00000e+00         1.21183e-06 1.58059e+00       CA_L
        1        1.21183e-06         2.45515e-06 7.37336e+01       CA_L
        2        2.45515e-06         3.73162e-06 6.25560e+01       CA_L
        3        3.73162e-06         5.04307e-06 2.34338e+02       CA_L
        4        5.04307e-06         6.39146e-06 7.09668e+02       CA_L
        5        6.39146e-06         7.77894e-06 1.30547e+03       CA_L
        6        7.77894e-06         9.20785e-06 1.78073e+03       CA_L
        7        9.20785e-06         1.06807e-05 2.02606e+03       CA_L
        8        1.06807e-05         1.22004e-05 2.07207e+03       CA_L
       9        1.22004e-05         1.37699e-05 2.06643e+03       CA_L
      10        1.37699e-05         1.53926e-05 2.57822e+03       CA_L
      11        1.53926e-05         1.70722e-05 2.57822e+03       CA_L
      12        1.70722e-05         1.88129e-05 2.77558e+03       CA_L
      13        1.88129e-05         2.06194e-05 2.77558e+03       CA_L
      14        2.06194e-05         2.24967e-05 2.68765e+03       CA_L
      15        2.24967e-05         2.44506e-05 2.68765e+03       CA_L
      16        2.44506e-05         2.64877e-05 2.68153e+03       CA_L
      17        2.64877e-05         2.86154e-05 2.68153e+03       CA_L
      18        2.86154e-05         3.08421e-05 2.54312e+03       CA_L
      19        3.08421e-05         3.31774e-05 2.54312e+03       CA_L
      20        3.31774e-05         3.56325e-05 2.43018e+03       CA_L
      21        3.56325e-05         3.82205e-05 2.43018e+03       CA_L
      22        3.82205e-05         4.09563e-05 2.29070e+03       CA_L
      23        4.09563e-05         4.38581e-05 2.29070e+03       CA_L
      24        4.38581e-05         4.69472e-05 2.14761e+03       CA_L
      25        4.69472e-05         5.02496e-05 2.14761e+03       CA_L
      26        5.02496e-05         5.37967e-05 2.00975e+03       CA_L
      27        5.37967e-05         5.76280e-05 2.00975e+03       CA_L
      28        5.76280e-05         6.17928e-05 1.87895e+03       CA_L
      29        6.17928e-05         6.63548e-05 1.87895e+03       CA_L
      30        6.63548e-05         7.13978e-05 1.74989e+03       CA_L
      31        7.13978e-05         7.70355e-05 1.74989e+03       CA_L
      32        7.70355e-05         8.34270e-05 1.62090e+03       CA_L
      33        8.34270e-05         9.08054e-05 1.62090e+03       CA_L
      34        9.08054e-05         9.95322e-05 1.48383e+03       CA_L
      35        9.95322e-05         1.10213e-04 1.48383e+03       CA_L
      36        1.10213e-04         1.23983e-04 1.33870e+03       CA_L
      37        1.23983e-04         1.43390e-04 1.33870e+03       CA_L
      38        1.43390e-04         1.76568e-04 1.11796e+03       CA_L
      39        1.76568e-04                 Inf 1.11796e+03       CA_L

MSMC-Plots

MSMC ggplot R • 5.5k views
ADD COMMENT
0
Entering edit mode

Those two plots look like they were generated from different data sets. So it's a little bit hard to figure out what you want exactly. But I think what you're after is essentially cropping so it's "zoomed" in. In that case us coord_cartesian(ylim=c(), xlim=c()) to set the upper and lower bounds.

https://stackoverflow.com/questions/25685185/limit-ggplot2-axes-without-removing-data-outside-limits-zoom

ADD REPLY
1
Entering edit mode

Hi Amar, it's the same data. I've added one population worth of data so you can see for yourself. It's not zooming that I'm after, it's to make the lines less "smoothed" and more abrupt when they go back in time as this is the norm for this type of plot in publications. Among other differences that are obvious from the plots.

ADD REPLY
0
Entering edit mode

Miles,

I'm attempting to do something similar but have two species for which I've generated psmc plots. I just want to plot them on the same axes to save some space in my figures and to allow direct comparison.

I don't typically use R but can if needed. My main question is this: how did you get the data structure you posted above?

time_index left_time_boundary right_time_boundary   lambda_00 Population
        0        0.00000e+00         1.21183e-06 1.58059e+00       CA_L
        1        1.21183e-06         2.45515e-06 7.37336e+01       CA_L

How did you get this information from the psmc file? Did you get this information from the psmc file? If not, from where did it come? If I can get my psmc data into a format like this, I could do what I need to but I don't know how to get that step completed.

Any help would be appreciated.

ADD REPLY
1
Entering edit mode

You're likely to get more help if you post your own question. Piggybacking on another one makes it tough for other users to find your question.

ADD REPLY
5
Entering edit mode
5.1 years ago
Haci ▴ 730

I think you are looking for geom_step(), which would generate a "staircase plot". Just replace geom_line() with geom_step().

ADD COMMENT
0
Entering edit mode

Thank you. I never realised this type of plot had a specific name.

ADD REPLY

Login before adding your answer.

Traffic: 1904 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