如何编写支持FASTA格式参数的R氨基酸分子量计算函数

The original code works great for raw amino acid strings, but it doesn’t account for FASTA’s structure (headers, multiple sequences, line breaks). Let’s fix that by building a function that parses FASTA input first, then calculates molecular weights for each sequence.

Modified Code #

# Monoisotopic amino acid weights (excluding water for peptide bonds)
aa_weights <- c(
  I = 131.1736, L = 131.1736, K = 146.1882, M = 149.2124, F = 165.19,
  T = 119.1197, W = 204.2262, V = 117.1469, R = 174.2017, H = 155.1552,
  A = 89.0935, N = 132.1184, D = 133.1032, C = 121.159, E = 147.1299,
  Q = 146.1451, G = 75.0669, P = 115.131, S = 105.093, Y = 181.1894
)

# Calculate molecular weights from FASTA input (file path or raw string)
calculate_fasta_mw <- function(fasta_input) {
  # Read input: handle both file paths and raw FASTA strings
  if (file.exists(fasta_input)) {
    fasta_lines <- readLines(fasta_input)
  } else {
    fasta_lines <- strsplit(fasta_input, "\n")[[1]]
  }

  # Parse FASTA: separate headers and sequences
  headers <- c()
  sequences <- c()
  current_seq <- ""

  for (line in fasta_lines) {
    line <- trimws(line)
    if (startsWith(line, ">")) {
      # Save the previous sequence if we're switching to a new header
      if (current_seq != "") {
        sequences <- c(sequences, current_seq)
        current_seq <- ""
      }
      # Store header (remove the leading ">")
      headers <- c(headers, substring(line, 2))
    } else if (line != "") {
      # Append sequence lines (convert to uppercase to match weight keys)
      current_seq <- paste0(current_seq, toupper(line))
    }
  }

  # Add the final sequence after loop ends
  if (current_seq != "") {
    sequences <- c(sequences, current_seq)
  }

  # Calculate MW for each sequence
  mw_results <- sapply(sequences, function(seq) {
    if (nchar(seq) == 0) return(NA)
    # Split sequence into individual amino acids
    aa_list <- unlist(strsplit(seq, ""))
    # Sum weights, subtract water lost for each peptide bond (18.01528 is precise H2O mass)
    sum(aa_weights[aa_list]) - ((nchar(seq) - 1) * 18.01528)
  })

  # Name results with FASTA headers (or default names if headers are missing)
  if (length(headers) == length(sequences)) {
    names(mw_results) <- headers
  } else {
    names(mw_results) <- paste0("Sequence_", 1:length(sequences))
  }

  return(mw_results)
}

Key Improvements #

• FASTA Parsing: Handles headers, multiple sequences, and line breaks in FASTA input.
• Flexible Input: Works with both raw FASTA strings and file paths.
• Robustness: Converts sequences to uppercase (so lowercase inputs work) and skips empty lines.
• Clear Output: Returns named results matching FASTA headers for easy reference.

Usage Examples #

# Example 1: Use a raw FASTA string
sample_fasta <- ">My_Protein_1
MAKELVING
>My_Protein_2
HELLOWORLD"

calculate_fasta_mw(sample_fasta)
# Output:
# My_Protein_1 My_Protein_2 
#      794.913     1126.296 

# Example 2: Read from a FASTA file
calculate_fasta_mw("path/to/your/sequences.fasta")

Optional Add-On #

If you want to filter out invalid amino acids (characters not in the weight list), add this line after splitting the sequence:

aa_list <- aa_list[aa_list %in% names(aa_weights)]

内容的提问来源于stack exchange，提问作者Mayank Singh Rajput