本文主要是介绍与PC1显著相关的基因 | p值计算，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
1. 相关系数的显著性
t=r*sqrt(n-2) / sqrt(1-r**2)
其中，统计量t符合自由度为 n-2 的t分布。
2. 与PC1显著相关的基因
就是求相关系数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为显著的基因即可。
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) 样本大小 nn <- 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] ** 2correlations <- 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 yestail(withPC1)
table(withPC1$sig)
# no yes
#1239 761
table(withPC1$sig, withPC1$Loading>0)
# FALSE TRUE
#no 853 386
#yes 650 111hist(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()
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.
- txtBlog::R/Rs1 统计
这篇关于与PC1显著相关的基因 | p值计算的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
