gbradburd/conStruct

Covariance figure blank

Closed this issue · 7 comments

ruthrivkin commented

Hi Gideon,

I'm running in to an issue with making the diagnostic figures that I'm hoping you can help with. When I run conStruct for any value of K and make.figs = TRUE, it creates a blank covariance figure. The pdf file is there, and the plot axes, but there is nothing in the plot. The trace figures are all created with no issues. When I try to make the figures after running conStruct (make.all.the.plots), it errors out saying that I'm using too much memory. That is very strange, because it gives me that error message even when I run it on the hpc, which has plenty of memory allocation for both computing and temp storage.

Any insight you have would be appreciated! I have included an example of the blank plot.

Thanks,
Ruth

pb.spatial.k2_model.fit.CIs.chain_1.pdf

gbradburd commented

Hi Ruth!

Hmm that's weird. How many samples do you have? And, can you show the R code for how you're running your analysis and making this plot?

ruthrivkin commented

I have 411 samples, with ~3,600 SNPs. I imputed missing SNPs (which were less than 1% anyway) just in case missingness was playing a factor in the missing plots. I am able to see all the trace plots and create the admixture and pie charts without any issues. It's just the covariance plot I can't seem to make. Since its not entirely clear what model is best explains patterns in the data (comparing K = 1-7), I was hoping looking at the covariance curves would help identify meaningful layers.

My conStruct code is as follows:

`#Load genetic data
library(conStruct)
load("ld.myConStructData.RData")

#Load coordinates
coord = read.delim("Subset_Coord.txt", header = FALSE)
names(coord) <- c("Lon", "Lat")
coordinates <- as.matrix(coord)

#Geographic distance
library(geosphere)
dist <- distm(coordinates)

my.run.k2 <- conStruct(spatial = TRUE,
K = 2,
freqs = freqs,
geoDist = dist,
coords = coordinates,
prefix = "pb.spatial.k2",
make.figs = TRUE)`

Then to make the plots, I'm running:

`#Check figs
load("pb.spatial.k2_data.block.Robj")
load("pb.spatial.k2_conStruct.results.Robj")

make.all.the.plots(conStruct.results = conStruct.results,
data.block = data.block,
prefix = "k2",
layer.colors = NULL)`

gbradburd commented

Ok I think I know where the issue is! I suspect it's breaking in the gtools::combinations call inside plot.model.fit.CIs. I think I can do the same operation in a different and more efficient way. Can you send me your data.block and conStruct.results Robjs as a zipped file, either via github or over email? While I'm testing possible solutions, it would also be helpful for you to try adjusting the stack limit in the shell (outside of R) before you try to make the plot, as detailed here.

ruthrivkin commented
gbradburd commented

Hi Ruth,

Thanks for sending that. I tried to run it on my laptop (2020 MB Air w/ M1 chip and 8GB ram, soft stack size set to 8176 kbytes) and, although it was slow, I was surprised to see that it worked fine! Nevertheless, I wanted to put together a version of this function that avoids this issue.

plot.model.fit.CIs <- function(data.block,conStruct.results){
	cov.range <- range(c(data.block$obsCov,
						conStruct.results$posterior$par.cov))
	ut <- upper.tri(data.block$geoDist,diag=TRUE)
	graphics::plot(data.block$geoDist[ut],data.block$obsCov[ut],
    	xlab = "geographic distance", 
        ylab = "covariance",
        main="Cov/geoDist",
        ylim = cov.range, type = "n")
	CIs <- apply(conStruct.results$posterior$par.cov,2:3,function(x){quantile(x,c(0.025,0.975))})
	invisible(
		lapply(1:data.block$N,
			function(i){
				lapply(i:data.block$N,
					function(j){
						graphics::segments(x0 = data.block$geoDist[i,j],
											y0 = CIs[1,i,j],
											x1 = data.block$geoDist[i,j],
											y1 = CIs[2,i,j],
											col = grDevices::adjustcolor(1,0.1),
											lwd=1.5)
					})
			}))
	graphics::points(data.block$geoDist[ut],data.block$obsCov[ut],col=2,pch=20,cex=0.8)
	graphics::legend(x="topright",legend=c("observed","95% CI"),pch=c(19,NA),lty=c(NA,1),col=c(2,"gray"))
	return(invisible("plotted"))
}

If you load your data.block and conStruct.results .Robj files, attach the conStruct library, then define the function above so it's in R's working memory, you should be able to test it using:

plot.model.fit.CIs(data.block,conStruct.results[[1]])

Please give that a try and let me know how it goes!

ruthrivkin commented
gbradburd commented

ok great! I'll add it to the next release. I'm going to close this comment but please reopen if other issues crop up!