GN020540 is the strongest single population difference in 13.4: the only gene with a significant gene-body difference (Qui vs Rio), and covered by seven 1 kb DMRs from the promoter to the last exon, all lower in Rio. This notebook looks at the locus directly:

per-oyster methylation and read depth along the locus;
what the gene is (annotation) and whether it is expressed (ONT RNA-seq, 09);
SNPs from 13.7 in the region;
the two haplotype assemblies (HapA, HapB) at the locus;
methylation at the other genomic copies of the repeat inside the intron, to tell a locus-specific (cis) change from a repeat-family-wide one.

Upstream: 13.2 (strand-merged CpG calls), 13.4, 13.5 (Swiss-Prot hits), 13.7 (SNPs). The haplotype assemblies are downloaded to data/ if absent (as in 12-compare-genomes.py). RNA-seq BAMs and counts are read from gannet, where bu.sh backs up output/09-align-v1-clean/.

Samples: the 9 independent oysters in data/13-wgbs-samples.tsv (duplicated libraries excluded; see README Known issues).

Where to run: a compute node, inside the lab R container (~10 min, 16 cores for minimap2). Command-line tools are called by absolute path from the srlab conda envs, which run inside the container; paths are relative to code/.

Network. klone compute nodes have no internet access. The RNA-seq inputs (count matrix, per-library read counts from the remote BAMs) and the haplotype assemblies are fetched only if not already cached in out_dir / data/. Run code/15-GN020540-fetch.R once from a login node (lab R container) to fill the cache before knitting on a compute node.

library(tidyverse)

out_dir  <- "../output/15-GN020540-locus"
meth_dir <- "../output/13-wgbs/methylation"
genome   <- "../data/13-wgbs-genome/merged_out.fasta"
gff_file <- "../data/13-wgbs-genome/GN.gene.gff3"
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

gene_id  <- "GN020540"
chrom    <- "Chromosome_9B"
view     <- c(7180000, 7260000)   # wide window for depth
locus    <- c(7215000, 7234000)   # gene +/- ~6 kb
promoter_bp <- 2000

envs     <- "/mmfs1/gscratch/srlab/sr320/miniforge3/envs"
minimap2 <- file.path(envs, "myflow/bin/minimap2")
samtools <- file.path(envs, "wgbs-snp/bin/samtools")
bcftools <- file.path(envs, "wgbs-snp/bin/bcftools")
bedtools <- file.path(envs, "wgbs-snp/bin/bedtools")
threads  <- 16

gannet <- "https://gannet.fish.washington.edu/v1_web/owlshell/bu-github"
Sys.setenv(MINIMAP2 = minimap2, SAMTOOLS = samtools, BEDTOOLS = bedtools,
           OUT = out_dir, METH = meth_dir, GENOME = genome, CHROM = chrom,
           VLO = view[1], VHI = view[2], THREADS = threads)

Gene model and annotation

gff <- read_tsv(gff_file, comment = "#", col_names = FALSE, show_col_types = FALSE)
gene_rows <- gff %>% filter(str_detect(X9, paste0("(ID|Parent)=", gene_id, "(;|\\.|$)")))
gene  <- gene_rows %>% filter(X3 == "gene") %>% transmute(start = X4, end = X5, strand = X7)
exons <- gene_rows %>% filter(X3 == "exon") %>% transmute(start = X4, end = X5)
gene; exons

## # A tibble: 1 x 3
##     start     end strand
##     <dbl>   <dbl> <chr> 
## 1 7221497 7227996 +

## # A tibble: 4 x 2
##     start     end
##     <dbl>   <dbl>
## 1 7221497 7221664
## 2 7221806 7221874
## 3 7227212 7227376
## 4 7227892 7227996

# Nearest other genes
gff %>% filter(X1 == chrom, X3 == "gene", X5 > gene$start - 50000, X4 < gene$end + 50000) %>%
  transmute(gene = str_extract(X9, "ID=[^;]+"), start = X4, end = X5)

## # A tibble: 1 x 3
##   gene          start     end
##   <chr>         <dbl>   <dbl>
## 1 ID=GN020540 7221497 7227996

annot <- read_tsv("../output/13-wgbs/13.5-gene-annotation/GN_gene_annotation.tsv",
                  show_col_types = FALSE)
annot %>% filter(gene_id == !!gene_id) %>%
  select(gene_id, entry_name, protein_name, pident, aln_len, qlen, slen, evalue)

## # A tibble: 1 x 8
##   gene_id  entry_name  protein_name          pident aln_len  qlen  slen   evalue
##   <chr>    <chr>       <chr>                  <dbl>   <dbl> <dbl> <dbl>    <dbl>
## 1 GN020540 NU107_MOUSE Nuclear pore complex~   54.5     167   168   926 2.36e-49

The predicted protein (qlen) is compared with the full-length Swiss-Prot hit (slen): a much shorter query means a partial gene model or a gene fragment.

Methylation and depth per oyster

samples <- read_tsv("../data/13-wgbs-samples.tsv", show_col_types = FALSE) %>%
  select(sample, population)
writeLines(samples$sample, file.path(out_dir, "samples.txt"))

Per-oyster CpG calls in the wide window, and each oyster’s genome-wide mean depth per covered CpG (to express local depth relative to the genome).

while read -r s; do
  (
    f="${METH}/${s}.CpG_report.merged_CpG_evidence.cov.gz"
    zcat "$f" | awk -F'\t' -v OFS='\t' -v c="$CHROM" -v lo="$VLO" -v hi="$VHI" -v s="$s" \
      '$1 == c && $2 >= lo && $2 <= hi {print s, $2, $5, $6}' > "${OUT}/cpg_${s}.tsv"
    zcat "$f" | awk -F'\t' -v s="$s" '{n++; d += $5 + $6} END {print s "\t" d / n}' \
      > "${OUT}/depth_${s}.tsv"
  ) &
done < "${OUT}/samples.txt"
wait
cat "${OUT}"/cpg_*.tsv  > "${OUT}/cpg_window.tsv"
cat "${OUT}"/depth_*.tsv > "${OUT}/genome_depth.tsv"
wc -l "${OUT}/cpg_window.tsv"

## cat: ../output/15-GN020540-locus/cpg_window.tsv: input file is output file
## 9777 ../output/15-GN020540-locus/cpg_window.tsv

cpg <- read_tsv(file.path(out_dir, "cpg_window.tsv"), col_names = c("sample", "pos", "nC", "nT"),
                show_col_types = FALSE) %>%
  left_join(samples, by = "sample") %>%
  mutate(cov = nC + nT)
gdepth <- read_tsv(file.path(out_dir, "genome_depth.tsv"), col_names = c("sample", "genome_depth"),
                   show_col_types = FALSE)

region_of <- function(pos) case_when(
  pos < gene$start - promoter_bp ~ "upstream flank",
  pos < gene$start ~ "promoter (2 kb)",
  pos <= gene$end ~ "gene body",
  TRUE ~ "downstream flank")

by_region <- cpg %>%
  filter(between(pos, locus[1], locus[2])) %>%
  mutate(region = factor(region_of(pos), levels = c("upstream flank", "promoter (2 kb)",
                                                    "gene body", "downstream flank"))) %>%
  group_by(population, sample, region) %>%
  summarise(cpgs = n(), meth_pct = 100 * sum(nC) / sum(cov), mean_depth = mean(cov),
            .groups = "drop") %>%
  left_join(gdepth, by = "sample") %>%
  mutate(rel_depth = mean_depth / genome_depth)
write_tsv(by_region, file.path(out_dir, "methylation_by_region.tsv"))

by_region %>% select(sample, region, meth_pct) %>%
  pivot_wider(names_from = region, values_from = meth_pct) %>%
  mutate(across(where(is.numeric), ~ round(.x, 1)))

## # A tibble: 9 x 5
##   sample `upstream flank` `promoter (2 kb)` `gene body` `downstream flank`
##   <chr>             <dbl>             <dbl>       <dbl>              <dbl>
## 1 Pum1               66.7              76.9        36.9                9.1
## 2 Qui1                0                92.1        61.1               15.6
## 3 Qui2               50                45.8        37.1               16.8
## 4 Qui3               50                68.1        35.7                6.7
## 5 Qui4               50                45          40.5               12.1
## 6 Rio1                0                23.4         2.7                0.8
## 7 Rio2                0                14.4         2.3                1.3
## 8 Rio3              100                13.3         4.3                1.5
## 9 Rio4               60                11.2         3.4                3.6

Per-CpG states in the gene body: are oysters near 0 / 100% (allelic or fixed) or intermediate?

cpg %>%
  filter(between(pos, gene$start, gene$end), cov >= 3) %>%
  mutate(state = cut(100 * nC / cov, c(-1, 20, 80, 100), labels = c("<=20%", "20-80%", ">80%"))) %>%
  count(population, sample, state) %>%
  group_by(sample) %>% mutate(frac = round(n / sum(n), 2)) %>% select(-n) %>%
  pivot_wider(names_from = state, values_from = frac, values_fill = 0)

## # A tibble: 9 x 5
## # Groups:   sample [9]
##   population sample `<=20%` `20-80%` `>80%`
##   <chr>      <chr>    <dbl>    <dbl>  <dbl>
## 1 Pum        Pum1      0.19     0.77   0.04
## 2 Qui        Qui1      0        0.59   0.41
## 3 Qui        Qui2      0.17     0.76   0.07
## 4 Qui        Qui3      0.27     0.63   0.1 
## 5 Qui        Qui4      0.24     0.69   0.07
## 6 Rio        Rio1      0.96     0.04   0   
## 7 Rio        Rio2      0.95     0.05   0   
## 8 Rio        Rio3      0.95     0.05   0   
## 9 Rio        Rio4      0.94     0.06   0

bins <- cpg %>%
  filter(between(pos, locus[1], locus[2])) %>%
  mutate(bin = floor(pos / 500) * 500 + 250) %>%
  group_by(population, sample, bin) %>%
  summarise(meth = 100 * sum(nC) / sum(cov), cov = sum(cov), .groups = "drop") %>%
  filter(cov >= 5)
dmr <- read_tsv("../output/13-wgbs/13.4-methylation/DMR_Qui_vs_Rio.tsv", show_col_types = FALSE) %>%
  filter(seqnames == chrom, end >= locus[1], start <= locus[2])

p <- ggplot(bins, aes(bin / 1e6, meth, colour = population, group = sample)) +
  geom_rect(data = dmr, inherit.aes = FALSE,
            aes(xmin = start / 1e6, xmax = end / 1e6, ymin = -Inf, ymax = Inf),
            fill = "grey90") +
  geom_rect(data = exons, inherit.aes = FALSE,
            aes(xmin = start / 1e6, xmax = end / 1e6, ymin = -8, ymax = -3), fill = "black") +
  annotate("segment", x = gene$start / 1e6, xend = gene$end / 1e6, y = -5.5, yend = -5.5) +
  geom_line(alpha = 0.8) + geom_point(size = 0.8) +
  scale_colour_manual(values = c(Pum = "grey40", Qui = "#1b9e77", Rio = "#d95f02")) +
  labs(x = paste(chrom, "(Mb)"), y = "CpG methylation (%), 500 bp bins",
       title = paste(gene_id, "methylation per oyster"),
       subtitle = "Grey = 13.4 DMR tiles; black = exons") +
  theme_bw()
p

ggsave(file.path(out_dir, paste0(gene_id, "_methylation.png")), p, width = 9, height = 4.5, dpi = 150)

Read depth relative to each oyster’s genome-wide mean, in 4 kb bins over the wide window. Depth well above 1 in every oyster points to sequence present more than once in the genome but once in the assembly; a Rio-only drop would point to a deletion.

depth_bins <- cpg %>%
  mutate(bin_kb = floor(pos / 4000) * 4) %>%
  group_by(sample, bin_kb) %>%
  summarise(depth = mean(cov), .groups = "drop") %>%
  left_join(gdepth, by = "sample") %>%
  mutate(rel = round(depth / genome_depth, 2)) %>%
  select(sample, bin_kb, rel) %>%
  pivot_wider(names_from = bin_kb, values_from = rel)
write_tsv(depth_bins, file.path(out_dir, "relative_depth_4kb.tsv"))
depth_bins

## # A tibble: 9 x 21
##   sample `7180` `7184` `7188` `7192` `7196` `7200` `7204` `7208` `7212` `7216`
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
## 1 Pum1     1.03   1.34   1.8    0.35   0.61   0.41   0.17   0.35   0.23   0.52
## 2 Qui1     1.3    1.4    0.89   0.89   0.46   0.64   0.18   0.35   0.24   0.97
## 3 Qui2     0.91   1.21   1.5    0.82   0.18   0.18   0.37   0.27   0.73   0.73
## 4 Qui3     1.49   1.26   1.66   0.58   0.37   0.41   0.17   0.28   0.33   1.08
## 5 Qui4     0.94   1.38   0.92   1.16   0.58   0.58   0.19   0.58   0.19   0.32
## 6 Rio1     1.16   1.75   1.03   0.47   0.4    0.37   0.47   0.23   0.39   0.27
## 7 Rio2     1.24   1.23   1.24   0.36   0.33   0.46   0.18   0.18   0.18   1.37
## 8 Rio3     1.04   1.27   1.34   1.02   0.3    0.71   0.41  NA      0.2    0.51
## 9 Rio4     1.34   1.49   1.27   0.57   0.42   0.54   0.16   0.24   0.41   1.55
## # i 10 more variables: `7220` <dbl>, `7224` <dbl>, `7228` <dbl>, `7232` <dbl>,
## #   `7236` <dbl>, `7240` <dbl>, `7244` <dbl>, `7248` <dbl>, `7252` <dbl>,
## #   `7256` <dbl>

SNPs in the region (13.7)

vcf_q20 <- "../output/13-wgbs/13.7-snps/wgbs_snps.q20.vcf.gz"
vcf_all <- "../output/13-wgbs/13.7-snps/wgbs_snps.all.vcf.gz"
region  <- sprintf("%s:%d-%d", chrom, locus[1], locus[2])
q20 <- system2(bcftools, c("view", "-H", "-r", region, vcf_q20), stdout = TRUE)
all_var <- system2(bcftools, c("query", "-r", region, "-i", "'ALT!=\".\"'", "-f", "'%QUAL\\n'", vcf_all),
                   stdout = TRUE)
tibble(region = region, snps_qual20 = length(q20),
       candidate_variant_sites = length(all_var),
       max_qual = suppressWarnings(max(as.numeric(all_var))))

## # A tibble: 1 x 4
##   region                        snps_qual20 candidate_variant_sites max_qual
##   <chr>                               <int>                   <int>    <dbl>
## 1 Chromosome_9B:7215000-7234000           0                     911       18

Expression (ONT RNA-seq, 09)

StringTie counts from 09 (current pipeline), and direct read counts at the gene from each library’s BAM, read remotely through its index. The RNA-seq oysters are different individuals from the WGBS oysters.

counts_file <- file.path(out_dir, "09_gene_count_matrix.csv")
if (!file.exists(counts_file))
  download.file(file.path(gannet, "project-ostrea-chil/output/09-align-v1-clean/stringtie_quant/gene_count_matrix.csv"),
                counts_file, quiet = TRUE)
counts <- read_csv(counts_file, show_col_types = FALSE)
libsize <- colSums(select(counts, -gene_id))
counts %>% filter(str_detect(gene_id, gene_id)) %>%
  pivot_longer(-gene_id, names_to = "library", values_to = "count") %>%
  mutate(cpm = round(1e6 * count / libsize[library], 1))

## # A tibble: 496,008 x 4
##    gene_id             library         count   cpm
##    <chr>               <chr>           <dbl> <dbl>
##  1 GN013935|GN013935.1 Pum_barcode02       0     0
##  2 GN013935|GN013935.1 Pum_barcode12_1     0     0
##  3 GN013935|GN013935.1 Pum_barcode13       0     0
##  4 GN013935|GN013935.1 Pum_barcode13_1     0     0
##  5 GN013935|GN013935.1 Qui_barcode01       0     0
##  6 GN013935|GN013935.1 Qui_barcode01_1     0     0
##  7 GN013935|GN013935.1 Qui_barcode02_1     0     0
##  8 GN013935|GN013935.1 Qui_barcode12       0     0
##  9 GN013935|GN013935.1 Rio_barcode03       0     0
## 10 GN013935|GN013935.1 Rio_barcode03_1     0     0
## # i 495,998 more rows

# How common is "zero in half the libraries" among expressed genes?
m <- as.matrix(select(counts, -gene_id))
expressed <- rowSums(m) >= 100
c(expressed_genes = sum(expressed), zero_in_6_plus = sum(rowSums(m[expressed, ] == 0) >= 6))

## expressed_genes  zero_in_6_plus 
##           10850            3022

libs <- setdiff(names(counts), "gene_id")
bam_url <- function(l) file.path(gannet, "project-ostrea-chil/output/09-align-v1-clean",
                                 paste0(l, "_combined_clean_porechop.sorted.bam"))
gene_region <- sprintf("%s:%d-%d", chrom, gene$start, gene$end)
reads_file <- file.path(out_dir, "rnaseq_reads_at_gene.tsv")
if (file.exists(reads_file)) {
  reads <- read_tsv(reads_file, show_col_types = FALSE)
} else {
old <- setwd(out_dir)   # samtools caches the remote .bai files in the working directory
reads <- map_dfr(libs, function(l) {
  idx <- system2(samtools, c("idxstats", bam_url(l)), stdout = TRUE)
  mapped <- sum(as.numeric(map_chr(strsplit(idx, "\t"), 3)))
  n <- as.integer(system2(samtools, c("view", "-c", "-F", "0x904", "-q", "20",
                                      bam_url(l), gene_region), stdout = TRUE))
  tibble(library = l, population = substr(l, 1, 3), mapped = mapped, gene_reads = n,
         per_million = round(1e6 * n / mapped, 1))
})
setwd(old)
write_tsv(reads, reads_file)
}
reads

## # A tibble: 12 x 5
##    library         population mapped gene_reads per_million
##    <chr>           <chr>       <dbl>      <dbl>       <dbl>
##  1 Pum_barcode02   Pum        368127          3         8.1
##  2 Pum_barcode12_1 Pum        333678          4        12  
##  3 Pum_barcode13   Pum        205659          0         0  
##  4 Pum_barcode13_1 Pum        204083          4        19.6
##  5 Qui_barcode01   Qui        265974          2         7.5
##  6 Qui_barcode01_1 Qui        197503          2        10.1
##  7 Qui_barcode02_1 Qui        398045          1         2.5
##  8 Qui_barcode12   Qui        360704          1         2.8
##  9 Rio_barcode03   Rio        349603          1         2.9
## 10 Rio_barcode03_1 Rio        333746          2         6  
## 11 Rio_barcode14   Rio        336004          5        14.9
## 12 Rio_barcode14_1 Rio        340277          2         5.9

reads %>% group_by(population) %>% summarise(mean_per_million = mean(per_million))

## # A tibble: 3 x 2
##   population mean_per_million
##   <chr>                 <dbl>
## 1 Pum                    9.93
## 2 Qui                    5.72
## 3 Rio                    7.42

Haplotype assemblies

The merged reference’s Chromosome_9B is HapB’s Chromosome_9B. The 50 kb around the gene, and separately the high-depth intron segment, are aligned to HapA, HapB and the merged assembly.

hap_urls <- c(HapA = file.path(gannet, "OCEAN/docs/jbrowse/data/HapA/Och_HapA_assembly.fa"),
              HapB = file.path(gannet, "OCEAN/docs/jbrowse/data/HapB/Och_HapB_assembly.fa"))
hap_files <- c(HapA = "../data/Och_HapA_assembly.fa", HapB = "../data/Och_HapB_assembly.fa")
for (h in names(hap_files)) {
  if (!file.exists(hap_files[h])) {
    download.file(hap_urls[h], paste0(hap_files[h], ".part"), mode = "wb", quiet = TRUE)
    file.rename(paste0(hap_files[h], ".part"), hap_files[h])
  }
}
file.size(hap_files) / 1e6

## [1]  802.7907 1150.8297

"$SAMTOOLS" faidx "$GENOME" "${CHROM}:7200000-7250000" > "${OUT}/locus50k.fa"
"$SAMTOOLS" faidx "$GENOME" "${CHROM}:7222000-7227000" > "${OUT}/intron_segment.fa"
for t in HapA:../data/Och_HapA_assembly.fa HapB:../data/Och_HapB_assembly.fa merged:"$GENOME"; do
  n=${t%%:*}; f=${t#*:}
  [[ -s "${OUT}/locus_vs_${n}.paf" ]] || \
    "$MINIMAP2" -t "$THREADS" -x asm20 -c --secondary=yes -N 200 -p 0.1 "$f" "${OUT}/locus50k.fa" \
      > "${OUT}/locus_vs_${n}.paf" 2> /dev/null
  [[ -s "${OUT}/intron_vs_${n}.paf" ]] || \
    "$MINIMAP2" -t "$THREADS" -x asm20 -c --secondary=yes -N 500 -p 0.05 "$f" "${OUT}/intron_segment.fa" \
      > "${OUT}/intron_vs_${n}.paf" 2> /dev/null
done
wc -l "${OUT}"/*.paf

##    360 ../output/15-GN020540-locus/intron_vs_HapA.paf
##    501 ../output/15-GN020540-locus/intron_vs_HapB.paf
##    501 ../output/15-GN020540-locus/intron_vs_merged.paf
##     68 ../output/15-GN020540-locus/locus_vs_HapA.paf
##     17 ../output/15-GN020540-locus/locus_vs_HapB.paf
##     17 ../output/15-GN020540-locus/locus_vs_merged.paf
##   1464 total

read_paf <- function(f) {
  read_tsv(f, col_names = FALSE, show_col_types = FALSE, col_select = 1:13) %>%
    transmute(qstart = X3, qend = X4, strand = X5, target = X6, tstart = X8, tend = X9,
              identity = round(X10 / X11, 3), aln_len = X11, mapq = X12,
              type = str_remove(X13, "tp:A:"))
}
offset <- 7200000   # query coordinates are relative to this; gene at q 21,497-27,996

# Primary placement of the 50 kb locus in each haplotype
map_dfr(c("HapA", "HapB"), ~ read_paf(file.path(out_dir, paste0("locus_vs_", .x, ".paf"))) %>%
          filter(type == "P") %>% mutate(assembly = .x)) %>%
  mutate(q = paste0(qstart + offset, "-", qend + offset)) %>%
  select(assembly, q, strand, target, tstart, tend, identity, aln_len)

## # A tibble: 0 x 8
## # i 8 variables: assembly <chr>, q <chr>, strand <chr>, target <chr>,
## #   tstart <dbl>, tend <dbl>, identity <dbl>, aln_len <dbl>

# Self-similar blocks upstream of the gene in HapB (tandem repeat array)
read_paf(file.path(out_dir, "locus_vs_HapB.paf")) %>%
  filter(type == "S", target == chrom, aln_len >= 5000) %>%
  arrange(qstart) %>% select(qstart, qend, tstart, tend, identity, aln_len)

## # A tibble: 0 x 6
## # i 6 variables: qstart <dbl>, qend <dbl>, tstart <dbl>, tend <dbl>,
## #   identity <dbl>, aln_len <dbl>

# Copies of the intron segment (>= 2 kb aligned)
map_dfr(c("HapA", "HapB", "merged"), ~ read_paf(file.path(out_dir, paste0("intron_vs_", .x, ".paf"))) %>%
          filter(aln_len >= 2000) %>%
          summarise(copies = n(), copies_95pct = sum(identity >= 0.95)) %>% mutate(assembly = .x))

## # A tibble: 3 x 3
##   copies copies_95pct assembly
##    <int>        <int> <chr>   
## 1    223           78 HapA    
## 2    368          125 HapB    
## 3    369          132 merged

Other copies of the intronic repeat

If the Rio loss were a property of the repeat family (e.g. genome-wide repeat methylation differing by population), the other copies would also be lower in Rio.

awk -F'\t' -v OFS='\t' -v c="$CHROM" \
  '$11 >= 2000 && !($6 == c && $8 < 7230000 && $9 > 7215000) {print $6, $8, $9}' \
  "${OUT}/intron_vs_merged.paf" | sort -k1,1 -k2,2n | "$BEDTOOLS" merge -i - > "${OUT}/repeat_copies.bed"
cut -f1,2 "${GENOME}.fai" | sort -k1,1 > "${OUT}/genome.sorted.sizes"
while read -r s; do
  (
    zcat "${METH}/${s}.CpG_report.merged_CpG_evidence.cov.gz" \
      | awk -F'\t' -v OFS='\t' '{print $1, $2 - 1, $2, $5, $6}' \
      | "$BEDTOOLS" intersect -sorted -g "${OUT}/genome.sorted.sizes" -a - -b "${OUT}/repeat_copies.bed" -wa -wb \
      | awk -F'\t' -v OFS='\t' -v s="$s" '{print s, $6 ":" $7, $4, $5}' > "${OUT}/repeat_${s}.tsv"
  ) &
done < "${OUT}/samples.txt"
wait
wc -l < "${OUT}/repeat_copies.bed"

## 368

rep <- map_dfr(samples$sample, ~ read_tsv(file.path(out_dir, paste0("repeat_", .x, ".tsv")),
                                          col_names = c("sample", "copy", "nC", "nT"),
                                          show_col_types = FALSE)) %>%
  left_join(samples, by = "sample")

per_copy <- rep %>% group_by(population, sample, copy) %>%
  summarise(calls = sum(nC + nT), meth = 100 * sum(nC) / calls, .groups = "drop") %>%
  filter(calls >= 20)

per_copy %>% group_by(population, sample) %>%
  summarise(copies = n(), mean_per_copy = round(mean(meth), 1),
            copies_below_10pct = sum(meth < 10), .groups = "drop")

## # A tibble: 9 x 5
##   population sample copies mean_per_copy copies_below_10pct
##   <chr>      <chr>   <int>         <dbl>              <int>
## 1 Pum        Pum1      321          18.9                216
## 2 Qui        Qui1      301          20.6                202
## 3 Qui        Qui2      307          19.9                217
## 4 Qui        Qui3      305          20.2                210
## 5 Qui        Qui4      297          20.5                204
## 6 Rio        Rio1      320          19.8                219
## 7 Rio        Rio2      309          18.8                214
## 8 Rio        Rio3      294          21.4                198
## 9 Rio        Rio4      305          20.9                202

# The GN020540 copy itself, for comparison (intron segment only)
cpg %>% filter(between(pos, 7222000, 7227000)) %>%
  group_by(population, sample) %>%
  summarise(meth_pct = round(100 * sum(nC) / sum(cov), 1), .groups = "drop")

## # A tibble: 9 x 3
##   population sample meth_pct
##   <chr>      <chr>     <dbl>
## 1 Pum        Pum1       33.3
## 2 Qui        Qui1       53.8
## 3 Qui        Qui2       30.9
## 4 Qui        Qui3       28.2
## 5 Qui        Qui4       31.7
## 6 Rio        Rio1        2.9
## 7 Rio        Rio2        2  
## 8 Rio        Rio3        4.8
## 9 Rio        Rio4        3.1

Summary

Rio oysters are near-unmethylated across GN020540 (promoter to last exon); Quilhua and Pum1 are partly methylated, with no overlap between the groups.
The gene model is a fragment (qlen vs slen above), and the intron holds a ~4 kb interspersed repeat with hundreds of genomic copies, which explains the doubled read depth and the intermediate methylation there.
No confident SNP lies in the region; ONT expression is too sparse to compare.
HapA carries an identical copy of the gene; the haplotypes differ by a tandem array just upstream, which short bisulfite reads cannot genotype.
Other copies of the intronic repeat are methylated the same in every population, so the Rio loss is specific to this locus.

sessionInfo()

## R version 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
## [1] C
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.9.4 forcats_1.0.0   stringr_1.5.1   dplyr_1.1.4    
##  [5] purrr_1.0.4     readr_2.1.5     tidyr_1.3.1     tibble_3.2.1   
##  [9] ggplot2_3.5.2   tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] bit_4.6.0         gtable_0.3.6      jsonlite_2.0.0    crayon_1.5.3     
##  [5] compiler_4.4.3    tidyselect_1.2.1  parallel_4.4.3    jquerylib_0.1.4  
##  [9] textshaping_1.0.0 systemfonts_1.2.2 scales_1.3.0      yaml_2.3.10      
## [13] fastmap_1.2.0     R6_2.6.1          labeling_0.4.3    generics_0.1.3   
## [17] knitr_1.50        munsell_0.5.1     bslib_0.9.0       pillar_1.10.2    
## [21] tzdb_0.5.0        rlang_1.1.5       utf8_1.2.4        cachem_1.1.0     
## [25] stringi_1.8.7     xfun_0.52         sass_0.4.9        bit64_4.6.0-1    
## [29] timechange_0.3.0  cli_3.6.4         withr_3.0.2       magrittr_2.0.3   
## [33] digest_0.6.37     grid_4.4.3        vroom_1.6.5       hms_1.1.3        
## [37] lifecycle_1.0.4   vctrs_0.6.5       evaluate_1.0.3    glue_1.8.0       
## [41] farver_2.1.2      ragg_1.3.3        colorspace_2.1-1  rmarkdown_2.29   
## [45] tools_4.4.3       pkgconfig_2.0.3   htmltools_0.5.8.1

15 - GN020540: a Rio-specific loss of methylation

2026-10-10