---
title: "09-Ptua-Gene-Methylation"
output: html_document
---

# Gene Methylation Analysis in Ptua


## gene file

```{bash}
head ../data/Pocillopora_meandrina_HIv1.genes.gff3

```




```{bash}
cut -f3 ../data/Pocillopora_meandrina_HIv1.genes-validated.gff3 | \
grep -v "^#" | \
sort | uniq -c | sort -nr
```

```{bash}
gawk '$3=="gene" {
    split($9, a, ";");
    for(i in a){
        if(a[i] ~ /^ID=/){
            sub("ID=", "", a[i]);
            gene_id = a[i];
        }
    }
    print $1 "\t" $4-1 "\t" $5 "\t" gene_id;
}' ../data/Pocillopora_meandrina_HIv1.genes-validated.gff3 > ../data/Pocillopora_meandrina_HIv1.GENESonly-validated.bed
```

```{bash}
head ../data/Pocillopora_meandrina_HIv1.GENESonly-validated.bed
wc -l ../data/Pocillopora_meandrina_HIv1.GENESonly-validated.bed
```


```{bash}
mkdir -p ../data/10xbedgraph
cd ../data/10xbedgraph
wget -r -l1 --no-parent -nd -A '*10x.bedgraph' \
     'https://gannet.fish.washington.edu/metacarcinus/E5/Ptuahiniensis/20250821_meth_Ptua/'
```

```{bash}
head ../data/10xbedgraph/POC-40-TP1_10x.bedgraph
```

```{bash}
for f in ../data/10xbedgraph/*_10x.bedgraph; do
    sort -k1,1 -k2,2n "$f" > ${f%.bedgraph}.sorted.bedgraph
done
```



```{bash}
head ../data/10xbedgraph/POC-40-TP1_10x.sorted.bedgraph
```


```{bash}
mkdir -p ../output/09-Ptua-Gene-Methylation/

sort -k1,1 -k2,2n ../data/Pocillopora_meandrina_HIv1.GENESonly-validated.bed \
  > ../data/Pocillopora_meandrina_HIv1.GENESonly-validated.sorted.bed

for f in ../data/10xbedgraph/*_10x.sorted.bedgraph; do
    sample=$(basename "$f" _10x.sorted.bedgraph)

    bedtools map \
        -a ../data/Pocillopora_meandrina_HIv1.GENESonly-validated.sorted.bed \
        -b "$f" \
        -c 4 \
        -o mean,count \
        > ../output/09-Ptua-Gene-Methylation/gene_methylation_${sample}.tsv
done
```


```{bash}
head ../output/09-Ptua-Gene-Methylation/gene_methylation_POC-40-TP*.tsv
```
Lets filter so all column 6 (count) is greater than 10 - meaning at least 10 CpG sites were covered in that gene.

```{bash}
for f in ../output/09-Ptua-Gene-Methylation/gene_methylation*.tsv; do
    out="../output/09-Ptua-Gene-Methylation/$(basename "${f%.tsv}_gt10.tsv")"
    awk '$6 > 10' "$f" > "$out"
done
```


```{bash}
head ../output/09-Ptua-Gene-Methylation/gene_methylation_POC-40-TP*_gt10.tsv
```


```{r}
library(tidyverse)

# Directory containing all filtered tsv files
indir <- "../output/09-Ptua-Gene-Methylation"

# List all filtered files
files <- list.files(indir, pattern = "gene_methylation.*_gt10.tsv", full.names = TRUE)

# Read and combine
df_list <- lapply(files, function(f) {
  sample_name <- basename(f) |> 
    sub("gene_methylation_", "", x = _) |> 
    sub("_gt10.tsv", "", x = _)

  read_tsv(f, col_names = FALSE) |>
    select(X4, X5) |>          # X4 = gene_id, X5 = methylation percentage
    rename(gene_id = X4, !!sample_name := X5)
})

# Full join all files by gene_id
merged <- reduce(df_list, full_join, by = "gene_id")

# Filter: keep only genes with no missing values across samples
merged_complete <- merged |> 
  filter(if_all(-gene_id, ~ !is.na(.x)))

# Save the filtered version
write_tsv(merged_complete, "../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_complete.tsv")
```

```{bash}
head ../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_complete.tsv
wc -l ../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_complete.tsv
```


```{r}
library(tidyverse)

# Directory containing all filtered tsv files
indir <- "../output/09-Ptua-Gene-Methylation"

# List all filtered files
files <- list.files(indir, pattern = "gene_methylation.*_gt10.tsv", full.names = TRUE)

# Read and combine
df_list <- lapply(files, function(f) {
  sample_name <- basename(f) |> 
    sub("gene_methylation_", "", x = _) |> 
    sub("_gt10.tsv", "", x = _)

  read_tsv(f, col_names = FALSE) |>
    select(X4, X5) |>
    rename(gene_id = X4, !!sample_name := X5)
})

# Full join all files by gene_id
merged <- reduce(df_list, full_join, by = "gene_id")

# ---- NEW PART: keep genes present in >=75% of samples ----
merged_75 <- merged %>% 
  mutate(
    # fraction of samples with non-NA values (ignores gene_id column)
    frac_nonNA = rowMeans(!is.na(across(-gene_id)))
  ) %>%
  filter(frac_nonNA >= 0.75) %>% 
  select(-frac_nonNA)

# Save 75% filtered matrix
write_tsv(merged_75, "../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_75pct.tsv")
```

```{bash}
wc -l ../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_75pct.tsv
head ../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_75pct.tsv
```






```{r}
library(tidyverse)
library(pheatmap)

# ------------------------------------------------------------
# Load data
# ------------------------------------------------------------
df <- read_tsv("../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_75pct.tsv")

mat <- df %>%
  column_to_rownames("gene_id") %>%
  as.matrix()

# ------------------------------------------------------------
# Step 1 — Impute NA with row means
# ------------------------------------------------------------
row_means <- rowMeans(mat, na.rm = TRUE)

for (i in seq_len(nrow(mat))) {
  if (is.na(row_means[i])) row_means[i] <- 0  # all-NA rows
  mat[i, is.na(mat[i, ])] <- row_means[i]
}

# ------------------------------------------------------------
# Step 2 — Remove zero-variance rows
# ------------------------------------------------------------
keep <- apply(mat, 1, function(x) sd(x) > 0)
mat2 <- mat[keep, , drop = FALSE]

# ------------------------------------------------------------
# Step 3 — Remove rows still containing Inf/NaN
# ------------------------------------------------------------
mat2 <- mat2[apply(mat2, 1, function(x) all(is.finite(x))), , drop = FALSE]

# ------------------------------------------------------------
# Heatmap 1 — Scaled (same as before)
# ------------------------------------------------------------
pheatmap(
  mat2,
  scale = "row",
  clustering_method = "ward.D2",
  show_rownames = FALSE,
  main = "Gene-level Methylation Patterns (Scaled)"
)

# ------------------------------------------------------------
# Heatmap 2 — Actual methylation values (Unscaled)
# ------------------------------------------------------------
pheatmap(
  mat2,
  scale = "none",
  clustering_method = "ward.D2",
  show_rownames = FALSE,
  main = "Gene-level Methylation Values (Raw)",
  color = colorRampPalette(c("navy","green", "firebrick"))(100)  # optional nicer palette
)
```





```{r}
library(tidyverse)
library(RColorBrewer)

# -------------------------
# Load data
# -------------------------
df <- read_tsv("../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_75pct.tsv")

# -------------------------
# Convert to matrix
# -------------------------
mat <- df %>%
  column_to_rownames("gene_id") %>%
  as.matrix()

# -------------------------
# Step 1: Impute NA with row means
# -------------------------
row_means <- rowMeans(mat, na.rm = TRUE)

for (i in seq_len(nrow(mat))) {
  if (is.na(row_means[i])) row_means[i] <- 0  # if entire row is NA
  mat[i, is.na(mat[i, ])] <- row_means[i]
}

# -------------------------
# Step 2: Remove zero-variance genes 
# (these blow up PCA scaling)
# -------------------------
keep <- apply(mat, 1, function(x) sd(x) > 0)
mat2 <- mat[keep, ]

# -------------------------
# Step 3: Remove genes containing any Inf/NaN
# -------------------------
mat2 <- mat2[apply(mat2, 1, function(x) all(is.finite(x))), ]

# -------------------------
# Step 4: Transpose (samples as rows)
# -------------------------
mat2_t <- t(mat2)

# -------------------------
# Extract sample group
# -------------------------
sample_groups <- str_extract(rownames(mat2_t), "POC-[0-9]+")

ngroups <- length(unique(sample_groups))

colors <- brewer.pal(max(3, ngroups), "Dark2")[1:ngroups]
group_colors <- setNames(colors, unique(sample_groups))

# -------------------------
# Step 5: Run PCA safely
# -------------------------
pca <- prcomp(mat2_t, scale. = TRUE)

# -------------------------
# Step 6: Build plotting DF
# -------------------------
pca_df <- data.frame(
  PC1 = pca$x[,1],
  PC2 = pca$x[,2],
  sample = rownames(pca$x),
  group = sample_groups
)

# -------------------------
# Step 7: Plot PCA
# -------------------------
ggplot(pca_df, aes(PC1, PC2, color = group, label = sample)) +
  geom_point(size = 4) +
  geom_text(vjust = -0.7, size = 3) +
  scale_color_manual(values = group_colors) +
  theme_minimal(base_size = 14) +
  labs(
    title = "PCA of Gene-Level Methylation",
    color = "Sample Group"
  )
```

```{r}
library(tidyverse)
library(pheatmap)

# ------------------------------------------------
# Load and convert to matrix
# ------------------------------------------------
mat <- read_tsv("../output/09-Ptua-Gene-Methylation/Ptua-gene-methylation_75pct.tsv") %>%
  column_to_rownames("gene_id") %>%
  as.matrix()

# ------------------------------------------------
# Step 1 — Impute NA with row means
# ------------------------------------------------
row_means <- rowMeans(mat, na.rm = TRUE)

for (i in seq_len(nrow(mat))) {
  if (is.na(row_means[i])) row_means[i] <- 0  # handles all-NA rows
  mat[i, is.na(mat[i, ])] <- row_means[i]
}

# ------------------------------------------------
# Step 2 — Remove zero-variance genes
# (these make correlation produce NA columns)
# ------------------------------------------------
keep <- apply(mat, 1, function(x) sd(x) > 0)
mat2 <- mat[keep, , drop = FALSE]

# ------------------------------------------------
# Step 3 — Remove any rows still containing Inf/NaN
# ------------------------------------------------
mat2 <- mat2[apply(mat2, 1, function(x) all(is.finite(x))), , drop = FALSE]

# ------------------------------------------------
# Step 4 — Compute correlation matrix safely
# ------------------------------------------------
cor_mat <- cor(mat2, use = "pairwise.complete.obs")

# ------------------------------------------------
# Step 5 — Remove NA rows/columns in correlation matrix (if any remain)
# ------------------------------------------------
finite_cols <- apply(cor_mat, 2, function(x) all(is.finite(x)))
cor_mat2 <- cor_mat[finite_cols, finite_cols]

# ------------------------------------------------
# Step 6 — Plot sample–sample correlation heatmap
# ------------------------------------------------
pheatmap(cor_mat2,
         main = "Sample–Sample Correlation",
         clustering_method = "average",
         fontsize = 12)
```


