added mk.image function

[GabrieleNocchi]

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()
    }
}