哦呦,听说过PascalX吗?这可不是什么普通的Python库,他可以用于对 GWAS 汇总统计进行高精度基因和通路评分。这小小的Python包不仅考虑了SNP之间的关系,还能玩出花样,基于χ2分布随机变量,精确到让你的双眼发直。最初的通路评分算法“Pascal”使用平方(归一化)效应大小的总和来对基因或通路进评分。但是,随着目前 GWAS 所用样本量的不断增加,信号强度也在不断增加。Pascal已经完全不够用了,以至于无法再用原始 "Pascal "软件中的双精度算法来准确计算,从而导致顶级基因之间的排名无法确定。“然后“PascalX”就登场了,像是个技术大咖,用C++多精度来搞定CDF计算,解决了那些分数疯长的问题。通过C外来函数接口,重新实现了之前Python里的原始精确基因和通路评分。而且,这货不止如此,还有一堆新的基因和通路评分方法,它允许使用我们最近提出的基因一致性和比率富集测试使用两组 GWAS 汇总统计数据进行交叉评分。它还实现了不仅通过基因转录位点周围的窗口选择SNP的可能性,而且还实现了通过其他标记数据选择SNP的可能性。- 多精度加权 cdf 计算(近似或正好最多 100 位)
- 通过自己的索引SNP数据库快速随机访问参考面板基因组数据
跑代码时卡顿、电脑不给力让人抓狂!找果叔试用稳定高速的服务器,让分析顺畅无比!
代码学不会?bug 频繁出现,束手无策?实操生信分析课程赶快学起来!滴滴果叔领取体验课程哦~
线上课程教学
课题设计、定制生信分析
生信服务器
加微信备注99领取使用
PascalX 可以通过命令行界面使用,也可以用作 python3 库。后者释放了PascalX的全部潜力。命令行界面可通过调用./pascalx.php来调用。执行 ./pascalx 命令可显示可用选项和设置。请注意,脚本具有全局设置和特定命令设置。后者可以在指定所需的全局(位置)参数和命令选择( 或 )后显示,例如:genescoringxscoring./pascalx ensemble.txt demo/EUR.simulated out.txt genescoring -h
第一个参数指定基因注释文件(如果指定文件不存在,将自动从 ensembl biomart 下载)。第二个参数设置要使用的参考面板(如果尚未导入,将自动从相应的 .vcf 文件导入)。第三个参数设置要存储结果的文件,第四个参数设置要执行的操作可以通过在文件前添加一个包含基因模块的文件来进行通路评分,例如-pw./pascalx -pw pathways.txt ensemble.txt refpanel out.txt genescoring gwasA.tsv.gz
有关在路径或 GWAS 文件中指定数据列的参数,请参考选项 -h
导入 genescorer 类并初始化基于基因评分器的总和:卡方from PascalX import genescorer
Scorer = genescorer.chi2sum()
请注意,使用 genescorer 的默认选项。具体而言,基因的起始和结束位置被延长了一个窗口,方差截止值设置为window=50000,varcutoff=0.99可通过执行以下文件夹中的脚本下载和转换 1000 基因组计划参考数据(对于 GRCh37,将 38 替换为 37)。bash get1KGGRCh38.sh pathtostore/ EUR 4 tped
第三个参数指定要使用的 CPU 内核数量。plink 转换后的文件将保存在文件夹中,文件名为 . 如果想保留所有来源的样本,请用 . 请注意,执行脚本需要 Linux 操作系统。例如,你可以在 PascalX 多平台(Docker)运行时.pathtostore/EUR.1KG.GRCh38.chr#EURALL 中运行。Scorer.load_refpanel('path/EUR.1KG.GRCh38',keepfile=None)
注意:如果尚未导入相应的参考数据,PascalX 将尝试从 .tped 文件或 .tped 文件中导入数据。对于 .tped 文件,基因型信息必须以 gzip 压缩的 1-2 编码 plink tped 文件形式提供。以下 plink 选项可以完成这项工作:. 默认情况下,PascalX 只使用一个 CPU 核心进行导入。你可以通过设置选项.filename.chr#.tped.gzfilename.chr#.vcf.gzpath/---recode 12 transposeparallel= 来增加使用的核心数量。为了将等位基因信息导入参照面板,必须使用原始 .vcf 文件进行导入。为此,请将上述导入脚本选项替换为 .vcf 文件。需要注意的是,如果只想在导入 .vcf 时保留一部分样本,则必须设置.tpedvcfkeepfile= 选项。如果要使用 ensembl 基因注释,可从 BioMart 自动下载,如下所示。from PascalX.genome import genome
G = genome()
G.get_ensembl_annotation('your_filename.tsv',genetype='protein_coding, lncRNA',version='GRCh38')
在该选项中,所有有效的 ensembl 基因类型都可以逗号分隔的字符串形式提供。PascalX 可通过加载的参照面板中的位置重叠将基因与变异 id 匹配。两个数据集必须基于相同的注释版本(例如都是 GRCh38)!Scorer.load_genome('path/filename')
也可以加载自己的注释。请参阅 Genescorer API 文档,了解需要设置的选项。可以通过 method.load_mapping 方法直接使用 SNP 到基因的映射,而不是使用位置注释。Scorer.load_GWAS('path/filename',rscol=0,pcol=1,a1col=None,a2col=None,header=False)
参数指定变异 ID 栏和 p 值栏。这两列分别包含候补等位基因和参考等位基因。如果参照面板不包含等位基因信息(.tped 导入),则将这两列都设为 。指定第一行是否为标题。文件可以是原始文本,也可以是 gzip 压缩文件,文件名以 .rscol=pcol=a1col=a2col=Noneheader=.gz 结尾。 请注意,加载的 GWAS SNP 可通过以下命令可视化:Scorer.plot_genesnps('AOAH',mark_window=True,show_correlation=True);
按上述方法加载基因注释、参考面板和 GWAS 后,可按如下方法计算基因得分。#对加载注释中的所有基因进行评分
R = Scorer.score_all()
#只给 21 号和 22 号染色体上的基因打分
R = Scorer.score_chr(chrs=[21,22])
#对基因 WDR12 和 FARP2 进行评分
R = Scorer.score(['WDR12','FARP2'])
R = [R_SUCCESS,R_FAIL,R_TOTALFAIL]
R_SUCCESS = [ ['Symbol',p-value,NSNP],...]
R_FAIL = [ ['Symbol',[infos]] ,...]
R_TOTALFAIL = [ ['Symbol','Reason'] ,...]
这里是得分成功的基因列表、因得分算法不收敛而得分失败的基因列表以及因其他原因(如无可用 SNPs)而得分失败的基因列表。默认情况下,genescorer 使用鞍点近似法计算 CDF()。如需精确计算,建议使用()自动选择最合适的算法和精度。可以使用皮尔逊算法自动重新评分一次()。请注意,如果迭代次数足够多,Ruben 在最大精度下最终会收敛。不过,如果最大特征值和最小特征值之间的比率很大,收敛速度可能会很慢。在这种情况下,通过基因分析器的参数来减少保持的方差往往会有所帮助。请注意,可以使用 Genescorer.method='saddle'method='auto'R_FAILautorescore=Truevarcutoff=R 对结果进行手动重新评分。Scorer.save_scores('filename')
Scorer.plot_Manhattan(R[0])
PascalX 提供两种不同的路径评分器。如上所述,两者都需要一个完全初始化的基因扫描器。定义基因评分器,为 GWAS 评分或加载评分基因。请注意,已保存的基因评分可通过以下方式读入:Scorer.load_scores('filename')
路径评分器的初始化过程如下。首先通过以下方式导入路径评分类from PascalX import pathway
Pscorer = pathway.chi2rank(Scorer)
等级评分器通过等级统一基因 p 值分布,并通过反变换将 p 值汇总为分布式随机变量。Pscorer = pathway.chi2perm(Scorer)
由于计算成本的原因,随机生成的基因集中的基因不会被融合。一般来说,我们建议只使用基于等级的评分方法。 基因模块/通路集可通过以下命令从以制表符分隔的文件中加载M = Pscorer.load_modules('filename.tsv',ncol=0,fcol=2)
ncol= 是包含模块名称的一列,也是包含基因符号的第一列。假定其他成员基因在随后的列中。 RESULT = Pscorer.score(M)
RESULT = [ ['name',[genes],[gene p-values],p-value],...]
from PascalX import xscorer
X = xscorer.zsum(leftTail=False)
X.load_genome('path/filename')
请注意,上述基因评分的默认初始化设置为测试侧,对应于反一致性和一致性。与标准基因评分器一样,必须加载基因注释。X.load_GWAS('path/filenameA',name='GWAS A',rscol=0,pcol=1,bcol=2,a1col=None,a2col=None,header=False)
X.load_GWAS('path/filenameB',name='GWAS B',rscol=0,pcol=1,bcol=2,a1col=None,a2col=None,header=False)
在 GWAS 数据加载例程中,我们还必须通过参数为每个要加载的 GWAS 设置一个名称,并指定包含原始贝塔值的列。X.matchAlleles('GWAS A','GWAS B')
X.jointlyRank('GWAS A','GWAS B')
R = X.score_all(E_A='GWAS A',E_B='GWAS B')
X.plot_genesnps('TRIM26','GWAS A','GWAS B',mark_window=True,show_correlation=True);
X = xscorer.rsum(leftTail=False)
from PascalX.genexpr import genexpr
GE = genexpr()
GE.load_genome('yourannotationfile')
GE.get_GTEX_expr('yourGTEXfilename')
GE.load_expr('GTEX/yourGTEXfilename')
组织表达(在 TPM 中)可以通过以下方式可视化基因列表GE.plot_genexpr(['AOAH','CSMD1','CDH13'],tzscore=True,cbar_pos=(0.0, 0.1, 0.01, 0.5))
PascalX 以与通路评分类似的方式测试组织富集。富集测试是在基因列表上进行的。将邻近的基因与元基因融合,并从原始 GTEx 读取计数中重新计算元基因的 TPM 值。基因TPM值通过排序进行均匀化,并通过逆cdf转换为分布式随机变量。根据具有 (# 个基因) 自由度的分布来检验总和。要测试通路的富集,请提供基因成员。要检测GWAS富集,请提供重要基因列表。R = GE.chi2rank([ ['PathwayName',['AOAH','CSMD1',' CDH13'] ] ])
R = ({'PathwayName': {'Tissue': pval,... },... }, FAILS, GENES, [{'Tissue': [pvalg1,...]}] )
PascalX 支持外部 SNP 基因碱基信息,例如 eQTL 或染色质相互作用数据集,以补充默认的基于距离的 SNP 基因关联。这需要包含 rsid(变体 ID)和与其关联的基因 Ensembl ID 的外部数据集。首先,对于标准基因评分,必须初始化基因评分器,并加载参考panel、ensembl基因注释和GWAS SNP。from PascalX import genescorer
Scorer = genescorer.chi2sum()
Scorer.load_refpanel("path/filename")
Scorer.load_genome("path/filename")
Scorer.load_GWAS("path/filename",rscol=0,pcol=1,header=False)
Scorer.load_mapping("path/filename", rcol=0, gcol=1, delimiter=',', header=True, joint=True)
from PascalX import pathway
Pscorer = pathway.chi2rank(Scorer)
M = Pscorer.load_modules("filename.tsv",ncol=0,fcol=2)
RESULT = Pscorer.score(M)
“PascalX”,通过C++多精度实现CDF计算解决了这个问题,通过C外来函数接口接口,重新实现了Lamparter等人在Python中的原始精确基因和通路评分。与原来的“Pascal”相比, “PascalX”实现提供了一些进一步的技术改进,即内在并行化、线性代数运算的 GPU 加速可能性以及具有快速随机访问的新数据存储模型。此外,“PascalX”还包括几种新的基因和通路评分方法。具体来说,它允许使用我们最近提出的基因一致性和比率富集测试使用两组 GWAS 汇总统计数据进行交叉评分。它还实现了不仅通过基因转录位点周围的窗口选择SNP的可能性,而且还实现了通过其他标记数据选择SNP的可能性。最后,大家如果对生信分析感兴趣但还不熟悉,又想尝试一下处理自己的数据,不妨试一下果叔开发的生信云平台哦,一键出图,一键导出CNS级别的Figture!!赶快去试试吧!!点击 http://www.biocloudservice.com/home.html不会分析还想用生信工具助力发文咋办?有这顾虑的朋友,想一步到位就带着想法来,不论是代码实操还是在线文章结果复现,果叔照样能提供,还有大家都想要的服务器,找果叔获取就对了!1.PascalX: a Python library for GWAS gene and pathway enrichment tests(https://doi.org/10.1093/bioinformatics/btad296)2.https://github.com/BergmannLab/PascalX3.https://bergmannlab.github.io/PascalX/