# Install and load the Biostrings package if you haven't already
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("Biostrings")

library(Biostrings)

# --- Step 1: Define Sample Biological Sequences ---
# For DNA sequences
dna_sequences <- DNAStringSet(c(
  "ATGCGGTAGCTGA",
  "CCGTAGCATCGG",
  "GATCGTACGTAGCATG"
))

# For Protein sequences
protein_sequences <- AAStringSet(c(
  "MPEPAVLGSQAYP",
  "KRASALGGVPES",
  "YFGSPWEEVAHK"
))

cat("\n--- Sample DNA Sequences ---\n")
print(dna_sequences)
cat("\n--- Sample Protein Sequences ---\n")
print(protein_sequences)

# --- Step 2: Basic Feature Extraction (Length) ---
# Length is a fundamental numerical feature for any sequence.

dna_lengths <- width(dna_sequences)
protein_lengths <- width(protein_sequences)

cat("\n--- DNA Sequence Lengths ---\n")
print(dna_lengths)
cat("\n--- Protein Sequence Lengths ---\n")
print(protein_lengths)

# --- Step 3: Extracting Compositional Features (GC Content for DNA, Amino Acid Counts for Protein) ---

# Calculate GC content for DNA sequences
gc_content_dna <- letterFrequency(dna_sequences, letters = c("G", "C"), as.prob = TRUE) %>%
  rowSums()

cat("\n--- DNA GC Content ---\n")
print(gc_content_dna)

# Calculate amino acid frequencies for protein sequences
# This gives a matrix where each row is a sequence and columns are amino acids
aa_frequencies_protein <- letterFrequency(protein_sequences, letters = AA_ALPHABET, as.prob = TRUE)

cat("\n--- Protein Amino Acid Frequencies (first 5 columns) ---\n")
print(head(aa_frequencies_protein[, 1:5])) # Display only first 5 columns for brevity

# --- Step 4: Combine into a Data Frame for R Analysis ---
# For DNA sequences
dna_features_df <- data.frame(
  Sequence = as.character(dna_sequences),
  Length = dna_lengths,
  GC_Content = gc_content_dna
)

cat("\n--- Combined DNA Features Data Frame ---\n")
print(dna_features_df)

# For Protein sequences
protein_features_df <- data.frame(
  Sequence = as.character(protein_sequences),
  Length = protein_lengths
)
# Add amino acid frequencies to the protein features data frame
protein_features_df <- cbind(protein_features_df, aa_frequencies_protein)

cat("\n--- Combined Protein Features Data Frame (first 7 columns) ---\n")
print(head(protein_features_df[, 1:7])) # Display only first 7 columns for brevity


# We now have numerical features ready for statistical analysis in R.
# This forms the foundation for more complex feature engineering.

# Ensure Biostrings is loaded
library(Biostrings)

# Re-define sequences for clean execution
dna_sequences <- DNAStringSet(c(
  "ATGCGGTAGCTGA",
  "CCGTAGCATCGG",
  "GATCGTACGTAGCATG"
))
protein_sequences <- AAStringSet(c(
  "MPEPAVLGSQAYP",
  "KRASALGGVPES",
  "YFGSPWEEVAHK"
))

# --- Step 1: K-mer Frequencies (Compositional) ---
# K-mers are contiguous sequences of length k.
# For DNA: Dinucleotide frequencies (k=2)
# oligonucleotideFrequency computes frequencies of k-mers
dna_dimer_freq <- oligonucleotideFrequency(dna_sequences, width = 2, as.prob = TRUE)
cat("\n--- DNA Dinucleotide Frequencies ---\n")
print(dna_dimer_freq)

# For Proteins: Di-peptide frequencies (k=2)
# aminoAcidFrequency can also be used for k-mers but oligonucleotideFrequency is more generic
# We'll demonstrate specific dipeptide counts for clarity, though it can be very high-dimensional
# Function to calculate dipeptide frequencies for a single protein sequence
get_dipeptide_freq <- function(seq_str) {
  dipeptides <- Biostrings::oligonucleotideFrequency(AAStringSet(seq_str), width = 2, as.prob = TRUE)
  return(as.vector(dipeptides))
}

# Apply to all protein sequences and combine into a matrix/dataframe
protein_dipeptide_list <- lapply(as.character(protein_sequences), get_dipeptide_freq)
# Pad with NAs if lengths are different due to sequence variations (rare for fixed k-mer length)
# For dipeptides, the number of features is fixed (20*20 = 400), so padding is not needed if all sequences are long enough.
protein_dipeptide_mat <- do.call(rbind, protein_dipeptide_list)
colnames(protein_dipeptide_mat) <- Biostrings::oligonucleotideFrequency(AAStringSet("AA"), width=2) %>% colnames()

cat("\n--- Protein Di-peptide Frequencies (first 5 columns) ---\n")
print(head(protein_dipeptide_mat[, 1:5]))

# --- Step 2: Physicochemical Features for Proteins ---
# Molecular Weight (MW) and Isoelectric Point (pI) are key physicochemical properties.
# These are often derived from the sequence composition and properties of amino acids.

# Calculate Molecular Weight
# Biostrings provides a direct way to calculate MW for AAStringSet objects
protein_mw <- letterFrequency(protein_sequences, AA_ALPHABET, as.prob = FALSE) %*%
              AA_MASS # AA_MASS is a predefined vector of amino acid masses

cat("\n--- Protein Molecular Weights ---\n")
print(protein_mw)

# Calculate Isoelectric Point (pI)
# This requires a more complex algorithm. We'll simulate using a simplified proxy or external package for actual calculation.
# For a true pI calculation, you would typically use dedicated bioinformatics packages/tools
# For demonstration, we'll use a placeholder or a simplified charge calculation based on acidic/basic residues.
# A common approach uses 'seqinr' or 'Peptides' package for actual pI, but we will demonstrate a simplified proxy.

# Simplified charge estimation (proxy for pI consideration) - count basic (K,R) vs acidic (D,E) residues
get_charge_proxy <- function(seq_str) {
  basic_count <- sum(letterFrequency(AAStringSet(seq_str), letters=c("K", "R")))
  acidic_count <- sum(letterFrequency(AAStringSet(seq_str), letters=c("D", "E")))
  return(basic_count - acidic_count) # Net charge proxy
}

protein_charge_proxy <- sapply(as.character(protein_sequences), get_charge_proxy)
cat("\n--- Protein Charge Proxy (Simplified) ---\n")
print(protein_charge_proxy)

# --- Step 3: Combine all features into a unified data frame ---
protein_features_engineered_df <- data.frame(
  Sequence = as.character(protein_sequences),
  Length = width(protein_sequences),
  Molecular_Weight = protein_mw,
  Charge_Proxy = protein_charge_proxy
)

# Add di-peptide frequencies. Ensure column names are unique.
protein_features_engineered_df <- cbind(protein_features_engineered_df, protein_dipeptide_mat)

cat("\n--- Engineered Protein Features Data Frame (first 10 columns) ---\n")
print(head(protein_features_engineered_df[, 1:10]))

# This refined data frame is now ready for more sophisticated statistical analysis.

# Ensure Biostrings is loaded
library(Biostrings)

# Re-define sequences for clean execution
protein_sequences <- AAStringSet(c(
  "MPEPAVLGSQAYP",
  "KRASALGGVPES",
  "YFGSPWEEVAHK"
))

# --- Step 1: One-Hot Encoding for Short Sequence Segments ---
# One-hot encoding is suitable for fixed-length segments where position matters.
# Let's take a fixed-length segment (e.g., first 5 amino acids) for demonstration.
# This generates a sparse matrix where each position and amino acid is a binary feature.

encode_one_hot <- function(seq_vec, alphabet) {
  max_len <- max(nchar(seq_vec))
  encoded_list <- list()
  for (i in seq_along(seq_vec)) {
    s <- strsplit(seq_vec[i], split = "")[[1]]
    row_encoding <- numeric(length(alphabet) * max_len)
    names(row_encoding) <- paste0("Pos", rep(1:max_len, each = length(alphabet)), "_", rep(alphabet, max_len))
    for (j in seq_along(s)) {
      aa_idx <- match(s[j], alphabet)
      if (!is.na(aa_idx)) {
        col_idx <- (j - 1) * length(alphabet) + aa_idx
        row_encoding[col_idx] <- 1
      }
    }
    encoded_list[[i]] <- row_encoding
  }
  do.call(rbind, encoded_list)
}

# Use first 5 AAs for one-hot encoding example
short_proteins <- substr(as.character(protein_sequences), 1, 5)
aa_alphabet <- AA_ALPHABET # Use the predefined Biostrings alphabet

one_hot_features <- encode_one_hot(short_proteins, aa_alphabet)

cat("\n--- One-Hot Encoded Features (first 5 amino acids) ---\n")
print(head(one_hot_features[, 1:10])) # Display first 10 columns for brevity

# --- Step 2: Physicochemical Property Matrices ---
# Transform sequences into numerical matrices based on amino acid properties.
# Use Kyte-Doolittle hydrophobicity scale as an example.

# Define a simple hydrophobicity scale (Kyte-Doolittle example values)
hydrophobicity_scale <- c(
  'A'=1.8, 'R'=-4.5, 'N'=-3.5, 'D'=-3.5, 'C'=2.5, 'Q'=-3.5, 'E'=-3.5, 'G'=-0.4, 'H'=-3.2, 
  'I'=4.5, 'L'=3.8, 'K'=-3.9, 'M'=1.9, 'F'=2.8, 'P'=-1.6, 'S'=-0.8, 'T'=-0.7, 'W'=-0.9, 
  'Y'=-1.3, 'V'=4.2
)

# Function to convert a protein sequence into a hydrophobicity vector
get_hydrophobicity_vector <- function(seq_str, scale) {
  s <- strsplit(seq_str, split = "")[[1]]
  return(unname(scale[s]))
}

# Apply to all protein sequences
hydrophobicity_vectors_list <- lapply(as.character(protein_sequences), get_hydrophobicity_vector)

# Pad vectors to the maximum length for matrix creation
max_protein_len <- max(width(protein_sequences))

padded_hydrophobicity_mat <- do.call(rbind, lapply(hydrophobicity_vectors_list, function(vec) {
  c(vec, rep(0, max_protein_len - length(vec))) # Pad with 0s for missing positions
}))

colnames(padded_hydrophobicity_mat) <- paste0("Hydro_Pos", 1:max_protein_len)

cat("\n--- Padded Hydrophobicity Vectors ---\n")
print(padded_hydrophobicity_mat)

# --- Step 3: Derived Structural Proxies (Simplified Helix Propensity) ---
# This is a highly simplified example. True secondary structure prediction 
# involves complex algorithms (e.g., Chou-Fasman, GOR, neural networks).
# We'll create a very basic proxy based on known helix-forming amino acids.

helix_propensity_scale <- c(
  'A'=1.42, 'L'=1.21, 'M'=1.45, 'E'=1.53, 'K'=1.16, 'H'=1.24, 'Q'=1.11, 'R'=1.21, 
  'C'=0.66, 'D'=0.76, 'G'=0.57, 'N'=0.73, 'P'=0.59, 'S'=0.79, 'T'=0.82, 'V'=1.06, 
  'W'=1.08, 'Y'=0.61, 'F'=1.13, 'I'=1.00 # Placeholder for non-specific values
)

get_helix_propensity_sum <- function(seq_str, scale) {
  s <- strsplit(seq_str, split = "")[[1]]
  propensities <- unname(scale[s])
  return(sum(propensities, na.rm = TRUE))
}

protein_helix_propensity_sum <- sapply(as.character(protein_sequences), get_helix_propensity_sum, scale=helix_propensity_scale)

cat("\n--- Summed Helix Propensity Proxy ---\n")
print(protein_helix_propensity_sum)

# --- Step 4: Combine into a final data frame ---
advanced_protein_features_df <- data.frame(
  Sequence = as.character(protein_sequences),
  OneHot_Sample = paste(one_hot_features[,1:5] %*% rep(1,5)), # Simplistic representation for df
  Sum_Hydrophobicity = rowSums(padded_hydrophobicity_mat),
  Helix_Propensity_Sum = protein_helix_propensity_sum
)

# For illustrative purposes, we're not adding all one-hot and hydrophobicity vectors directly due to dimensionality.
# In a real scenario, you'd add the full matrices or use dimensionality reduction techniques.
# Here we just show a few summary features or subset for clarity.

cat("\n--- Advanced Protein Features Data Frame ---\n")
print(advanced_protein_features_df)

# Ensure Biostrings is loaded and previous data frames exist (e.g., protein_features_engineered_df)
library(Biostrings)

# Re-create a sample protein features data frame with some numerical features
# For this part, we'll create a synthetic dataframe to ensure we have numerical data to work with
set.seed(123)
sample_proteins <- AAStringSet(c(
  "MPEPAVLGSQAYP", "KRASALGGVPES", "YFGSPWEEVAHK", "LGFSYQARAS", "VPLEGGATY"
))

sample_df <- data.frame(
  Sequence = as.character(sample_proteins),
  Length = width(sample_proteins),
  Molecular_Weight = c(1400, 1350, 1500, 1100, 1250) + rnorm(5, sd=50),
  Charge_Proxy = c(-2, 3, 0, 1, -1) + rnorm(5, sd=0.5),
  Hydrophobicity_Sum = c(10, -5, 15, 8, 20) + rnorm(5, sd=2),
  DiPeptide_AA_Ratio = c(0.1, 0.05, 0.12, 0.08, 0.09) + rnorm(5, sd=0.01)
)

cat("\n--- Original Sample Feature Data Frame ---\n")
print(sample_df)

# --- Step 1: Normalization and Scaling ---
# Crucial for algorithms sensitive to feature scales (e.g., PCA, SVM, neural networks).
# We will apply Z-score normalization (mean=0, sd=1) to numerical columns.

numerical_cols <- c("Length", "Molecular_Weight", "Charge_Proxy", "Hydrophobicity_Sum", "DiPeptide_AA_Ratio")

# Apply scale() function for Z-score normalization
scaled_features <- scale(sample_df[, numerical_cols])

# Replace original numerical columns with scaled ones in a new data frame
normalized_df <- sample_df
normalized_df[, numerical_cols] <- scaled_features

cat("\n--- Z-score Normalized Feature Data Frame ---\n")
print(normalized_df)

# --- Step 2: Simple Feature Selection (e.g., based on variance or correlation) ---
# Identify and potentially remove features with very low variance or high correlation.
# For demonstration, let's assume we want to remove 'DiPeptide_AA_Ratio' if its variance is below a threshold.

# Calculate variance of scaled features
feature_variances <- apply(scaled_features, 2, var)
cat("\n--- Variances of Scaled Features ---\n")
print(feature_variances)

# Threshold example: remove features with variance less than 0.1 (arbitrary for demonstration)
low_variance_features <- names(feature_variances[feature_variances < 0.1])
cat("\n--- Features with low variance (removed) ---\n")
print(low_variance_features)

# Create a data frame without low variance features
filtered_df_variance <- normalized_df[, !(colnames(normalized_df) %in% low_variance_features)]
cat("\n--- Feature Data Frame after Variance Filtering ---\n")
print(filtered_df_variance)

# --- Step 3: Data Splitting for Model Validation ---
# Essential for unbiased model evaluation.
# Split data into training and testing sets (e.g., 70% train, 30% test).

library(caret) # Useful for data splitting

# Assuming a 'target_variable' exists (e.g., protein activity score)
# For this example, let's add a dummy target variable.
filtered_df_variance$Target_Activity <- c(10, 25, 15, 30, 20)

# Create an index for training set
train_index <- createDataPartition(filtered_df_variance$Target_Activity, p = 0.7, list = FALSE)

# Create training and testing sets
training_set <- filtered_df_variance[train_index, ]
testing_set <- filtered_df_variance[-train_index, ]

cat("\n--- Training Set (first 2 rows) ---\n")
print(head(training_set, 2))
cat("\n--- Testing Set (first 2 rows) ---\n")
print(head(testing_set, 2))

# This process ensures features are prepared for robust statistical modeling.