python - 在 python 中使用超几何测试
问题描述
我有两个基因列表,我计算它们之间的交集。
我需要计算 p 值以假设 - 这些列表的交集是偶然发生的。
我尝试使用 Fisher 精确测试(scipy 函数)来实现它。
请注意,我需要一个单面的 p 值。
我的代码:
def main(gene_path1, gene_path2, pop_size):
genes1 = pd.read_csv(gene_path1, sep='\n', header=None)
genes2 = pd.read_csv(gene_path2, sep='\n', header=None)
intersection = pd.merge(genes1, genes2, how='inner').drop_duplicates([0])
len_genes1 = genes1[0].count()
len_genes2 = genes2[0].count()
len_intersection = intersection[0].count()
oddsratio, pvalue = stats.fisher_exact([[len_genes1 - len_intersection, len_genes1], [len_genes2 - len_intersection, len_genes2]], alternative='less')
print(f'Genes1 len: {len_genes1}, Genes2 len: {len_genes2}, Intersection: {len_intersection}, pvalue: {pvalue}')
为了简单起见,我使用了一个数字列表(不是基因)。
由于它太长了,我不会复制整个文件,而是想象两个文件有很多随机数,由换行符分隔。
例如:
1
2
3
246
51451
...
问题是 - 我如何确定我正确指定了 Fisher 精确函数的参数?根据我要检查的假设是否正确?
我怀疑我做错了,但我不知道为什么。可能暗示出了什么问题 - 我知道人口规模应该是相关的,但我不确定在哪里使用它以及如何使用它。
任何线索或见解将不胜感激。
更新:
我试图以不同的方式实现它。
from scipy.stats import hypergeom as hg
import pandas as pd
def main(gene_path1, gene_path2, pop_size):
genes1 = pd.read_csv(gene_path1, sep='\n', header=None)
genes2 = pd.read_csv(gene_path2, sep='\n', header=None)
intersection = pd.merge(genes1, genes2, how='inner').drop_duplicates([0])
len_genes1 = genes1[0].count()
len_genes2 = genes2[0].count()
len_intersection = intersection[0].count()
pvalue = hg.cdf(int(len_intersection)-1, int(pop_size), int(len_genes1), int(len_genes2))
print(f'Genes1 len: {len_genes1}, Genes2 len: {len_genes2}, Intersection: {len_intersection}, p value: {pvalue})
我只是想知道我的论点是否在正确的地方,我该如何验证呢?
解决方案
这也应该有帮助: http: //pedagogix-tagc.univ-mrs.fr/courses/ASG1/practicals/go_statistics_td/go_statistics_td_2015.html
g = 75 ## Number of submitted genes
k = 59 ## Size of the selection, i.e. submitted genes with at least one annotation in GO biological processes
m = 611 ## Number of "marked" elements, i.e. genes associated to this biological process
N = 13588 ## Total number of genes with some annotation in GOTERM_BP_FAT.
n = N - m ## Number of "non-marked" elements, i.e. genes not associated to this biological process
x = 19 ## Number of "marked" elements in the selection, i.e. genes of the group of interest that are associated to this biological process
# Python
stats.hypergeom(M=N,
n=m,
N=k).sf(x-1)
# 4.989682834451419e-12
# R
phyper(q=x -1, m=m, n=n, k=k, lower.tail=FALSE)
# [1] 4.989683e-12
推荐阅读
- mysql - 使用带有 knex.js 的 leftjoin 仅将 where 条件应用于辅助表
- sql-server - 在 SQL Server 中使用 :out 和变量
- python - 在python中使用直方图颜色范围获取二进制图像
- mongodb - 从集合 mongoDB 中的数组中查询文档
- javascript - 如何使用 JPEG 图像生成器并将其转换为视频流?
- python - 使用字典键值对列表进行排序
- puppeteer - Puppeteer:page.goto waitUntil 选项是否包括打开的网络套接字?
- laravel - laravel base64_encode ajax 响应返回 0
- python - 寻找特定的顺序模式
- ios - 无法使用 delphi firemonkey iOS 应用程序在应用程序购买中进行测试