added mk.image function
GabrieleNocchi opened this issue · 0 comments
GabrieleNocchi commented
As I noticed it was not included in this package, I re-added the mk.image function (basically the heatmap of the SNPs with the individuals ordered by h_index). Note that in this heatmap I don't plot all markers, but only those included in the introgress output and with a cline p-value < 0.05. If you want to plot ALL the markers in the heatmap, just remove the argument loci.touse from the below code:
########################## heatmap of identified clines
# Read in the data from files
sep <- "\t"
gen.data <- read.table(file = "eatt_admix.txt", header = FALSE, sep = sep)
loci.data <- read.table(file = "eatt_loci.txt", header = TRUE, sep = sep)
p1 <- read.table(file = "eatt_p1data.txt", header = FALSE, sep = sep)
p2 <- read.table(file = "eatt_p2data.txt", header = FALSE, sep = sep)
# Prep data for analysis
count.matrix <- introgress::prepare.data(
admix.gen = gen.data,
loci.data = loci.data,
parental1 = p1,
parental2 = p2,
pop.id = FALSE,
ind.id = FALSE
)
hi.index <-
introgress::est.h(introgress.data = count.matrix,
loci.data = loci.data,
fixed = FALSE)
my_loci <- read.table("./introgress_output/EATT_summaryData.txt", h = F, skip =1)
my_loci$V4 <- as.numeric(as.character(my_loci$V4))
my_loci <- my_loci[my_loci$V4 < 0.05,]
index <- match(my_loci$V1, loci.data$locus)
gab_plot(introgress.data=count.matrix, loci.data=loci.data, marker.order=NULL, hi.index=hi.index, loci.touse = index, ylab.image="Individuals", main.image="", xlab.h="Population 2 ancestry", pdf=TRUE, out.file="clines_heatmap.pdf")
gab_plot is a modified version of mk.image, as mk.image has some bugs if you don't plot all the loci in the original data:
gab_plot <- function (introgress.data = NULL, loci.data = NULL, marker.order = NULL,
hi.index = NULL, ind.touse = NULL, loci.touse = NULL,
ylab.image = "Individuals", main.image = "", xlab.h = "population 2 ancestry",
col.image = NULL, pdf = TRUE, out.file = "image.pdf")
{
# Check if hi.index is a data frame, if yes extract the second column
if (is.data.frame(hi.index) == TRUE)
hi.index <- hi.index[, 2]
# Ensure hi.index is numeric
if (!is.numeric(hi.index)) {
stop("Please provide a numeric vector of hybrid indexes (hi.index).")
}
# Check for introgress.data and loci.data
if (is.null(introgress.data) == TRUE)
stop("error, introgress.data was not supplied")
if (is.null(loci.data) == TRUE)
stop("error, loci.data was not supplied")
# If introgress.data is a list, extract the second element
if (is.list(introgress.data) == TRUE)
count.matrix <- introgress.data[[2]]
else count.matrix <- introgress.data
# Set up PDF output if requested
if (pdf == TRUE) {
pdf(file = out.file, height = 9.5, width = 7)
}
# Handling linkage groups and marker count if loci.data has more than two columns
if (dim(loci.data)[2] > 2) {
lg <- sort(unique(as.numeric(loci.data[, 3])))
marker.n <- numeric(length(lg))
for (i in 1:length(lg)) {
marker.n[i] <- sum(loci.data[, 3] == lg[i])
}
}
# Default to all loci if loci.touse is not specified
if (is.null(loci.touse)) {
loci.touse <- 1:(dim(count.matrix)[1])
}
# Default to all individuals if ind.touse is not specified
if (is.null(ind.touse)) {
ind.touse <- 1:(dim(count.matrix)[2])
}
# Default color scale
if (is.null(col.image)) {
col.image <- c("#A1D99B", "#41AB5D", "#005A32")
}
# Handle marker order if provided
if (is.null(marker.order)) {
marker.order <- loci.touse
}
else if (is.character(marker.order) || is.factor(marker.order)) {
tmp <- numeric(length(marker.order))
for (i in 1:length(marker.order)) {
tmp[i] <- which(loci.data == as.character(marker.order)[i])
}
marker.order <- tmp
}
# Set layout
nf <- layout(matrix(c(1, 2), 1, 2, byrow = TRUE), c(6, 1), c(1, 1), respect = FALSE)
par(mar = c(5, 4, 1, 1))
# Fix for non-sequential loci: use seq_along(loci.touse) for consistent column sizes
image(x = seq_along(loci.touse), y = ind.touse,
z = count.matrix[marker.order, order(hi.index[ind.touse])],
col = col.image, breaks = c(-0.1, 0.9, 1.9, 2.1),
xlab = "", ylab = ylab.image, main = main.image, axes = FALSE)
# Add markers label
mtext("Markers", 1, line = 4)
# Y-axis for individuals
axis(2)
# X-axis with marker labels and gray vertical lines separating linkage groups
if (dim(loci.data)[2] == 2) {
axis(1, at = seq(0.5, length(loci.touse) + 0.5, 1), labels = FALSE)
axis(1, at = seq_along(loci.touse), tick = FALSE,
labels = loci.data[, 1][marker.order], cex.axis = 0.3, las = 2, line = 1)
} else {
axis(1, at = seq_along(loci.touse), tick = FALSE,
labels = loci.data[, 1][marker.order], cex.axis = 0.3, las = 2, line = 1)
axis(1, at = c(0, cumsum(marker.n) + 0.5), labels = FALSE)
abline(v = c(0, cumsum(marker.n) + 0.5), col = "lightgray", lwd = 0.6)
axis(1, at = diff(c(0, cumsum(marker.n)))/2 +
c(0, cumsum(marker.n)[-length(marker.n)] + 0.5),
tick = FALSE, labels = lg, cex.axis = 0.5, line = 0)
}
# Draw box around the heatmap
box()
# Side plot for hybrid index
par(mar = c(5, 0, 1, 1))
n.inds <- length(ind.touse)
plot(sort(hi.index[ind.touse]), 1:n.inds, axes = FALSE, xlab = "", ylab = "",
ylim = c(1 + 0.0355 * n.inds, n.inds - 0.0355 * n.inds),
cex = 0.6, type = "n", xlim = 0:1)
mtext("Proportion", 1, line = 3, cex = 0.6)
mtext(xlab.h, 1, line = 4, cex = 0.6)
abline(v = 0.5, col = "lightgray")
lines(sort(hi.index[ind.touse]), 1:n.inds)
axis(1, at = c(0, 0.5, 1), cex.axis = 0.6)
box()
# Close PDF device if writing to file
if (pdf == TRUE) {
dev.off()
}
}