GWAS - PRS多基因風(fēng)險(xiǎn)評分計(jì)算學(xué)習(xí)筆記

一、安裝PRSice（mac版）

經(jīng)試驗(yàn)我覺得直接從git hub中下載對應(yīng)的安裝包是最快的：https://github.com/choishingwan/PRSice，下載之后解壓，解壓文件如圖所示

  • PRSice_mac為軟件的運(yùn)行腳本
  • PRSice.R是執(zhí)行腳本的封裝
  • TOY開頭的是數(shù)據(jù)集，分為BASE和TARGET兩部分
    *接下來可以先用這個數(shù)據(jù)集測試一下

跟plink的安裝一樣，打開terminal，cd到解壓的文件夾，./PRSice_mac如果顯示下面的畫面表現(xiàn)為安裝成功：



之后下載R包

Rscript PRSice.R --dir .

Rscript 表示用R腳本來調(diào)用該軟件
dir參數(shù)為指定R包ggplot2安裝的路徑
因?yàn)榻Y(jié)果展示會調(diào)用ggplot2畫圖進(jìn)行可視化，如果你的R中已經(jīng)安裝了這個包，這個參數(shù)可以不要!!!!(我的R中有，所以后文中的此條命令我均不需要運(yùn)行)
**因?yàn)橹苯舆\(yùn)行該命令總?cè)菀讏?bào)錯，因此也可以通過R studio來安裝這個包：點(diǎn)擊面板中的package，選擇ggplot2點(diǎn)擊install，如果太慢的話選擇清華的鏡像就會快一點(diǎn)

二、文件準(zhǔn)備

上面安裝成功的圖中可以看到PRSice的注釋信息，需要準(zhǔn)備base和target兩個文件

  • Base Dataset：以空格分割的SNP關(guān)聯(lián)分析文件，為.assoc結(jié)尾的。關(guān)聯(lián)分析的教程見前面。
    必須要有的幾列是：SNP, A1, OR or BETA（對于分類型表型提供OR值）, P，其他幾列不是必須要的，如下圖：


    BASE

    如果表頭沒有以上述名稱命名，必須用參數(shù)--chr, --A1, --A2, --stat, --snp, --bp, --se, --pvalue指明表頭。
    參數(shù)可以這樣設(shè)置：--snp SNP --chr CHR --bp BP --A1 A1 --A2 A2 --stat OR --se SE --pvalue P
    或者：--snp 0 --chr 1 --bp 2 --A1 3 --A2 4 --stat 5 --se 6 --pvalue 7 --index

  • Target Dataset 是plink產(chǎn)生的二進(jìn)制文件.bed, .bim, 和 .fam file，要經(jīng)過質(zhì)控QC，過程翻看之前的教程。

三、PRSice運(yùn)行

#二元性狀用的是OR值，命令如下
Rscript PRSice.R --dir . \ 
  --prsice ./PRSice \
  --base TOY_BASE_GWAS.assoc \
  --target TOY_TARGET_DATA \
  --thread 1 \
  --stat OR \
  --binary-target T

#數(shù)量性狀用的是BETA值，命令如下
Rscript PRSice.R --dir . \
  --prsice ./PRSice \
  --base TOY_BASE_GWAS.assoc \
  --target TOY_TARGET_DATA \
  --thread 1 \
  --stat BETA \
  --beta \
  --binary-target F
  • --dir參數(shù)指定可執(zhí)行文件的路徑；（安裝過ggplot2的可以直接忽略掉這一條命令）
  • base參數(shù)指定base data（驗(yàn)證樣本）的關(guān)聯(lián)分析結(jié)果，包含每個遺傳變異的效應(yīng)值和 p 值。
  • target參數(shù)指定target data（測試樣本）的分型結(jié)果，指 target 樣本二進(jìn)制 plink 格式的基因型數(shù)據(jù)前綴，也支持 BGEN 格式
  • thread參數(shù)指定線程數(shù)
  • stat指定base data關(guān)聯(lián)分析結(jié)果中的效應(yīng)值，有OR和BETA兩種取值
  • binary-target參數(shù)指定target data中分型結(jié)果是否為二分類的表型，T=ture
    ·
  • 如果下載下來的數(shù)據(jù)第六行并不是按我們的需求寫上是表型的，為了后面的關(guān)聯(lián)分析，要手工定義第六行是實(shí)驗(yàn)組，即需要把.fam 文件中的第6行，全部變成2。操作：
#第一步： mv操作,重命名.fam文件為tmp文件，這樣可以在上tmp上面修改
mv xxxx.fam tmp
#第二步：awk操作，把第六行全部變成2，然后這個東西再重新寫入.fam文件
awk '{print $1,$2,$3,$4,$5,$6=2}' tmp > xxxx.fam
  • 在上述情況下，或者如果fam中沒有表型文件，需要另外再創(chuàng)建一個表型文件：
    表型文件一般表型文件為txt格式，表型文件有三列，分別為：
    FID（Family ID）、IID（Individual ID）、Pheno（Phenotype）
    在命令語句中加入【--pheno xxx.txt】
  • 用bar-level設(shè)定你想計(jì)算的閾值，或者PRSice會默認(rèn)計(jì)算upper和lower間一定間隔的p值下的所有PRS，命令語句中加入【--all-score】命令可以得到這所有的PRS結(jié)果

四、運(yùn)行結(jié)果解釋

運(yùn)行完畢后，PRSice文件夾中應(yīng)該有如下的文件，兩個圖和幾個文本文件


  • PRSice_BARPLOT_.png 顯示了不同閾值得到的多基因評分的對應(yīng) R2 分布。柱狀圖最高的點(diǎn)表示模型最優(yōu)，如該圖顯示的是P值閾值為0.4463時(shí)，模型最優(yōu)，R方=0.0520082（P=4.7*10-18）這些信息也可以在.summary文件中查到 。


    BARPLOT

  • PRSice_HIGH-RES_PLOT_ .png條形圖顯示了許多不同的 p值閾值，以及對應(yīng) R2 的 p 值（黑點(diǎn)），綠線為趨勢線。y 軸最高點(diǎn)代表最優(yōu)解best-fit PRS
    這里最佳的P值閾值為綠色最高點(diǎn)處，此時(shí)P值的閾值為0.4463


    high-res

    這兩個圖均表明，許多影響 base 樣本性狀的 SNP 可用于預(yù)測 target 樣品中的性狀。 注意，這兩個性狀可以相同也可以不同。 如果使用相同的性狀，則預(yù)測值與該性狀的遺傳性（以及 base 樣本的樣本量）有關(guān)。 如果分析不同的性狀，則預(yù)測值與兩個性狀之間的遺傳重疊有關(guān)。無論哪種方式，多基因風(fēng)險(xiǎn)評分分析通常都表明，寬松 p 值閾值的模型通常比更嚴(yán)格閾值的模型有更好的預(yù)測能力，這表明許多統(tǒng)計(jì)學(xué)上無關(guān)緊要的 SNP 在多基因性狀上具有預(yù)測價(jià)值。

  • PRSice.prsice文件包含：在給定不同閾值的P值以后，符合要求的SNP數(shù)量（Num_SNP）,SNP的解釋度（R2），回歸系數(shù)
  • PRSice.best文件包含F(xiàn)ID,IID,是否回歸，PRS值。這個文件計(jì)算的是每個個體最優(yōu)的PRS值。
  • PRSice.summary文件，包含表型，P的閾值，PRS的解釋方差，所有變量的解釋方差，協(xié)變量的解釋方差，回歸系數(shù)，P值，該閾值下的SNP數(shù)量。 這個文件給出的是該表型下最優(yōu)的模型。
  • Clumping ：識別并選擇每個LD block中最顯著的 SNP（即 p 值最低）以進(jìn)行進(jìn)一步分析。這樣可以減少 SNP 之間的相關(guān)性，同時(shí)保留具有最強(qiáng)統(tǒng)計(jì)證據(jù)的 SNP。

附1：OR值與log OR值的轉(zhuǎn)換

下載PGC數(shù)據(jù)庫關(guān)于精神分裂癥的數(shù)據(jù)https://www.med.unc.edu/pgc/data-index/
下載完成為.gz文件，解壓語句：

tar -xzvf file.tar.gz      #解壓xxxtar.gz文件
gunzip FileName.gz    #解壓xxx.gz文件

我下載的為第二種文件，解壓完成后得到一個文稿，可以用txt打開
打開之后檢查一下文稿里面是OR值 或者 log OR值（計(jì)算PRS需要用OR），可以用R進(jìn)行轉(zhuǎn)換：

dat <- read table("xxxx.txt",header=T)
dat$OR <- exp(dat$LogOR)
write.table(dat,"xxxtransformed.txt",quote=F, row.names=F)
q()
  • 文本文件可以打開txt查看

附2:萬能查看文件語句

用txt可以打開任何一個文件查看，什么后綴的都可以但是?。。。。?！排列不整齊?。。?！進(jìn)行后續(xù)篩選非常不方便！?。。。?！
這個時(shí)候sort -o可以幫你

sort -o happy.txt sad.assoc

上面的命令就可以把關(guān)聯(lián)分析得到的sad.assoc文件 另存為 happy.txt，然后就可以用EXCEL打開了，盡情查看吧。

參考資料

http://www.360doc.com/content/19/1224/13/68068867_881784568.shtml
https://choishingwan.github.io/PRSice/step_by_step/
http://m.itdecent.cn/p/636048889b2a
http://m.itdecent.cn/p/656573127d11
https://www.cnblogs.com/chenwenyan/p/10686136.html
http://m.itdecent.cn/c/7df7c15887bd

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容