R_scripts.txt

### These are the R scripts used in the article "Genomic identification of Wolbachia reveals widespread signals of codivergence in the butterfly Eudaminae subfamily" 


##################### Drop Tip #####################
# Use according to the trees you want to obtain for mantel tests

library(ape)

setwd("C:/Users/pedro/OneDrive/Documents/PhD/Projects/Eudaminae_MacroWolbachia/PairwiseDistances_R")


# Import the tree file in Newick format
# Replace "your_treefile.treefile" with the name of your file
treeSuperA <- read.tree("wolbConcatenatedAB.newick")

# Visualize the tree to identify tip names
plot(treeSuperA)

# Print the tip names (to see which ones you might want to remove)
print(treeSuperA$tip.label)

#### This is for only ONE sample
# Remove one or more tips (replace "TipName" with the name of the tip to remove)
treeSuperA <- drop.tip(treeSuperA, "EMW233.B")

# Remove one or more tips (replace "TipName" with the name of the tip to remove)
treeSuperA <- drop.tip(treeSuperA, c("Ob_Wba", "wAlbB", "wAu", "wBm", "wBol1", "wCle", "wDacA", "wDacB", "wDi", "wFol", "wGmm", "wHa", "wMel", "wMelPop", "wNfe", "wNfla", "wNleu", "wNo", "wNpa", "wOneA1", "wOo", "wOv", "wPip", "wPip_Mol", "wRi", "wSpc", "wStri", "wSuzi", "wTpre", "wUni", "wVitA", "wVitB", "wWb", "EMW012.B", "EMW020.B", "EMW075.B", "EMW079.B", "EMW080.B", "EMW081.B", "EMW082.B", "EMW088.B", "EMW090.B", "EMW109.B", "EMW113.B", "EMW117.B", "EMW118.B", "EMW119.B", "EMW120.B", "EMW124.B", "EMW125.B", "EMW138.B", "EMW145.B", "EMW150.B", "EMW151.B", "EMW152.B", "EMW158.B", "EMW163.B", "EMW171.B", "EMW175.B", "EMW176.B", "EMW179.B", "EMW182.B", "EMW185.B", "EMW186.B", "EMW191.B", "EMW192.B", "EMW193.B", "EMW200.B", "EMW205.B", "EMW222.B", "EMW224.B", "EMW229.B", "EMW233.B", "EMW244.B", "EMW250.B", "EMW256.B", "EMW257.B", "EMW270.B", "EMW273.B , .B", "EMW169.B", "EMW110.B", "EMW170.B", "EMW149.B", "EMW142.B", "EMW043.B", "EMW260.B", "EMW252.B", "EMW235.B", "EMW246.B", "EMW207.B", "EMW217.B", "EMW246.A", "EMW203.A", "EMW237.A", "EMW065.A", "EMW169.A", "EMW076.A", "EMW249.A", "EMW251.A", "EMW279.A", "EMW100.A", "EMW067.A", "EMW192.B", "EMW193.B", "EMW192.A", "EMW193.A", "EMW205.A", "EMW205.B", "EMW229.B", "EMW206.A", "EMW233.A"))

# To remove multiple tips, use a vector of names
# Example: new_tree <- drop.tip(tree, c("Tip1", "Tip2"))

# Visualize the modified tree
plot(treeSuperA)

# Save the modified tree back to a Newick file
write.tree(treeSuperA, file = "PairwiseWolbSuperA.newick")

############### Cophenetic Phylo Part 

# Load the ape package if not already loaded
library(ape)

# Load the modified tree from the file "new_tree.tree"
modTreeSuperA <- read.tree("PairwiseWolbSuperA.newick")

# Compute the cophenetic distance matrix
cophDistSuperA <- cophenetic.phylo(modTreeSuperA)

# View the distance matrix
print(cophDistSuperA)

# If you want to save the distance matrix to a file (e.g., CSV)
write.csv(cophDistSuperA, file = "pairwiseDistWolbA.csv")


##################### Mantel tests, Heatmaps, and PACo #####################

install.packages("vegan")

library(vegan)

setwd("C:\\Users\\pedro\\OneDrive\\Documents\\PhD\\Projects\\Eudaminae_MacroWolbachia\\PairwiseDistances_R\\PariwiseAndHeatmap_feb2025")
# Load the two distance matrices
# Replace these with your actual matrices (Wolbachia and butterfly)
WolbMatrixA <- as.matrix(read.csv("WolbSuperAMatrix.csv", row.names = 1))
ButtMatrixA <- as.matrix(read.csv("ButtSuperAMatrix.csv", row.names = 1))

#Check rows and columns
rownames(WolbMatrixA)
colnames(WolbMatrixA)
rownames(ButtMatrixA)
colnames(ButtMatrixA)


# Transform/Ensure the row and column names of both matrices match
WolbMatrixA <- WolbMatrixA[rownames(ButtMatrixA), colnames(ButtMatrixA)]

#Check
all(rownames(WolbMatrixA) == rownames(ButtMatrixA))

# Perform the Mantel test
mantel_resultA <- mantel(WolbMatrixA, ButtMatrixA, method = "pearson", permutations = 9999)

# Print the results
print(mantel_resultA)


########### HEATMAP


install.packages("pheatmap")
library(pheatmap)
library(grid)

# Compute absolute differences between the two matrices
diff_matrixA <- abs(ButtMatrixA - WolbMatrixA)

# Plot heatmap of the difference matrix
pheatmap(diff_matrixA,
         main = "Heatmap of Absolute Differences Between Matrices",
         color = colorRampPalette(c("white", "blue"))(50),
         clustering_distance_rows = "euclidean",
         clustering_distance_cols = "euclidean",
         clustering_method = "average")

# Improved heatmap with better aesthetics for publication
heatmapA <- grid.grabExpr(pheatmap(diff_matrixA,
         main = "Heatmap of Absolute Differences - Supergroup B",
         color = colorRampPalette(c("white", "blue"))(100),  # Light Blue to Tomato Red,  # Diverging color palette
         clustering_distance_rows = "euclidean",  # Euclidean distance for rows
         clustering_distance_cols = "euclidean",  # Euclidean distance for columns
         clustering_method = "average",  # Average linkage for clustering
         fontsize = 13,  # Smaller font size for readability on smaller plots
         fontsize_row = 22,  # Smaller row label font size
         fontsize_col = 22,  # Smaller column label font size
         cellwidth = 20,  # Smaller cell width for better compactness
         cellheight = 20,  # Smaller cell height for better compactness
         border_color = "gray",  # Border color for cells
         legend = TRUE,  # Show the legend
         legend_title = "Difference",  # Legend title
         annotation_legend = TRUE,  # Add annotation to the legend if applicable
         annotation_names_row = TRUE,  # Show row names
         annotation_names_col = TRUE,  # Show column names
         angle_col = 90,  # Rotate column labels for better fit
         cluster_rows = TRUE,  # Cluster rows
         cluster_cols = TRUE,  # Cluster columns
         show_rownames = TRUE,  # Show row names
         show_colnames = TRUE,  # Show column names
         width = 6,  # Width of the saved plot (in inches)
         height = 6  # Height of the saved plot (in inches)
))

# Alternatively, save the heatmap to a file with custom size using ggsave
ggsave("SuperAHeatMap.png", plot = heatmapA, width = 15, height = 15, units = "in")
ggsave("SuperAHeatMap.svg", plot = heatmapA, width = 22, height = 22, units = "in")

# Add subtitle
grid.text("Phylogenetic Trees: Butterflies (Rows) | Wolbachia (Columns)", 
          x = 0.5, y = 0.92, gp = gpar(fontsize = 12, fontface = "italic"))


########### PACo

# Load required packages
install.packages(c("vegan", "ape"))
library(vegan)
library(ape)


setwd("C:\\Users\\pedro\\OneDrive\\Documents\\PhD\\Projects\\Eudaminae_MacroWolbachia\\PairwiseDistances_R\\PariwiseAndHeatmap_feb2025")
# import distance matrices from dropped tip trees
# Read host and parasite distance matrices from CSV
buttMatrixAPaco <- as.matrix(read.csv("ButtSuperAMatrix.csv", row.names=1))
wolbMatrixAPaco <- as.matrix(read.csv("WolbSuperAMatrix.csv", row.names=1))

# matrices in the same order
# Convert to distance objects
buttDistA <- as.dist(buttMatrixAPaco)
wolbDistA <- as.dist(wolbMatrixAPaco)

# run PACo
paco_resultA <- procrustes(buttDistA, wolbDistA, symmetric = TRUE)

# visualize
plot(paco_resultA, main="Procrustes Analysis (PACo)")

# test for significance
protest_resultA <- protest(buttDistA, wolbDistA, permutations = 999)

# Print results
print(protest_resultA)



##################### GLM to test correlation between percent covered and extracted loci #####################
# This is just one example of how the test was done, especifically testing for correlation of the distribution of extracted loci for samples above and below a threshold of 30% evenness
# of coverage i.e., distribution of reads mapped along Wolbachia reference genome. This test was also done for complete BUSCOs.

# Load necessary packages

install.packages(c("ggplot2", "MASS", "mgcv", "sandwich", "lmtest", "pscl", "AER"))

library(ggplot2)
library(MASS)  # For glm.nb (Negative Binomial GLM)
library(mgcv)  # For Generalized Additive Models
library(sandwich)  # For robust standard errors
library(lmtest)  # For hypothesis testing
library(pscl)  # For zero-inflated models

# Load data
setwd("C:\\Users\\pedro\\OneDrive\\Documents\\PhD\\Projects\\Eudaminae_MacroWolbachia\\PercentCoveredStats")
# Read your data
dataB <- read.csv("percentCoveredB.csv")  # Replace with your actual file

# Create the threshold indicator variable
dataB$AboveThresholdB <- ifelse(dataB$PercentCoveredB >= 30, 1, 0)

# Fit the Quasi-Poisson Model
quasi_poisson_modelExtractedB <- glm(ExtractedContigsB ~ AboveThresholdB, family = quasipoisson, data = dataB)

# Model summaries
cat("Quasi-Poisson Model Summary:\n")
summary(quasi_poisson_modelExtractedB)

# Check for Overdispersion (Pearson Chi-Square Test)
overdispersion_test_quasiB <- sum(residuals(quasi_poisson_modelB, type = "pearson")^2) / quasi_poisson_modelB$df.residual
cat("\nQuasi-Poisson Overdispersion Test (Pearson Chi-Square):", overdispersion_test_quasiB, "\n")

# If the result is much greater than 1, it indicates overdispersion.
# For a better diagnostic, we can use the "dispersiontest" from the AER package:
# install.packages("AER")
# library(AER)
# dispersiontest(quasi_poisson_modelA)

# Residual Diagnostics for each model
par(mfrow = c(2, 2))

# Quasi-Poisson Model
plot(quasi_poisson_modelA)