码迷,mamicode.com
首页 > Web开发 > 详细

使用metaCCA进行单/多个SNP与多表型的典型相关性分析

时间:2020-01-28 17:46:42      阅读:192      评论:0      收藏:0      [点我收藏+]

标签:多表   传统   keep   type   ade   err   个人   rename   http   

欢迎来到"bio生物信息"的世界

新年的第一篇更文。

祝大家新春快乐!身体健康!

18号回家以后,经历了如下过程。

20号 喉咙痛

21号 喉咙痛

22号喉咙痛 咳嗽

23-24号 咳嗽

25号 咳嗽为主 鼻塞 夜间咳嗽加剧

26号 咳嗽为主 鼻塞 流鼻涕 夜间咳嗽加剧

27号 咳嗽为主 鼻塞 流鼻涕 打喷嚏 夜间咳嗽加剧

28号 咳嗽为主 鼻塞 流鼻涕 打喷嚏 夜间咳嗽加剧

换了多种药物,均不见效。

索性也少关注新型病毒肺炎的事了,缓解些焦虑感。

好了,以上是交代我这段时间没更文的原因。

主要是懒,其次是生病。


这次给大家介绍一个R包metaCCA

1 metaCCA有什么用

传统的GWAS分析都是基于多个SNP与单表型的关联分析,找出显著的候选位点。

随着复杂表型的不断涌现,这种分析慢慢的有些局限。

比如,已有多个研究表明,精神分裂症、抑郁症、自闭症、双向情感障碍等疾病存在着遗传上的相关性,也就是说,他们有着共享的风险基因,但这些共享的风险基因是什么?是一个位点还是多个位点?一个基因座还是多个基因座?共享基因是如何共同影响多种不同的疾病?

这些问题都没有很好的被解释清楚。

因此,通过metaCCA就能解决这个问题。

其功能有两个。

第一、分析单个SNP与多种表型的相关性;

第二、分析多个SNP与多种表型的相关性;

2 metaCCA思想

metaCCA的思想类似于降维。

将多维的基因型数据和多维的表型数据各自进行降维,再计算降维之后的基因型数据和表型数据之间的相关性。

metaCCA提出了两个算法:metaCCA和metaCCA+。

这两个算法的适用范围不大一样。

诸位进行分析前,需要对自己的数据有一定的了解,再酌情选择合适的算法。

适用范围不一样表现在:

metaCCA算法适用于基因型相关系数矩阵XX和表型数据相关系数矩阵YY的计算结果比较精确的情况。举例来说,基因型相关系数矩阵是通过GWAS summary数据对应参考人群基因组评估出来的。

metaCCA+算法适用于基因型相关系数矩阵XX和表型数据相关系数矩阵YY计算结果不是很准的情况下。举例来说,参考人群的基因型数据只有1000个位点,样本数只有10个人,而GWAS summary数据就有几百万个位点,几十万个样本,在这种情况下,GWAS summary数据通过参考人群基因组评估出来的基因型相关系数矩阵XX和表型数据相关系数矩阵YY就不是很准确。此时强行使用metaCCA算法就容易有假阳性结果。这种情况下,最好是用metaCCA+算法。

3 metaCCA安装

if (!requireNamespace("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install("metaCCA")

这一步如果产生如下错误:

Error in read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
无法打开链结
此外: Warning messages:
1: In download.file(url, destfile, method, mode = "wb", ...) :
downloaded length 409600 != reported length 1002514
2: In unzip(zipname, exdir = dest) : 从zip文件中抽取1时出了错
3: In read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
无法打开压缩文件‘metaCCA/DESCRIPTION‘,可能是因为‘No such file or directory‘

解决办法:先番蔷(自己意会)。再重新安装就可以了。

我没有尝试过在不借助工具的情况下,解决这个报错,如果你有办法,欢迎告知,谢谢。

4 基于单个SNP与多种表型的相关性分析

进行单个SNP与多表型的相关性分析,需要准备一个输入文件,这里将输入文件命名微S_XY_full_study

S_XY_full_study输入文件的格式如下:

技术图片

每一个位点占一行,文件列分别微位点的两个基因型(allele_0、allele_1)、不同表型在该位点的效应值(trait1_b,trait2_b,……)和误差(trait1_se, trait2_se, ……)。

注意,文件顺序是位点、位点基因型1、位点基因型2、表型1效应值、表型1效应值误差、表型2效应值、表型2效应值误差、表型3效应值、表型3效应值误差、……、表型n效应值、表型n效应值误差;

基因型由“A ”,“C”, “G”,“T”组成;

准备好输入文件好后,输入以下命令,即可得到单个SNP与多表型的相关性结果

metaCCA_res1 =metaCcaGp( nr_studies = 1,S_XY =list( S_XY_full_study ),std_info = 0,S_YY =list( estimateSyy( S_XY_full_study ) ),N = N )

N指的是GWAS summary的样本数;

生成的结果如下所示:

技术图片

总共有三列。

第一列为每个位点与多种表型的典型相关性值;

第二列为P值;

第三列为每个表型的权重值;

P值的显著性阈值为0.05/SNP数量;

如果只想一个位点一个位点放进去计算,则需要加上analysis_type = 1 和 SNP_id的参数,命令如下:

metaCCA_res1 =metaCcaGp( nr_studies = 1,S_XY =list( S_XY_full_study ),std_info = 0,S_YY =list( estimateSyy( S_XY_full_study ) ),N = N,analysis_type = 1,SNP_id =‘rs666‘) )

rs666指的是想单独计算的SNP位点;

5 多个SNP与多种表型的典型相关性分析

5.1 计算基因型相关系数矩阵

计算基因型相关系数矩阵前需要准备几个数据,第一个是下载GWAS summary对应人群的参考基因组数据,假如是中国人群,则需要下载中国人群的公共数据库,这里以千人基因组为例。

计算基因型相关系数矩阵需要有PLINK软件的使用基础。

1000G为千人基因组包含个体层次的基因型数据;

SNP为一个文件,包含GWAS summary文件的SNP ID号,通常以rs开头,每一个SNP为一行;

CHB_sample为一个文件,包含中国人群的样本ID,文件有两列,分别为FID和IID,详细参考PLINK的输入文件格式;

S_XX_study为输出文件的文件名;

准备好后,输入如下命令:

plink2 --bfile 1000G --extract SNP --keep CHB_sample --r2 inter-chr with-freqs --ld-window-r2 0 --make-bed --out S_XX_study

输出文件S_XX_study即为基因型的相关系数矩阵;

5.2 计算多个SNP与多种表型的典型相关性

准备好输入文件S_XX_studyS_XY_full_study后,假定想研究基因A上的五个SNP位点‘rs10‘,‘rs80‘,‘rs140‘,‘rs170‘,‘rs172‘与多种表型的典型相关结果,则输入以下命令:

metaCCA_res3 =metaCcaGp( nr_studies = 1,S_XY =list( S_XY_full_study ),std_info = 0 ,S_YY =list( estimateSyy( S_XY_full_study ) ),N = N,analysis_type = 2,SNP_id =c(‘rs10‘,‘rs80‘,‘rs140‘,‘rs170‘,‘rs172‘),S_XX =list( S_XX_study))

N指的是GWAS summary的样本数;

生成的结果如下所示:

技术图片

总共有四列。

第一列为‘rs10‘,‘rs80‘,‘rs140‘,‘rs170‘,‘rs172‘与多种表型的典型相关性值;

第二列为P值;

第三列为每个表型的权重值;

第四列为每个SNP位点的权重值;

P值的显著性阈值为0.05/metaCCA_res3行数;

好了,今天的内容就讲到这,感兴趣的看原文文献:
https://academic.oup.com/bioinformatics/article/32/13/1981/1742730

祝各位假期愉快!学业有成。

使用metaCCA进行单/多个SNP与多表型的典型相关性分析

标签:多表   传统   keep   type   ade   err   个人   rename   http   

原文地址:https://www.cnblogs.com/chenwenyan/p/12238275.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!