好多同學(xué)詢問(wèn)孟德爾隨機(jī)化的問(wèn)題，我再來(lái)嘗試著梳理一遍，希望對(duì)大家有所幫助，首先看下圖1分鐘，盯著看將下圖印在腦海中：
[圖片上傳失敗...(image-9bf17c-1668168748272)]
上圖是工具變量（不知道工具變量請(qǐng)翻之前的文章）的模式圖，明確一個(gè)點(diǎn)：我們做孟德爾的時(shí)候感興趣的是x和y的關(guān)系，也就是小b，但是我們直接去跑x對(duì)y的回歸肯定是不對(duì)的，因?yàn)橛泻芏嗟腢，因此我們借助工具變量G（關(guān)于工具變量我們之前的文章有詳細(xì)的解釋，請(qǐng)自行查閱），去估計(jì)我們感興趣的小b。
現(xiàn)在有天然良好的工具變量G，也就是我們的基因變量，此時(shí)有上面的圖，再次重申：我們感興趣的，最終希望得到準(zhǔn)確估計(jì)的值是小b，按照上圖我們應(yīng)該有GY的關(guān)系是ab，GX的關(guān)系是a，于是乎b可以寫成ab/a，就是我們感興趣的b可以換一種思路得到，如下：
[圖片上傳失敗...(image-363f74-1668168748272)]
上面的式子要跑通的話，我們需要知道G-Y的關(guān)系和G-X的關(guān)系。
但是我們GY也就是基因和結(jié)局的關(guān)系已經(jīng)有人給我們研究好了，我們可以直接去GWAS里面找研究好的summarydata拿來(lái)用就行。
但是我們的的GX也就是基因和暴露的關(guān)系也已經(jīng)有人給我們研究好了，我們可以直接去GWAS里面找研究好的summarydata拿來(lái)用就行。
也就是說(shuō)，通過(guò)孟德爾隨機(jī)化，我們完全可以毫不費(fèi)力地估計(jì)出我們需要的小b，也就是暴露和結(jié)局的關(guān)系----就是今天要再次給大家介紹的孟德爾隨機(jī)化研究。
思路就是這么清晰。就是這么清晰。搞不明白的同學(xué)再多讀幾遍。
術(shù)語(yǔ)解析
為了幫助大家理解思想，在孟德爾隨機(jī)化的實(shí)操中有幾個(gè)術(shù)語(yǔ)得提點(diǎn)一波：
連鎖不平衡（linkage disequilibrium）：剛剛講我們可以有很多的基因結(jié)局/暴露的關(guān)系的，就是GWAS里面好些基因可以用，這個(gè)時(shí)候我們不希望基因之間有相關(guān)（會(huì)造成double counting，使得結(jié)果偏倚）：
[圖片上傳失敗...(image-96faae-1668168748272)]
我們實(shí)際做的時(shí)候，模式是像上圖，snp之間你說(shuō)不相干就不相干？當(dāng)兩個(gè)位點(diǎn)的不同等位基因的關(guān)聯(lián)頻率高于或低于獨(dú)立隨機(jī)關(guān)聯(lián)的條件下的期望頻率，這種情況是客觀存在的，此時(shí)時(shí)這些工具變量之間相關(guān)性就叫連鎖不平衡，其大小可以用LD r方來(lái)表示，這個(gè)指標(biāo)也是我們?cè)诓僮鲿r(shí)需要設(shè)定的指標(biāo)之一。
水平基因多效性（Horizontal Pleiotropy）：理解這個(gè)概念先看下圖：
[圖片上傳失敗...(image-b6893b-1668168748272)]
意思是我的理想的情況是通過(guò)ab/a的操作估計(jì)出b，但是看上圖，是不是免不了會(huì)出現(xiàn)f這條路徑，如果出現(xiàn)了f，我們的基因和結(jié)局之間的關(guān)系就是f+ab，此時(shí)，我用原來(lái)的方法估計(jì)的就不是b了，而是b+f/a了，就不對(duì)了（始終記住我們關(guān)心的是b）。
但是如果我的基因變量很多，從而有很多的f，如果所有f的期望均值為0，那么最后我們匯總一下得到的結(jié)果也基本上就是b了，無(wú)傷大雅。但是就怕所有的f都是一邊偏向的（都大于0或都小于0），此時(shí)就有問(wèn)題了，叫做定向多效性directional pleiotropy，這也是為什么我們最后要做漏斗圖的原因。
就是通過(guò)漏斗圖一看都是所有的工具變量都是呈漏斗分布的，就說(shuō)明沒(méi)有偏向，這個(gè)時(shí)候我們認(rèn)為定向多效性都被沖掉了，不影響。
好，解釋了上面的一些術(shù)語(yǔ)之后，我們實(shí)操一波。
實(shí)操
最基本的例子：BMI on CHD的例子，我想看一下BMI作為暴露，CHD作為結(jié)局的mr，代碼就4條：
bmi_exp_dat <- extract_instruments(outcomes = 'ieu-a-2')
chd_out_dat <- extract_outcome_data(snps = bmi_exp_dat$SNP, outcomes = 'ieu-a-7')
dat <- harmonise_data(bmi_exp_dat, chd_out_dat)
res <- mr(dat)
結(jié)果如下，下圖中有不同方法出來(lái)的我們關(guān)心的小b：
[圖片上傳失敗...(image-38f4e7-1668168748272)]
這個(gè)就算做完了，就這么簡(jiǎn)單快速。
接下來(lái)就是敏感性分析，首先是各個(gè)工具變量的異質(zhì)性檢驗(yàn)：
mr_heterogeneity(dat)
運(yùn)行代碼后可以得到Cochran’s Q統(tǒng)計(jì)量
然后是水平基因多效性檢驗(yàn)，代碼如下：
mr_pleiotropy_test(dat)
運(yùn)行代碼可以得到egger_intercept
然后是單個(gè)SNP結(jié)果檢驗(yàn)，代碼如下：
res_single <- mr_singlesnp(dat)
運(yùn)行后可以得到每個(gè)SNP的小b
[圖片上傳失敗...(image-297c9d-1668168748272)]
然后是留一檢驗(yàn)，代碼如下：
mr_leaveoneout(dat)
接下來(lái)，論文中還會(huì)有幾個(gè)圖，首先是點(diǎn)圖，代碼如下：
mr_scatter_plot(res, dat)
[圖片上傳失敗...(image-9b736a-1668168748272)]
點(diǎn)圖是將同一個(gè)SNP對(duì)暴露的效果放在橫軸，對(duì)結(jié)局的效果放在縱軸，此時(shí)圖中的斜率就是我們的估計(jì)的小b。
然后是單個(gè)SNP效應(yīng)組合的森林圖用mr_forest_plot函數(shù)可以得到，mr_leaveoneout_plot可以得到留一分析的森林圖，mr_funnel_plot可以幫我們得到漏斗圖。
到這就出了所有需要報(bào)告的東西，做完了。
但是上面的流程有很多的前提，比如你得知道暴露和結(jié)局的GWASid才能進(jìn)行下去，GWAS又有很多，比如你直接用上面的代碼的話其實(shí)是MR Base GWAS catalog里面的GWAS，當(dāng)然你還可以選別的，或者用自己找來(lái)的最新的GWAS都是可以的。
[圖片上傳失敗...(image-6cc548-1668168748272)]
第一步首先是在相應(yīng)的GWAS中找到暴露的summary data：
那么有那些GWAS可以供我們使用呢？我們可以直接把GWAS的目錄調(diào)出來(lái)瞅瞅，代碼如下：
data(gwas_catalog)
運(yùn)行后大約可以得到15萬(wàn)個(gè)全基因組關(guān)聯(lián)研究的數(shù)據(jù)，截圖如下：
[圖片上傳失敗...(image-a8269d-1668168748272)]
那么對(duì)我們而言，我們現(xiàn)在需要找到我們關(guān)心的暴露對(duì)應(yīng)的GWAS，比如我現(xiàn)在要找與“blood”表型相關(guān)的GWAS，我可以寫出如下代碼：
exposure_gwas <- subset(gwas_catalog, grepl("Blood", Phenotype_simple)）
上面的代碼相當(dāng)于只用Phenotype_simple這一列做篩選，當(dāng)然你也可以結(jié)合其它的列比如人群，比如作者，比如地區(qū)等等，都是可以的。
選好暴露相關(guān)的GWAS之后要做的就是進(jìn)一步確定基因工具變量和暴露的強(qiáng)度，在論文中一般是這么描述：First, relevance assumption was met considering that all SNPs have reached genome-wide significance (p < 5 × 10?8)
具體的操作如下：
exposure_gwas<-exposure_gwas[exposure_gwas$pval<5*10^-8,]
通過(guò)上面的步驟保證我們的基因工具變量一定是和暴露強(qiáng)相關(guān)。
然后就是將準(zhǔn)備好的暴露的GWAS數(shù)據(jù)形成可以用來(lái)做MR分析的數(shù)據(jù)格式，需要用到format_data()函數(shù)：
exposure_data<-format_data(exposure_gwas)
此時(shí)的exposure_data大概長(zhǎng)這樣：
[圖片上傳失敗...(image-9dab47-1668168748272)]
可以看到有很多個(gè)基因工具變量SNP，這個(gè)時(shí)候我們需要考慮連鎖不平衡（linkage disequilibrium）：
exposure_data<-clump_data(exposure_data, clump_r2 = 0.001）
上面的代碼中clump_r2則是設(shè)定的容許相關(guān)性，到這兒我們算是手動(dòng)地將工具變量都篩出來(lái)了，解決了找工具變量的問(wèn)題，還有一個(gè)方法是自動(dòng)篩選工具變量，比如我暴露是bmi，我可以寫出如下代碼：
subset(ao, grepl("body mass", trait))
[圖片上傳失敗...(image-91f404-1668168748272)]
運(yùn)行后我知道我可以選的gwasid是ieu-b-40，這個(gè)時(shí)候我也可以自動(dòng)提取出工具變量，這兩種方法達(dá)到的目的都是一樣的：
extract_instruments('ieu-b-40')
然后依照工具變量進(jìn)行結(jié)局的summary estimates的提取，提取結(jié)局的summary data也需要是需要知道GWASid的對(duì)吧，比如我現(xiàn)在關(guān)心的結(jié)局是收縮壓，我就可以寫出如下代碼：
outcome_gwas <- subset(ao, grepl("Systolic", trait))
運(yùn)行后我就可以知道所有的和收縮壓相關(guān)的gwasid了，我選一個(gè)最新的，比如我就選下面的2021年的：
[圖片上傳失敗...(image-95c7a7-1668168748272)]
看圖我們知道它對(duì)于的id是ieu-b-5075,我就這么寫：
outcome_data <- extract_outcome_data(
snps = exposure_data$SNP, outcomes = "ieu-b-5075")
后續(xù)通過(guò)合并直接做mr分析就可以，流程就沒(méi)有不同了。
小結(jié)
今天給大家寫了孟德爾隨機(jī)話的實(shí)操，文章圖示例來(lái)自【中文孟德爾隨機(jī)化】英國(guó)布里斯托大學(xué)MRC-IEU《R語(yǔ)言做孟德爾隨機(jī)化》第一章：用MRBase網(wǎng)頁(yè)工具和R包TwoSampleMR做兩樣本孟德爾隨機(jī)化_嗶哩嗶哩_bilibili
