search

xiaoming-liu / stairway-plot-v2

The stairway plot is a method for inferring detailed population demographic history using the site frequency spectrum (SFS) from DNA sequence data.
Other
31 stars 4 forks source link

How can I plot the demographic history for several populations in one figure? #2

Open lonemen opened 3 years ago

lonemen commented 3 years ago

Dear Dr.Liu Thank you for providing this great software. I wonder How I could plot the demographic history for several populations in one figure? Thank you!

lucasrocmoreira commented 3 years ago

I am also interested in knowing the answer to this question? Is there an easy implementation of plotting multiple populations/species, like in PSMC for example?

Thanks, Lucas

cistarsa commented 3 years ago

did anyone figure this out?

lucasrocmoreira commented 3 years ago

I believe I had to manually plot them in R, it was not too hard. Here is my code, I hope it helps.

Alaska <- read.table("Alaska.final.summary",header=T)
Eastern <- read.table("Eastern.final.summary",header=T)
Northwest <- read.table("Northwest.final.summary",header=T)
Rockies <- read.table("Rockies.final.summary",header=T)

pdf("PP-stairway2.plot.pdf",height = 6.25,width = 10)

plot(Alaska$year/1000,Alaska$Ne_median/1000, log=c("xy"), type="n", xlab="Time (1k years ago)", ylab="Effective Population Size (1k individuals)",xlim=c(4,1000),ylim=c(20,30000))

lines(Alaska$year/1000,Alaska$Ne_median/1000,type="s",col="orange",lwd = 5)
lines(Alaska$year/1000,Alaska$Ne_2.5./1000,type="s",col="orange",lty=3)
lines(Alaska$year/1000,Alaska$Ne_97.5./1000,type="s",col="orange",lty=3)

lines(Eastern$year/1000,Eastern$Ne_median/1000,type="s",col="green",lwd = 5)
lines(Eastern$year/1000,Eastern$Ne_2.5./1000,type="s",col="green",lty=3)
lines(Eastern$year/1000,Eastern$Ne_97.5./1000,type="s",col="green",lty=3)

lines(Northwest$year/1000,Northwest$Ne_median/1000,type="s",col="red",lwd = 5)
lines(Northwest$year/1000,Northwest$Ne_2.5./1000,type="s",col="red",lty=3)
lines(Northwest$year/1000,Northwest$Ne_97.5./1000,type="s",col="red",lty=3)

lines(Rockies$year/1000,Rockies$Ne_median/1000,type="s",col="blue",lwd = 5)
lines(Rockies$year/1000,Rockies$Ne_2.5./1000,type="s",col="blue",lty=3)
lines(Rockies$year/1000,Rockies$Ne_97.5./1000,type="s",col="blue",lty=3)

legend(400,28000,legend = c("Alaska","Eastern","Northwest","Rockies","95% CI"),col=c("orange","green","red","blue","black"),lty=c(1,1,1,1,3),lwd=c(5,5,5,5,1),cex=0.8)

dev.off()
cistarsa commented 3 years ago

thank you! very helpful code.

chenqing-1996 commented 3 years ago

I believe I had to manually plot them in R, it was not too hard. Here is my code, I hope it helps.

Alaska <- read.table("Alaska.final.summary",header=T)
Eastern <- read.table("Eastern.final.summary",header=T)
Northwest <- read.table("Northwest.final.summary",header=T)
Rockies <- read.table("Rockies.final.summary",header=T)

pdf("PP-stairway2.plot.pdf",height = 6.25,width = 10)

plot(Alaska$year/1000,Alaska$Ne_median/1000, log=c("xy"), type="n", xlab="Time (1k years ago)", ylab="Effective Population Size (1k individuals)",xlim=c(4,1000),ylim=c(20,30000))

lines(Alaska$year/1000,Alaska$Ne_median/1000,type="s",col="orange",lwd = 5)
lines(Alaska$year/1000,Alaska$Ne_2.5./1000,type="s",col="orange",lty=3)
lines(Alaska$year/1000,Alaska$Ne_97.5./1000,type="s",col="orange",lty=3)

lines(Eastern$year/1000,Eastern$Ne_median/1000,type="s",col="green",lwd = 5)
lines(Eastern$year/1000,Eastern$Ne_2.5./1000,type="s",col="green",lty=3)
lines(Eastern$year/1000,Eastern$Ne_97.5./1000,type="s",col="green",lty=3)

lines(Northwest$year/1000,Northwest$Ne_median/1000,type="s",col="red",lwd = 5)
lines(Northwest$year/1000,Northwest$Ne_2.5./1000,type="s",col="red",lty=3)
lines(Northwest$year/1000,Northwest$Ne_97.5./1000,type="s",col="red",lty=3)

lines(Rockies$year/1000,Rockies$Ne_median/1000,type="s",col="blue",lwd = 5)
lines(Rockies$year/1000,Rockies$Ne_2.5./1000,type="s",col="blue",lty=3)
lines(Rockies$year/1000,Rockies$Ne_97.5./1000,type="s",col="blue",lty=3)

legend(400,28000,legend = c("Alaska","Eastern","Northwest","Rockies","95% CI"),col=c("orange","green","red","blue","black"),lty=c(1,1,1,1,3),lwd=c(5,5,5,5,1),cex=0.8)

dev.off()

Hi lucasrocmoreira， I‘m confused about your code of plot, why do you use " log=c("xy") "? I find it wrong not to use " log=c("xy") " when plot, but I don;t what the reason is for it? Looking forward to your reply.

Best wishes.

lucasrocmoreira commented 3 years ago

The log=c("xy") argument just puts both axes on the log scale. You can choose to transform just the y or just the x axis, as well. I hope that makes sense.

chenqing-1996 commented 3 years ago

Hi lucasrocmoreira， Thank you very much! I can understand what you mean. Thanks again.

Best wishes.