歡迎來到"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_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
祝各位假期愉快！學業有成。