当前位置：首页 > news >正文

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

news 来源：原创 2024/9/21 22:31:03

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 统计