lvzhaopeng的数字花园

R语言源码那些事（一） KSEA激酶分析

5分钟阅读

缘由

不懂原理的生信人，不是一个好的调包侠。
这个系列不是教程系列 而是原理系列，教程地址见Ref

碎碎念

我的梦想，值得我本人去争取，我今天的生活，绝不是我昨天生活的冷淡抄袭。——《红与黑》

Question 答案见底部

  1. KSEA激酶分析回答了什么问题?
  2. KSEA激酶的核心原理是什么?

Ref

  1. The KSEA App 原始出处
  2. 使用KSEA App进行激酶底物富集分析（KSEA） 使用教程 友链

输入文件

KSData

PX

代码

source <- function(KSData, PX, NetworKIN, NetworKIN.cutoff, m.cutoff,
                   p.cutoff)
{ 
# 1. 整理数据
  if (length(grep(";", PX$Residue.Both)) == 0) {
    new <- PX
    colnames(new)[c(2, 4)] <- c("SUB_GENE", "SUB_MOD_RSD")
    new$log2FC <- log2(abs(as.numeric(as.character(new$FC))))
    new <- new[complete.cases(new$log2FC), ] # 删除log2FC为NA的行
  } else {
    double <- PX[grep(";", PX$Residue.Both), ]
    residues <- as.character(double$Residue.Both)
    residues <- as.matrix(residues, ncol = 1)
    split <- strsplit(residues, split = ";")
    x <- sapply(split, length)
    single <- data.frame(
      Protein = rep(double$Protein, x),
      Gene = rep(double$Gene, x),
      Peptide = rep(double$Peptide, x),
      Residue.Both = unlist(split),
      p = rep(double$p, x),
      FC = rep(double$FC, x)
    )
    new <- PX[-grep(";", PX$Residue.Both), ]
    new <- rbind(new, single)
    colnames(new)[c(2, 4)] <- c("SUB_GENE", "SUB_MOD_RSD")
    new$log2FC <- log2(abs(as.numeric(as.character(new$FC))))
    new <- new[complete.cases(new$log2FC), ]
  }
#2. 选择数据率
  if (NetworKIN == TRUE) {
    KSData.filtered <- KSData[grep("[a-z]", KSData$Source), ]
    KSData.filtered <- KSData.filtered[(KSData.filtered$networkin_score >=
      NetworKIN.cutoff), ]
  } else {
    KSData.filtered <- KSData[grep("PhosphoSitePlus", KSData$Source), ]
  }
#3. 合并数据库和鉴定的位点
  # 合并的逻辑，认为底物名称和位点位置完全一致
  KSData.dataset <- merge(KSData.filtered, new)
  KSData.dataset <- KSData.dataset[order(KSData.dataset$GENE), ]
  KSData.dataset$Uniprot.noIsoform <- sapply(
    KSData.dataset$KIN_ACC_ID,
    function(x) unlist(strsplit(as.character(x), split = "-"))[1]
  )
  KSData.dataset.abbrev <- KSData.dataset[, c(
    5, 1, 2, 16:19,
    14
  )]
 
  # 最终使用的数据框包含 激酶基因 底物基因 底物 肽段 FC log2FC source列
  colnames(KSData.dataset.abbrev) <- c(
    "Kinase.Gene", "Substrate.Gene",
    "Substrate.Mod", "Peptide", "p", "FC", "log2FC", "Source"
  )
  KSData.dataset.abbrev <- KSData.dataset.abbrev[order(
    KSData.dataset.abbrev$Kinase.Gene,
    KSData.dataset.abbrev$Substrate.Gene, KSData.dataset.abbrev$Substrate.Mod,
    KSData.dataset.abbrev$p
  ), ]
#计算相关指标
  #按照分组进行计算同一个激酶 底物 修饰位点 来源计算取均值
  KSData.dataset.abbrev <- aggregate(
    log2FC ~ Kinase.Gene +
      Substrate.Gene + Substrate.Mod + Source,
    data = KSData.dataset.abbrev,
    FUN = mean
  )
  KSData.dataset.abbrev <- KSData.dataset.abbrev[order(KSData.dataset.abbrev$Kinase.Gene), ]
  kinase.list <- as.vector(KSData.dataset.abbrev$Kinase.Gene)
  #计算每个激酶的底物数
  kinase.list <- as.matrix(table(kinase.list)) 
  #计算每个激酶平均log2FC
  Mean.FC <- aggregate(log2FC ~ Kinase.Gene,
    data = KSData.dataset.abbrev,
    FUN = mean
  ) 
  Mean.FC <- Mean.FC[order(Mean.FC[, 1]), ]
  Mean.FC$mS <- Mean.FC[, 2]
  Mean.FC$Enrichment <- Mean.FC$mS / abs(mean(new$log2FC, na.rm = T))
  Mean.FC$m <- kinase.list
#进行富集或者匮乏分析
  ###核心原理代码###
  Mean.FC$z.score <- ((Mean.FC$mS - mean(new$log2FC, na.rm = T)) *
    sqrt(Mean.FC$m)) / sd(new$log2FC, na.rm = T)
  Mean.FC$p.value <- pnorm(-abs(Mean.FC$z.score)) # 计算单尾概率
  Mean.FC$FDR <- p.adjust(Mean.FC$p.value, method = "fdr")
  Mean.FC.filtered <- Mean.FC[(Mean.FC$m >= m.cutoff), -2]
  Mean.FC.filtered <- Mean.FC.filtered[order(Mean.FC.filtered$z.score), ]
 
#作图
  plot.height <- length(Mean.FC.filtered$z.score)^0.55
  Mean.FC.filtered$color <- "black"
  Mean.FC.filtered[(Mean.FC.filtered$p.value < p.cutoff) &
    (Mean.FC.filtered$z.score < 0), ncol(Mean.FC.filtered)] <- "blue"
  Mean.FC.filtered[(Mean.FC.filtered$p.value < p.cutoff) &
    (Mean.FC.filtered$z.score > 0), ncol(Mean.FC.filtered)] <- "red"
  tiff("KSEA Bar Plot.tiff", width = 6 * 300, height = 300 *
    plot.height, res = 300, pointsize = 13)
  par(mai = c(1, 1, 0.4, 0.4))
  barplot(as.numeric(Mean.FC.filtered$z.score),
    col = Mean.FC.filtered$color,
    border = NA, xpd = F, cex.names = 0.6, cex.axis = 0.8,
    xlab = "Kinase z-score", names.arg = Mean.FC.filtered$Kinase.Gene,
    horiz = T, las = 1
  )
  dev.off()
  write.csv(KSData.dataset.abbrev,
    file = "Kinase-Substrate Links.csv",
    quote = F, row.names = F
  )
  write.csv(Mean.FC[order(Mean.FC$Kinase.Gene), -ncol(Mean.FC)],
    file = "KSEA Kinase Scores.csv", quote = F, row.names = F
  )
}

Answer

  1. KSEA激酶分析回答了什么问题?
  • 哪些激酶的功能活性可能被激活或抑制了（基于底物磷酸化）
  1. KSEA激酶的核心原理是什么?
  • KSEA（Kinase-Substrate Enrichment Analysis）的核心原理是：通过比较特定激酶下游底物的磷酸化水平（log2FC）与全局平均水平，统计判断这些底物是否整体显著上调或下调，从而间接推断该激酶是否活性增强或抑制。它不依赖激酶自身是否被检测到，而是基于激酶-底物数据库和质谱磷酸化数据，利用 z-score 和 p 值评估每个激酶的调控活性变化趋势。

关系图谱

反向链接