与PC1显著相关的基因 | p值计算

1. 相关系数r的显著性

t统计量 t=r*sqrt(n-2) / sqrt(1-r**2)，
符合自由度为 n-2 的t分布。

2. 与PC1显著相关的基因

(1) 直接计算

就是求相关系数r=cor(PC1_score, Xk)，其中

• PC1_score 长度为样品总数，是PC1 的loading * 每个变量的scale后的值
• Xk是第k个变量在每个样品的值

然后由r计算t统计量，及对用的p值，见上文1。

对每个变量都这么计算p值，做p值矫正，取 p_val_adj <0.05 或 0.01为显著的基因即可。

(2) 推导公式简化计算

参考文献[1]，loadings和r成比例：r = sqrt(lambda) * W / sd(Xk)

• lambda 是PC1对应的特征值，也是特征向量所在方向的方差。
• W是权重，也就是loading
• Xk是原始autoscale后数据的第k个变量的值向量。
• 和硬算cor.test(PC1_score, scale.data.gene)结果是一致的，且节省很多算力。
• 假设 autoscale 后的数据的 sd=1，可以忽略分母。我测试发现假设不成立，有少量sd小于1，不知道Seurat怎么处理的… 另外：忽略分母会造成 显著基因数的略微减少。

3. 由Seurat对象求与某PC显著的基因

#' Calc p value of each gene with one given PC
#'
#' @param object Seurat object
#' @param pc_num like 1 for "PC_1"
#' @param p.adj.method default fdr
#' @param verbose default no output
#' @param p.val.threshold p value significant threshold
#'
#' @return
#' @export
#'
#' @examples
#' withPC1=SigGenesWithPC(scObj, pc_num = 2)
SigGenesWithPC=function(object, pc_num=1, p.adj.method="fdr", p.val.threshold=0.05, verbose=F){
  #rk=sqrt(lambda)*W / sd(Xk)
  # rk 是PC1和第k个变量的r值
  # lambda 是PC1对应的特征值
  # W是PC1对应的每个变量的 loading
  # Xk是第k个变量的autoscale后的值
  
  # 1)提取 PC1 的因子加载
  loadings <- object@reductions$pca@feature.loadings[,pc_num]  # PC1 的加载
  
  # 2) 样本大小 n
  n <- ncol(scale.data); n  # 样本大小
  
  # 3) autoscale 后的数据
  scale.data=object@assays$RNA@scale.data[object@assays$RNA@var.features,]
  
  # 4)计算相关系数 r 和
  lambda = object@reductions$pca@stdev[pc_num] ** 2
  correlations <- sqrt(lambda) * as.numeric(loadings)
  
  sd.arr=apply(scale.data, 1, sd)
  if(verbose) { hist(sd.arr, n=100) } #并不是都是1，why?
  
  for(i in 1:length(correlations)){
    correlations[i] = correlations[i] / sd.arr[i]
  }
  
  # 5) 计算 t 统计量
  t_statistics <- (correlations * sqrt(n - 2)) / sqrt(1 - correlations^2)
  
  # 6) 计算 p 值
  p_values <- 2 * pt(-abs(t_statistics), df = n - 2)
  if(verbose) { hist(p_values, n=100) }
  
  # 7)创建结果数据框
  dat.use <- data.frame(Variable = names(loadings), 
                        Loading = loadings, 
                        t_statistic = t_statistics, 
                        pc=pc_num,
                        p_val = p_values)
  # 8)p 值矫正
  dat.use$p_val_adj = p.adjust(dat.use$p_val, method = p.adj.method)
  
  # 筛选显著变量（例如 p < 0.05）
  dat.use$sig=ifelse(dat.use$p_val_adj < p.val.threshold, "yes", "no")
  
  # 9) order by Loading 
  dat.use=dat.use[order(-dat.use$Loading),]
  dat.use
}

# scObj 为pbmc3k的Seurat对象，v4。

withPC1=SigGenesWithPC(scObj, pc_num = 2, p.val.threshold = 0.01)
head(withPC1)
#         Variable   Loading t_statistic pc         p_val     p_val_adj sig
#CD79A       CD79A 0.1244749    34.49119  2 2.703414e-216 6.007586e-214 yes
#MS4A1       MS4A1 0.1130108    30.16768  2 1.561846e-172 2.402839e-170 yes
#TCL1A       TCL1A 0.1061570    27.79212  2 1.035920e-149 1.218729e-147 yes
#HLA-DQA1 HLA-DQA1 0.1031493    26.79132  2 2.104518e-140 2.338353e-138 yes
#HLA-DRA   HLA-DRA 0.1022300    26.49010  2 1.219714e-137 1.283910e-135 yes
#HLA-DQB1 HLA-DQB1 0.1007367    26.00524  2 3.128123e-133 2.979165e-131 yes

tail(withPC1)
table(withPC1$sig)
#  no  yes 
#1239  761
table(withPC1$sig, withPC1$Loading>0)
#    FALSE TRUE
#no    853  386
#yes   650  111

hist(withPC1$p_val_adj, n=100)

library(ggplot2)
ggplot(withPC1, aes(x=Variable, y=Loading, fill=sig))+
  geom_col()+
  scale_x_discrete(limits=withPC1$Variable, labels = NULL)+
  scale_fill_manual(values = c("red", "grey"), breaks = c("yes","no") )+
  ggtitle("p.adj<0.01")

# get sig genes: all
withPC1[which(withPC1$sig=="yes"), ]$Variable |> jsonlite::toJSON()
# get sig genes: pos sig
withPC1[which(withPC1$sig=="yes" & withPC1$Loading>0), ]$Variable |> jsonlite::toJSON()
# get sig genes: neg sig
withPC1[which(withPC1$sig=="yes" & withPC1$Loading<0), ]$Variable |> jsonlite::toJSON()

图中可见，负相关的基因个数比正相关的多。

富集分析结果

和PC1显著相关的761个基因：
其中，正相关：111，负相关：650个。

---

正相关：血液发育、细胞周期、体液免疫响应

---

负相关：细胞周期、有丝分裂G1/S转换、蛋白成熟、有丝分裂细胞周期

Ref

• [1] Yamamoto H., Fujimori T., Sato H., Ishikawa G., Kami K., Ohashi Y. (2014). “Statistical hypothesis testing of factor loading in principal component analysis and its application to metabolite set enrichment analysis”. BMC Bioinformatics, (2014) 15(1):51. 这文献写的公式有问题，1-6公式中根号括的过宽。文章声明可以用loading定义r，我认为不能省略系数 sqrt(lambda)
• txtBlog::R/Rs1 统计