【生信分析】免疫组库基础分析8-CDR3氨基酸使用频率

# 加载必要的可视化包
library(ggplot2)
# 定义标准20种氨基酸（生物化学标准）
standard_aas <- c("G", "A", "S", "C", "D", "P", "N", "T", "E", "V",
                  "Q", "H", "M", "I", "L", "K", "R", "F", "Y", "W")
# 创建包含10个虚拟CDR3序列及其频率的数据框
# 序列长度在8-20个氨基酸之间，频率总和为100%
seq_data <- data.frame(
  Sequence = c(
    "CASSLGNTGQLY-F",   # 序列1：长度13
    "CASSQETGQ*YF",     # 序列2：长度11
    "CASSEHGQLYF",     # 序列3：长度11
    "CASSLkQGNTIYF",   # 序列4：长度13
    "CASSQPGGQETQYF",  # 序列5：长度14
    "CASSLDVNTGQLYF",  # 序列6：长度13
    "CASSLEGANTGQFY",  # 序列7：长度13
    "CASSLGQWGNTIYF",  # 序列8：长度14
    "CASSQGGNKEQYF",   # 序列9：长度12
    "CASSLGMGQGRNTIYF" # 序列10：长度15
  ),
  Frequency = c(15, 12, 10, 8, 20, 5, 10, 8, 7, 5) # 各序列占比(%)
)
# 验证频率总和是否为100%
if (sum(seq_data$Frequency) != 100) {
  stop("错误：序列频率总和必须等于100%")
}
# [3] 核心计算函数（分步实现）
calculate_aa_stats <- function(seq_df, aas) {
  # 初始化结果数据框（包含所有统计维度）
  result <- data.frame(
    AminoAcid = aas,
    RawCount = 0,        # 原始出现次数
    RawFrequency = 0,    # 原始频率（不考虑克隆权重）
    WeightedCount = 0,   # 加权计数（考虑克隆频率）
    WeightedFrequency = 0,  # 最终加权频率
    stringsAsFactors = FALSE
  )
  # [3.1] 原始计数计算
  total_aa <- 0  # 总氨基酸计数
  for(seq in seq_df$Sequence) {
    chars <- strsplit(seq, "")[[1]]  # 拆分为单字符
    counts <- table(factor(chars, levels=aas))  # 强制统计所有氨基酸
    result$RawCount <- result$RawCount + as.numeric(counts)
    total_aa <- total_aa + sum(counts)  # 累计总氨基酸数
  }
  # 计算原始频率（百分比）
  result$RawFrequency <- round(result$RawCount / total_aa * 100, 2)
  # [3.2] 加权频率计算
  total_weighted <- 0  # 加权总和
  for(i in 1:nrow(seq_df)) {
    seq <- seq_df$Sequence[i]
    weight <- seq_df$Frequency[i] / 100  # 转换为权重系数
    chars <- strsplit(seq, "")[[1]]
    counts <- table(factor(chars, levels=aas))
    weighted_counts <- as.numeric(counts) * weight  # 当前序列的加权贡献
    result$WeightedCount <- result$WeightedCount + weighted_counts
    total_weighted <- total_weighted + sum(weighted_counts)
  }
  # 计算加权频率（百分比）
  result$WeightedFrequency <- round(result$WeightedCount / total_weighted * 100, 2)
  return(result)
}
aa_stats <- calculate_aa_stats(seq_data, standard_aas)  # 计算统计量
# [5] 结果输出（按加权频率降序）
final_result <- aa_stats[order(-aa_stats$WeightedFrequency), 
                         c("AminoAcid", "RawFrequency", "WeightedFrequency")]
print(final_result)
# [6] 专业可视化
ggplot(final_result, aes(x=reorder(AminoAcid, -WeightedFrequency))) +
  geom_col(aes(y=RawFrequency, fill="原始频率"), 
           width=0.35, position=position_nudge(x=-0.18)) +
  geom_col(aes(y=WeightedFrequency, fill="加权频率"), 
           width=0.35, position=position_nudge(x=0.18)) +
  geom_text(aes(y=RawFrequency, label=sprintf("%.1f%%", RawFrequency)),
            vjust=-0.5, size=3, nudge_x=-0.18) +
  geom_text(aes(y=WeightedFrequency, label=sprintf("%.1f%%", WeightedFrequency)),
            vjust=-0.5, size=3, nudge_x=0.18) +
  scale_fill_manual(values=c("#4E79A7", "#F28E2B")) +
  labs(title="CDR3区域氨基酸使用频率对比",
       subtitle="原始频率 vs 克隆加权频率",
       x="氨基酸", y="频率(%)", fill="统计方法") +
  theme_minimal(base_size=12) +
  theme(
    axis.text.x=element_text(angle=45, hjust=1, face="bold"),
    legend.position="top",
    panel.grid.major.x=element_blank()
  )