玩转GWAS基因和通路富集评分,PascalX让科研充满趣味!

学术   2024-11-16 19:00   上海  
前言
哦呦,听说过PascalX吗?这可不是什么普通的Python库,他可以用于对 GWAS 汇总统计进行高精度基因和通路评分。这小小的Python包不仅考虑了SNP之间的关系,还能玩出花样,基于χ2分布随机变量,精确到让你的双眼发直。最初的通路评分算法“Pascal”使用平方(归一化)效应大小的总和来对基因或通路进评分。但是,随着目前 GWAS 所用样本量的不断增加,信号强度也在不断增加。Pascal已经完全不够用了,以至于无法再用原始 "Pascal "软件中的双精度算法来准确计算,从而导致顶级基因之间的排名无法确定。“然后“PascalX”就登场了,像是个技术大咖,用C++多精度来搞定CDF计算,解决了那些分数疯长的问题。通过C外来函数接口,重新实现了之前Python里的原始精确基因和通路评分。而且,这货不止如此,还有一堆新的基因和通路评分方法,它允许使用我们最近提出的基因一致性和比率富集测试使用两组 GWAS 汇总统计数据进行交叉评分。它还实现了不仅通过基因转录位点周围的窗口选择SNP的可能性,而且还实现了通过其他标记数据选择SNP的可能性。
PascalX在富集评分上的优势:
  • 多精度加权 cdf 计算(近似或正好最多 100 位)
  • 染色体和/或基因的平行化
  • GPU 支持线性代数运算加速
  • 通过自己的索引SNP数据库快速随机访问参考面板基因组数据
  • 两个GWAS之间的基因一致性测试
  • 支持通过外部数据选择SNP    
  • 组织富集试验(实验性)

跑代码时卡顿、电脑不给力让人抓狂!找果叔试用稳定高速的服务器,让分析顺畅无比!

代码学不会?bug 频繁出现,束手无策?实操生信分析课程赶快学起来!滴滴果叔领取体验课程哦~


线上课程教学

课题设计、定制生信分析

生信服务器

加微信备注99领取使用



用法
PascalX 可以通过命令行界面使用,也可以用作 python3 库。后者释放了PascalX的全部潜力。
命令行界面可通过调用./pascalx.php来调用。执行 ./pascalx 命令可显示可用选项和设置。
./pascalx -h
请注意,脚本具有全局设置和特定命令设置。后者可以在指定所需的全局(位置)参数和命令选择( 或 )后显示,例如: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 genescorerScorer = 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 genomeG = 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 到基因的映射,而不是使用位置注释。
GWAS 统计摘要
要评分的 GWAS 统计摘要可通过以下方式导入
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 pathwayPscorer = 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 xscorerX = 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')
必须导入GTEx数据。对于自动导入,请调用
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,...]}] )         

SNP基因评分
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)
导入外部 SNP 基因数据集:
Scorer.load_mapping("path/filename", rcol=0, gcol=1, delimiter=',', header=True, joint=True)
基因评分    
R = Scorer.score_all()
通路评分
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
不会分析还想用生信工具助力发文咋办?有这顾虑的朋友,想一步到位就带着想法来,不论是代码实操还是在线文章结果复现,果叔照样能提供,还有大家都想要的服务器,找果叔获取就对了!

往期回顾

01

顶刊Cell的富贵还是让“多组学”接住了!樊荣团队开发新技术,天花板的代码分享!携单细胞测序玩转空间转录组,别错过!

02

刚逃过“三花淡奶”又陷入“预制菜”!“NHANES最新数据+多变量回归”实锤:超加工食品会加重骨质疏松!不爱做饭的亲,注意咯!

03

遇强则强!还是低估孟德尔随机化的实力了!山东第一医科大学团队搭配机器学习+多组学分析,便轻轻松松0实验发了6分+!

04

思路简直开挂,已经复现坐等毕业了!青岛大学的网络药理学研究咋就这么牛?看来“药” 做就做高级范,毕业才能so easy!

参考文献
1.PascalX: a Python library for GWAS gene and pathway enrichment tests(https://doi.org/10.1093/bioinformatics/btad296)
2.https://github.com/BergmannLab/PascalX
3.https://bergmannlab.github.io/PascalX/ 

生信果
生信入门、R语言、生信图解读与绘制、软件操作、代码复现、生信硬核知识技能、服务器等
 最新文章