欢迎来到"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输入文件的格式如下：
image
每一个位点占一行，文件列分别微位点的两个基因型（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的样本数；
生成的结果如下所示：
image
总共有三列。
第一列为每个位点与多种表型的典型相关性值；
第二列为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_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的样本数；
生成的结果如下所示：
image
总共有四列。
第一列为
'rs10','rs80','rs140','rs170','rs172'与多种表型的典型相关性值；第二列为P值；
第三列为每个表型的权重值；
第四列为每个SNP位点的权重值；
P值的显著性阈值为0.05/metaCCA_res3行数；
好了，今天的内容就讲到这，感兴趣的看原文文献：
https://academic.oup.com/bioinformatics/article/32/13/1981/1742730
祝各位假期愉快！学业有成。
网友评论