标签:多表 传统 keep type ade err 个人 rename http
欢迎来到"bio生物信息"的世界
新年的第一篇更文。
祝大家新春快乐!身体健康!
18号回家以后,经历了如下过程。
20号 喉咙痛
21号 喉咙痛
22号喉咙痛 咳嗽
23-24号 咳嗽
25号 咳嗽为主 鼻塞 夜间咳嗽加剧
26号 咳嗽为主 鼻塞 流鼻涕 夜间咳嗽加剧
27号 咳嗽为主 鼻塞 流鼻涕 打喷嚏 夜间咳嗽加剧
28号 咳嗽为主 鼻塞 流鼻涕 打喷嚏 夜间咳嗽加剧
换了多种药物,均不见效。
索性也少关注新型病毒肺炎的事了,缓解些焦虑感。
好了,以上是交代我这段时间没更文的原因。
主要是懒,其次是生病。
这次给大家介绍一个R包metaCCA
传统的GWAS分析都是基于多个SNP与单表型的关联分析,找出显著的候选位点。
随着复杂表型的不断涌现,这种分析慢慢的有些局限。
比如,已有多个研究表明,精神分裂症、抑郁症、自闭症、双向情感障碍等疾病存在着遗传上的相关性,也就是说,他们有着共享的风险基因,但这些共享的风险基因是什么?是一个位点还是多个位点?一个基因座还是多个基因座?共享基因是如何共同影响多种不同的疾病?
这些问题都没有很好的被解释清楚。
因此,通过metaCCA就能解决这个问题。
其功能有两个。
第一、分析单个SNP与多种表型的相关性;
第二、分析多个SNP与多种表型的相关性;
metaCCA的思想类似于降维。
将多维的基因型数据和多维的表型数据各自进行降维,再计算降维之后的基因型数据和表型数据之间的相关性。
metaCCA提出了两个算法:metaCCA和metaCCA+。
这两个算法的适用范围不大一样。
诸位进行分析前,需要对自己的数据有一定的了解,再酌情选择合适的算法。
适用范围不一样表现在:
metaCCA算法适用于基因型相关系数矩阵XX和表型数据相关系数矩阵YY的计算结果比较精确的情况。举例来说,基因型相关系数矩阵是通过GWAS summary数据对应参考人群基因组评估出来的。
metaCCA+算法适用于基因型相关系数矩阵XX和表型数据相关系数矩阵YY计算结果不是很准的情况下。举例来说,参考人群的基因型数据只有1000个位点,样本数只有10个人,而GWAS summary数据就有几百万个位点,几十万个样本,在这种情况下,GWAS summary数据通过参考人群基因组评估出来的基因型相关系数矩阵XX和表型数据相关系数矩阵YY就不是很准确。此时强行使用metaCCA算法就容易有假阳性结果。这种情况下,最好是用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‘
解决办法:先番蔷(自己意会)。再重新安装就可以了。
我没有尝试过在不借助工具的情况下,解决这个报错,如果你有办法,欢迎告知,谢谢。
进行单个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位点;
计算基因型相关系数矩阵前需要准备几个数据,第一个是下载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
即为基因型的相关系数矩阵;
准备好输入文件S_XX_study
和 S_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