something like this?

  ggplot(df1, aes((start + stop) / 2, hetero)) +
  facet_grid(. ~ chromosome) +
  geom_rect(aes(
    xmin = min((start + stop) / 2),
    xmax = max((start + stop) / 2),
    ymin = -Inf,
    ymax = +Inf,
    fill = factor(col)
  )) +
  scale_fill_manual(values = c("grey90", "white")) +
  theme_void() +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none"
  ) +
  geom_point(aes(colour = factor(col))) +
  geom_point(data=df1[df1$hetero>0.4, ], aes((start + stop) / 2, hetero), colour="red")+
  scale_colour_manual(values = c("steelblue", "lightblue"))

input:

df1 = data.frame(
  chromosome = as.factor(rep(
    1:10, 100, each = 10, replace = T
  )),
  start = round(1:10000),
  stop =  round(20:10019),
  hetero = sample(c(
    runif(9000, 0.05, 0.2), runif(1000, 0.2, 0.6)
  ))
)
df1$col = ifelse(as.numeric(df1$chromosome) %% 2 == 0, 0, 1)
> head(df1)
  chromosome start stop     hetero col
1          1     1   20 0.07051028   1
2          1     2   21 0.12652906   1
3          1     3   22 0.08135082   1
4          1     4   23 0.19199763   1
5          1     5   24 0.09772686   1
6          1     6   25 0.38813218   1

Sort answer already provided by @MAPK: use col = myFactorColumn.

Longer answer, prepare data for plotting manhattan, then plot, see example, important bits:

• we don't have to use actual bp positions when we have many points, we just sort on chr and bp, and get row_number.
• we mark chromosome on x-axis at start of each minimum position per chromosome.

Here is the code, see the comments in the code for more info:

# example data
set.seed(1)
df1 <- data.frame(chromosome = sample(1:23, 10000, replace = TRUE),
                  start = round(sample(1:100000, 10000), 0),
                  end = 1,
                  heterozygosity = runif(10000, -0.01, 0.40),
                  window = 1)

library(dplyr)

# plot data prep
plotDat <- df1 %>% 
  # sort by chr and start, then get 1:n for x-axis values (bp)
  arrange(chromosome, start) %>% 
  mutate(bp = row_number()) %>% 
  # get the x value where we need to put chr name, in this case at min(bp)
  group_by(chromosome) %>% 
  mutate(lab = if_else(min(start) == start, chromosome, NA_integer_),
         # here we assign chr colour. "%%" give remainder, ie, is it even or odd.
         myCol = if_else(chromosome %% 2 == 1, "black", "grey"),
         # adding custom settings for plot, based on your ifelse conditions
         myPch = if_else(heterozygosity > 0, 21, 22),
         myBg = if_else(heterozygosity > 0.1374, "red", "white")
         ) %>% 
  ungroup() %>% 
  # convert back to dataframe (I prefer to work with dataframes, instead of tibbles)
  data.frame()

# plot without x axis
plot(heterozygosity ~ bp,
     col = myCol, pch = myPch, bg = myBg,
     data = plotDat,
     xaxt = "n", bty = "l", xlab = "chr")
# add custom axis labels
axis(1, at = plotDat[ !is.na(plotDat$lab), "bp"],
     labels = plotDat[ !is.na(plotDat$lab), "chromosome"],
     las = 2)

You could simply use col = your-factors option for your factor. What's your factor here? Try something like this: https://stackoverflow.com/questions/7721262/colouring-plot-by-factor-in-r