歡迎來(lái)到"bio生物信息"的世界
新年的第一篇更文。
祝大家新春快樂(lè)！身體健康！
18號(hào)回家以后，經(jīng)歷了如下過(guò)程。
20號(hào) 喉嚨痛
21號(hào) 喉嚨痛
22號(hào)喉嚨痛 咳嗽
23-24號(hào) 咳嗽
25號(hào) 咳嗽為主 鼻塞 夜間咳嗽加劇
26號(hào) 咳嗽為主 鼻塞 流鼻涕 夜間咳嗽加劇
27號(hào) 咳嗽為主 鼻塞 流鼻涕 打噴嚏 夜間咳嗽加劇
28號(hào) 咳嗽為主 鼻塞 流鼻涕 打噴嚏 夜間咳嗽加劇
換了多種藥物，均不見(jiàn)效。
索性也少關(guān)注新型病毒肺炎的事了，緩解些焦慮感。
好了，以上是交代我這段時(shí)間沒(méi)更文的原因。
主要是懶，其次是生病。
這次給大家介紹一個(gè)R包metaCCA
1 metaCCA有什么用
傳統(tǒng)的GWAS分析都是基于多個(gè)SNP與單表型的關(guān)聯(lián)分析，找出顯著的候選位點(diǎn)。
隨著復(fù)雜表型的不斷涌現(xiàn)，這種分析慢慢的有些局限。
比如，已有多個(gè)研究表明，精神分裂癥、抑郁癥、自閉癥、雙向情感障礙等疾病存在著遺傳上的相關(guān)性，也就是說(shuō)，他們有著共享的風(fēng)險(xiǎn)基因，但這些共享的風(fēng)險(xiǎn)基因是什么？是一個(gè)位點(diǎn)還是多個(gè)位點(diǎn)？一個(gè)基因座還是多個(gè)基因座？共享基因是如何共同影響多種不同的疾病？
這些問(wèn)題都沒(méi)有很好的被解釋清楚。
因此，通過(guò)metaCCA就能解決這個(gè)問(wèn)題。
其功能有兩個(gè)。
第一、分析單個(gè)SNP與多種表型的相關(guān)性；
第二、分析多個(gè)SNP與多種表型的相關(guān)性；
2 metaCCA思想
metaCCA的思想類(lèi)似于降維。
將多維的基因型數(shù)據(jù)和多維的表型數(shù)據(jù)各自進(jìn)行降維，再計(jì)算降維之后的基因型數(shù)據(jù)和表型數(shù)據(jù)之間的相關(guān)性。
metaCCA提出了兩個(gè)算法：metaCCA和metaCCA+。
這兩個(gè)算法的適用范圍不大一樣。
諸位進(jìn)行分析前，需要對(duì)自己的數(shù)據(jù)有一定的了解，再酌情選擇合適的算法。
適用范圍不一樣表現(xiàn)在：
metaCCA算法適用于基因型相關(guān)系數(shù)矩陣XX和表型數(shù)據(jù)相關(guān)系數(shù)矩陣YY的計(jì)算結(jié)果比較精確的情況。舉例來(lái)說(shuō)，基因型相關(guān)系數(shù)矩陣是通過(guò)GWAS summary數(shù)據(jù)對(duì)應(yīng)參考人群基因組評(píng)估出來(lái)的。
metaCCA+算法適用于基因型相關(guān)系數(shù)矩陣XX和表型數(shù)據(jù)相關(guān)系數(shù)矩陣YY計(jì)算結(jié)果不是很準(zhǔn)的情況下。舉例來(lái)說(shuō)，參考人群的基因型數(shù)據(jù)只有1000個(gè)位點(diǎn)，樣本數(shù)只有10個(gè)人，而GWAS summary數(shù)據(jù)就有幾百萬(wàn)個(gè)位點(diǎn)，幾十萬(wàn)個(gè)樣本，在這種情況下，GWAS summary數(shù)據(jù)通過(guò)參考人群基因組評(píng)估出來(lái)的基因型相關(guān)系數(shù)矩陣XX和表型數(shù)據(jù)相關(guān)系數(shù)矩陣YY就不是很準(zhǔn)確。此時(shí)強(qiáng)行使用metaCCA算法就容易有假陽(yáng)性結(jié)果。這種情況下，最好是用metaCCA+算法。
3 metaCCA安裝
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("metaCCA")
這一步如果產(chǎn)生如下錯(cuò)誤：
Error in read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
無(wú)法打開(kāi)鏈結(jié)
此外: 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時(shí)出了錯(cuò)
3: In read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
無(wú)法打開(kāi)壓縮文件'metaCCA/DESCRIPTION'，可能是因?yàn)?No such file or directory'
解決辦法：先番薔（自己意會(huì)）。再重新安裝就可以了。
我沒(méi)有嘗試過(guò)在不借助工具的情況下，解決這個(gè)報(bào)錯(cuò)，如果你有辦法，歡迎告知，謝謝。
4 基于單個(gè)SNP與多種表型的相關(guān)性分析
進(jìn)行單個(gè)SNP與多表型的相關(guān)性分析，需要準(zhǔn)備一個(gè)輸入文件，這里將輸入文件命名微S_XY_full_study
S_XY_full_study輸入文件的格式如下：
每一個(gè)位點(diǎn)占一行，文件列分別微位點(diǎn)的兩個(gè)基因型（allele_0、allele_1）、不同表型在該位點(diǎn)的效應(yīng)值（trait1_b,trait2_b，……）和誤差(trait1_se, trait2_se, ……)。
注意，文件順序是位點(diǎn)、位點(diǎn)基因型1、位點(diǎn)基因型2、表型1效應(yīng)值、表型1效應(yīng)值誤差、表型2效應(yīng)值、表型2效應(yīng)值誤差、表型3效應(yīng)值、表型3效應(yīng)值誤差、……、表型n效應(yīng)值、表型n效應(yīng)值誤差;
基因型由“A ”，“C”, “G”，“T”組成；
準(zhǔn)備好輸入文件好后，輸入以下命令，即可得到單個(gè)SNP與多表型的相關(guān)性結(jié)果
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的樣本數(shù)；
生成的結(jié)果如下所示：
總共有三列。
第一列為每個(gè)位點(diǎn)與多種表型的典型相關(guān)性值；
第二列為P值；
第三列為每個(gè)表型的權(quán)重值；
P值的顯著性閾值為0.05/SNP數(shù)量；
如果只想一個(gè)位點(diǎn)一個(gè)位點(diǎn)放進(jìn)去計(jì)算，則需要加上analysis_type = 1 和 SNP_id的參數(shù)，命令如下：
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指的是想單獨(dú)計(jì)算的SNP位點(diǎn)；
5 多個(gè)SNP與多種表型的典型相關(guān)性分析
5.1 計(jì)算基因型相關(guān)系數(shù)矩陣
計(jì)算基因型相關(guān)系數(shù)矩陣前需要準(zhǔn)備幾個(gè)數(shù)據(jù)，第一個(gè)是下載GWAS summary對(duì)應(yīng)人群的參考基因組數(shù)據(jù)，假如是中國(guó)人群，則需要下載中國(guó)人群的公共數(shù)據(jù)庫(kù)，這里以千人基因組為例。
計(jì)算基因型相關(guān)系數(shù)矩陣需要有PLINK軟件的使用基礎(chǔ)。
1000G為千人基因組包含個(gè)體層次的基因型數(shù)據(jù)；
SNP為一個(gè)文件，包含GWAS summary文件的SNP ID號(hào)，通常以rs開(kāi)頭，每一個(gè)SNP為一行；
CHB_sample為一個(gè)文件，包含中國(guó)人群的樣本ID，文件有兩列，分別為FID和IID，詳細(xì)參考PLINK的輸入文件格式；
S_XX_study為輸出文件的文件名；
準(zhǔn)備好后，輸入如下命令：
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即為基因型的相關(guān)系數(shù)矩陣；
5.2 計(jì)算多個(gè)SNP與多種表型的典型相關(guān)性
準(zhǔn)備好輸入文件S_XX_study 和 S_XY_full_study后，假定想研究基因A上的五個(gè)SNP位點(diǎn)'rs10','rs80','rs140','rs170','rs172'與多種表型的典型相關(guān)結(jié)果，則輸入以下命令：
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的樣本數(shù)；
生成的結(jié)果如下所示：
總共有四列。
第一列為
'rs10','rs80','rs140','rs170','rs172'與多種表型的典型相關(guān)性值；第二列為P值；
第三列為每個(gè)表型的權(quán)重值；
第四列為每個(gè)SNP位點(diǎn)的權(quán)重值；
P值的顯著性閾值為0.05/metaCCA_res3行數(shù)；
好了，今天的內(nèi)容就講到這，感興趣的看原文文獻(xiàn)：
https://academic.oup.com/bioinformatics/article/32/13/1981/1742730
祝各位假期愉快！學(xué)業(yè)有成。
