翻譯自：https://choishingwan.github.io/PRS-Tutorial/
需要安裝：
1.R[1] (version 3.2.3+)2.PLINK 1.9[2]
多基因風險評分（Polygenic Risk Score）分析過程概覽。PRS 分析需要兩個輸入數據集：i）base data（GWAS）：全基因組範圍內遺傳變異的基因型-表型關聯的摘要統計信息（例如 beta，P值） ；ii）target data：目標樣本中個體的基因型和表型。基於 base data 得到的 SNP 效應值計算 target data 中樣本的 PRS。
教程分為四個部分：
1.Base Data 的質控2.Target Data 的質控3.計算和分析 PRS4.可視化 PRS 結果
示例數據下載地址（需翻牆）：
# base datahttps://drive.google.com/file/d/1RWjk49QNZj9zvJHc9X_wyZ51fdy6xQjv/view?usp=sharing# target datahttps://drive.google.com/file/d/1uhJR_3sn7RA8U5iYQbcmTp6vFdQiF4F2/view?usp=sharing我也把兩個數據集傳到了堅果雲上，方便大家下載：https://www.jianguoyun.com/p/DQ53W3sQ757XCBik57AD
Base Data 的質控PRS 分析的第一步就是得到 base data 的 GWAS 結果數據。查看 Height.gwas.txt.gz示例文件：
gunzip -c Height.gwas.txt.gz | headCHR BP SNP A1 A2 N SE P OR INFO MAF1 756604 rs3131962 A G 388028 0.00301666 0.4831710.997886915712657 0.890557941364774 0.3693895927649211 768448 rs12562034 A G 388028 0.00329472 0.8348081.00068731609353 0.895893511351165 0.3368457540962891 779322 rs4040617 G A 388028 0.00303344 0.428970.997603556067569 0.897508290615237 0.3773680109408141 801536 rs79373928 G T 388028 0.00841324 0.8089991.00203569922793 0.908962856432993 0.4832122453740951 808631 rs11240779 G A 388028 0.00242821 0.5902651.00130832511154 0.893212523690488 0.4504095589995871 809876 rs57181708 G A 388028 0.00336785 0.714751.00123165786833 0.923557624081969 0.4997439326567591 835499 rs4422948 G A 388028 0.0023758 0.7108840.999119752645202 0.906437735120596 0.4810160058161681 838555 rs4970383 A C 388028 0.00235773 0.1509930.996619945289758 0.907716506801574 0.3271640296727541 840753 rs4970382 C T 388028 0.00207377 0.1999670.99734567895614 0.914602590137255 0.498936220426316包含以下幾列：
•CHR：SNP 所在的染色體•BP：SNP 的座標•SNP：SNP ID，通常採用 rs-ID 的形式•A1：SNP 的 effect allele•A2：SNP 的 non-effect allele•N：用於計算效應值的樣本數•SE：效應量的標準差•P：SNP 基因型與表型關聯的 p 值•OR：SNP 的效應量（case-control）。若為連續的表型，那這一列通常為 beta•INFO：imputation information score•MAF：minor allele frequency
檢查遺傳度我們建議 base data 計算得到的 h2 > 0.05（遺傳度可用 LDSC 計算）。
Effect allele一些 GWAS 結果文件沒有明確哪個等位基因是效應等位基因，哪個是非效應等位基因。如果在計算 PRS 時進行了錯誤的假設，那麼 PRS 對 target data 的效應估算將是錯誤的。
檢查 GWAS 結果文件的完整性另一個常見的問題是，下載的 base data 文件可能在下載過程中損壞，這可能導致 PRS 軟體崩潰或在產生錯誤的結果。我們可用 md5sum檢查文件的完整性：
md5sum Height.gwas.txt.gz參考基因組我們還需要檢查 base data 和 target data 是否使用了相同的參考基因組。若使用了不同的參考基因組，則需要用 LiftOver[3] 進行坐標轉換，使他們保持一致。
GWAS 結果 QCBase data 和 target data 都應執行嚴格的質量控制步驟。低 MAF 和 INFO 的 SNPs 在執行下遊分析之前通常會被去除。我們建議移除 MAF < 1% 和 INFO < 0.8 的 SNP。這一步可通過以下代碼實現:
gunzip -c Height.gwas.txt.gz |\awk 'NR==1 || ($11 > 0.01) && ($10 > 0.8) {print}' |\gzip > Height.gz等位基因不匹配在 base data 和 target data 中不匹配的 SNP 可通過「鏈翻轉」進行匹配，例如某個 SNP 在 base data 中為 A/C，target data 中為 G/T，亦或者是一些不可解析的 SNP，例如在 base data 中為 C/G，target data 中為 C/T。大多數 PRS 軟體會對可解析的 SNPs 進行自動鏈翻轉，並刪除不可解析的 SNPs。
因為我們需要 target data 來知道哪些 SNPs 具有不匹配的等位基因，所以我們將在 target data 中執行這種鏈翻轉。
重複的 SNPs大多數 PRS 軟體不允許 base data 中存在重複的 SNPs，我們可以用下面的命令刪除它們：
# 輸出重複的 SNP IDgunzip -c Height.gz |\awk '{ print $3}' |\sort |\uniq -d > duplicated.snp# 刪除重複的 SNPgunzip -c Height.gz |\grep -vf duplicated.snp |\gzip - > Height.nodup.gz去除一些可能產生歧義的 SNPs可以使用以下方法保留無歧義的 SNPs:
gunzip -c Height.nodup.gz |\awk '!( ($4=="A" && $5=="T") || \ ($4=="T" && $5=="A") || \ ($4=="G" && $5=="C") || \ ($4=="C" && $5=="G")) {print}' |\ gzip > Height.QC.gz檢查樣本性別有時樣品標籤會出現錯誤，最典型的問題就是通過 plink 計算出的性別與記錄的性別有差。這部分請參閱 target data 的質控部分。
樣本重複由於示例的 target data 是模擬產生的，所以這裡的 base data 和 target data 之間沒有重疊的樣本。
親緣關係base data 和 target data 中有親緣關係的樣本可能導致過擬合的結果。此部分可參考 target data 的質控部分。
最終，Height.QC.gz 文件可用於下遊分析。
Target Data 的質控樣本量我們建議 target data 的樣本量大於 100。示例數據的樣本量是 503 人。
基因型數據質控plink \ --bfile EUR \ --maf 0.01 \ --hwe 1e-6 \ --geno 0.01 \ --mind 0.01 \ --write-snplist \ --make-just-fam \ --out EUR.QC主要參數的含義如下：
•--bfile：輸入基因型文件•--maf：刪除 MAF 小於 0.01 的 SNPs•--hwe：刪除 HWE p 值低於 1e-6 的 SNPs•--geno：排除大部分樣本中缺失的 SNPs•--mind：刪除基因型缺失率高的樣本
執行 prunning 來刪除高度相關的 SNPs:
plink \ --bfile EUR \ --keep EUR.QC.fam \ --extract EUR.QC.snplist \ --indep-pairwise 200 50 0.25 \ --out EUR.QC主要參數的含義如下：
•--keep：只保留 EUR.QC.fam 中存在的樣本•--extract：只保留 EUR.QC.snplist 中的 SNPs•--indep-pairwise：執行 prunning 的窗口大小為 200kb，每次步長為 50 個 SNP，並過濾掉 LD r2 大於 0.25 的 SNPs
這行命令將生成兩個文件1) EUR.QC.prune.in 和2) EUR.QC.prune.out。
下一步我們需要排除樣本雜合率均值偏離 ±3 SD 的個體。用 plink 計算雜合率：
plink \ --bfile EUR \ --extract EUR.QC.prune.in \ --keep EUR.QC.fam \ --het \ --out EUR.QC這行代碼會生成 EUR.QC.het 文件。用 R 代碼進行過濾:
library(data.table)# Read in filedat <- fread("EUR.QC.het")# Get samples with F coefficient within 3 SD of the population meanvalid <- dat[F<=mean(F)+3*sd(F) & F>=mean(F)-3*sd(F)] # print FID and IID for valid samplesfwrite(valid[,c("FID","IID")], "EUR.valid.sample", sep="\t") q() # exit R等位基因不匹配的 SNPs在 base data 和 target data 中等位基因不匹配的 SNPs，可通過將等位基因翻轉到它們的互補等位基因來解決。R 代碼如下。
1.首先讀入數據：
# magrittr allow us to do piping, which help to reduce the # amount of intermediate data typeslibrary(data.table)library(magrittr)# Read in bim file bim <- fread("EUR.bim") %>% # Note: . represents the output from previous step # The syntax here means, setnames of the data read from # the bim file, and replace the original column names by # the new names setnames(., colnames(.), c("CHR", "SNP", "CM", "BP", "B.A1", "B.A2")) %>% # And immediately change the alleles to upper cases .[,c("B.A1","B.A2"):=list(toupper(B.A1), toupper(B.A2))]# Read in summary statistic data (require data.table v1.12.0+)height <- fread("Height.QC.gz") %>% # And immediately change the alleles to upper cases .[,c("A1","A2"):=list(toupper(A1), toupper(A2))]# Read in QCed SNPsqc <- fread("EUR.QC.snplist", header=F)2.識別需要翻轉 DNA 鏈的 SNPs
# Merge summary statistic with targetinfo <- merge(bim, height, by=c("SNP", "CHR", "BP")) %>% # And filter out QCed SNPs .[SNP %in% qc[,V1]]# Function for calculating the complementary allelecomplement <- function(x){ switch (x, "A" = "T", "C" = "G", "T" = "A", "G" = "C", return(NA) )} # Get SNPs that have the same alleles across base and targetinfo.match <- info[A1 == B.A1 & A2 == B.A2, SNP]# Identify SNPs that are complementary between base and targetcom.snps <- info[sapply(B.A1, complement) == A1 & sapply(B.A2, complement) == A2, SNP]# Now update the bim filebim[SNP %in% com.snps, c("B.A1", "B.A2") := list(sapply(B.A1, complement), sapply(B.A2, complement))]3.識別需要在 target data 中重新編碼的 SNPs
# identify SNPs that need recodingrecode.snps <- info[B.A1==A2 & B.A2==A1, SNP]# Update the bim filebim[SNP %in% recode.snps, c("B.A1", "B.A2") := list(B.A2, B.A1)]# identify SNPs that need recoding & complementcom.recode <- info[sapply(B.A1, complement) == A2 & sapply(B.A2, complement) == A1, SNP]# Now update the bim filebim[SNP %in% com.recode, c("B.A1", "B.A2") := list(sapply(B.A2, complement), sapply(B.A1, complement))]# Write the updated bim filefwrite(bim, "EUR.QC.adj.bim", col.names=F, sep="\t")4.確定在 base 和 target 上有不同等位基因的 SNPs (通常是由於使用不同參考基因組或 Indel造成的)
mismatch <- bim[!(SNP %in% info.match | SNP %in% com.snps | SNP %in% recode.snps | SNP %in% com.recode), SNP]write.table(mismatch, "EUR.mismatch", quote=F, row.names=F, col.names=F)q() # exit R5.將 EUR.bim 替換為 EUR.QC.adj.bim：
# Make a back upmv EUR.bim EUR.bim.bkln -s EUR.QC.adj.bim EUR.bim注意：大多數 PRS 軟體將自動執行這種轉換，所以一般不需要這一步。
重複的 SNPs確保刪除 target data 中重複的 SNPs (示例 target data 是模擬產生的，因此不包括重複的 SNPs)。
檢查性別我們可根據 X 染色體雜合/純合率檢查數據集中記錄的性別與真實性別之間的差異。男性 X 染色體純合度估計值 > 0.8，女性 <0.2 。
在執行性別檢查之前，應先進行 pruning：
plink \ --bfile EUR \ --extract EUR.QC.prune.in \ --keep EUR.valid.sample \ --check-sex \ --out EUR.QC這將生成 EUR.QC.sexcheck 的文件，根據這個文件進行過濾，R 代碼如下：
library(data.table)# Read in filevalid <- fread("EUR.valid.sample")dat <- fread("EUR.QC.sexcheck")[FID%in%valid$FID]fwrite(dat[STATUS=="OK",c("FID","IID")], "EUR.QC.valid", sep="\t") q() # exit R樣本重疊由於 target data 是模擬的，所以這裡的 base data 和 target data 之間沒有重疊的樣本。
親緣關係Target data 中密切相關的個體可能導致過擬合，需要進行過濾。
plink \ --bfile EUR \ --extract EUR.QC.prune.in \ --keep EUR.QC.valid \ --rel-cutoff 0.125 \ --out EUR.QC生成最終 QC 後的 target data 文件plink \ --bfile EUR \ --make-bed \ --keep EUR.QC.rel.id \ --out EUR.QC \ --extract EUR.QC.snplist \ --exclude EUR.mismatch計算和分析 PRS用 plink 計算 PRS根據前面的質控步驟，我們得到了以下 六個文件：
•Height.QC.gz•EUR.QC.bed•EUR.QC.bim•EUR.QC.fam•EUR.height•EUR.covariate
更新效應量這裡為了簡化計算，我們使用 OR 的自然對數，這樣可以用求和來計算 PRS。
library(data.table)dat <- fread("Height.QC.gz")fwrite(dat[,BETA:=log(OR)], "Height.QC.Transformed", sep="\t")q() # exit R由於數值的四捨五入，使用 awk 來進行轉換可能會導致不太精確的結果。因此，我們建議在 R 中執行轉換，或者用 PRS 軟體直接執行轉換。
Clumping儘管原則上所有常見的 SNP 都可以用於 PRS 分析中，但習慣上先將 GWAS 結果 clump，然後再計算風險評分。所謂 clumping 就是識別並選擇每個 LD block 中最顯著的 SNP（即 p 值最低）以進行進一步分析。這樣可以減少 SNP 之間的相關性，同時保留具有最強統計證據的 SNP。
plink \ --bfile EUR.QC \ --clump-p1 1 \ --clump-r2 0.1 \ --clump-kb 250 \ --clump Height.QC.Transformed \ --clump-snp-field SNP \ --clump-field P \ --out EUR主要參數的含義如下：
•--clump-p1：index SNP 的 P 值，取 1 即是納入所有的 SNPs•--clump-r2：r2 > 0.1 的 SNP 將被刪除•--clump-kb：範圍取 250kb•--clump：base data 的 summary statistic 結果文件•--clump-snp-field：指定包含 SNP ID 的列名•--clump-field：指定 P 值的列名
上面的命令將得到 EUR.clumped。我們可以通過執行以下命令提取 index SNP ID:
awk 'NR!=1{print $3}' EUR.clumped > EUR.valid.snp如果你的 target data 很小(比如 n < 500) ，那麼你可以用千人基因組計劃的數據來計算 LD。確保使用最接近反映 base 樣本的人群數據。
計算 PRS在 plink 中我們可使用 --score 和 --q-score-range 計算 PRS。
這裡我們需要三個文件：
1.base data 文件：Height.QC.Transformed2.一個包含 SNP id 及其對應的 p 值的文件：
awk '{print $3,$8}' Height.QC.Transformed > SNP.pvalue1.一個包含不同 p 值閾值的文件：
echo "0.001 0 0.001" > range_listecho "0.05 0 0.05" >> range_listecho "0.1 0 0.1" >> range_listecho "0.2 0 0.2" >> range_listecho "0.3 0 0.3" >> range_listecho "0.4 0 0.4" >> range_listecho "0.5 0 0.5" >> range_list三列分別表示：
例如對於 0.05 的閾值，我們包括 p 值從 0 到 0.05 的所有 SNPs，包括 p 值等於 0.05 的 SNPs。
然後我們就可以用 plink 計算 PRS 了：
plink \ --bfile EUR.QC \ --score Height.QC.Transformed 3 4 12 header \ --q-score-range range_list SNP.pvalue \ --extract EUR.valid.snp \ --out EUR參數含義如下:
•--score Height.QC.Transformed 3 4 12 header：輸入 Height.QC.Transformed 文件，第三列為 SNP ID，第四列為 effective allele，第十二列為效應量，且該文件包含 header。•--q-score-range range_list SNP.pvalue：定義 PRS 閾值文件，和 SNP P 值文件。
上面的命令和 range_list 將生成7個文件:
1.EUR.0.5.profile2.EUR.0.4.profile3.EUR.0.3.profile4.EUR.0.2.profile5.EUR.0.1.profile6.EUR.0.05.profile7.EUR.0.001.profile
考慮人群分層因素人口結構是 GWAS 中混雜因素的主要來源，我們通常將主成分(PC)作為協變量進行校正。
用 plink 計算主成分：
# First, we need to perform prunningplink \ --bfile EUR.QC \ --indep-pairwise 200 50 0.25 \ --out EUR# Then we calculate the first 6 PCsplink \ --bfile EUR.QC \ --extract EUR.prune.in \ --pca 6 \ --out EUR我們可以用 LDSC 來挑選合適的主成分數量。使用不同數量的主成分進行校正，LDSC 截距最接近 1 的即是最佳的主成分數量。如果 base 樣本和 target 樣本是從世界各地不同的人群中收集的，PRS 分析的結果可能存在偏差。
尋找最佳的 PRS 閾值為了接近「最佳擬合」的 PRS，我們可在一定 P 值範圍內計算 PRS，然後選擇能夠解釋最高表型差異的 PRS。這可以通過 R來實現：
library(data.table)library(magrittr)p.threshold <- c(0.001,0.05,0.1,0.2,0.3,0.4,0.5)phenotype <- fread("EUR.height")pcs <- fread("EUR.eigenvec", header=F) %>% setnames(., colnames(.), c("FID", "IID", paste0("PC",1:6)) )covariate <- fread("EUR.cov")pheno <- merge(phenotype, covariate) %>% merge(., pcs)null.r2 <- summary(lm(Height~., data=pheno[,-c("FID", "IID")]))$r.squaredprs.result <- NULLfor(i in p.threshold){ pheno.prs <- paste0("EUR.", i, ".profile") %>% fread(.) %>% .[,c("FID", "IID", "SCORE")] %>% merge(., pheno, by=c("FID", "IID")) model <- lm(Height~., data=pheno.prs[,-c("FID","IID")]) %>% summary model.r2 <- model$r.squared prs.r2 <- model.r2-null.r2 prs.coef <- model$coeff["SCORE",] prs.result %<>% rbind(., data.frame(Threshold=i, R2=prs.r2, P=as.numeric(prs.coef[4]), BETA=as.numeric(prs.coef[1]), SE=as.numeric(prs.coef[2])))}print(prs.result[which.max(prs.result$R2),])q() # exit R用 PRSice-2 計算 PRS在運行 PRSice-2 之前，我們還需要一個協變量文件，整合協變量信息和主成分信息。R 代碼如下:
library(data.table)covariate <- fread("EUR.covariate")pcs <- fread("EUR.eigenvec", header=F)colnames(pcs) <- c("FID","IID", paste0("PC",1:6))cov <- merge(covariate, pcs)fwrite(cov,"EUR.cov", sep="\t")q()生成 EUR.cov 文件。
然後可以運行 PRSice-2，計算 PRS：
Rscript PRSice.R \ --prsice PRSice_linux \ --base Height.QC.gz \ --target EUR.QC \ --binary-target F \ --pheno EUR.height \ --cov EUR.cov \ --base-maf MAF:0.01 \ --base-info INFO:0.8 \ --stat OR \ --or \ --out EUR這行代碼會自動尋找最佳的參數並進行可視化。
可視化 PRS 結果我們可以用 ggplot2 進行可視化：
# ggplot2 is a handy package for plottinglibrary(ggplot2)# generate a pretty format for p-value outputprs.result$print.p <- round(prs.result$P, digits = 3)prs.result$print.p[!is.na(prs.result$print.p) & prs.result$print.p == 0] <- format(prs.result$P[!is.na(prs.result$print.p) & prs.result$print.p == 0], digits = 2)prs.result$print.p <- sub("e", "*x*10^", prs.result$print.p)# Initialize ggplot, requiring the threshold as the x-axis (use factor so that it is uniformly distributed)ggplot(data = prs.result, aes(x = factor(Threshold), y = R2)) + # Specify that we want to print p-value on top of the bars geom_text( aes(label = paste(print.p)), vjust = -1.5, hjust = 0, angle = 45, cex = 4, parse = T ) + # Specify the range of the plot, *1.25 to provide enough space for the p-values scale_y_continuous(limits = c(0, max(prs.result$R2) * 1.25)) + # Specify the axis labels xlab(expression(italic(P) - value ~ threshold ~ (italic(P)[T]))) + ylab(expression(paste("PRS model fit: ", R ^ 2))) + # Draw a bar plot geom_bar(aes(fill = -log10(P)), stat = "identity") + # Specify the colors scale_fill_gradient2( low = "dodgerblue", high = "firebrick", mid = "dodgerblue", midpoint = 1e-4, name = bquote(atop(-log[10] ~ model, italic(P) - value),) ) + # Some beautification of the plot theme_classic() + theme( axis.title = element_text(face = "bold", size = 18), axis.text = element_text(size = 14), legend.title = element_text(face = "bold", size = 18), legend.text = element_text(size = 14), axis.text.x = element_text(angle = 45, hjust = 1) )# save the plotggsave("EUR.height.bar.png", height = 7, width = 7)q() # exit R我們還可以看看 PRS 與表型之間的關係：
library(ggplot2)# Read in the filesprs <- read.table("EUR.0.3.profile", header=T)height <- read.table("EUR.height", header=T)sex <- read.table("EUR.cov", header=T)# Rename the sexsex$Sex <- as.factor(sex$Sex)levels(sex$Sex) <- c("Male", "Female")# Merge the filesdat <- merge(merge(prs, height), sex)# Start plottingggplot(dat, aes(x=SCORE, y=Height, color=Sex))+ geom_point()+ theme_classic()+ labs(x="Polygenic Score", y="Height")q() # exit R引用連結[1] R: https://www.r-project.org/
[2] PLINK 1.9: https://www.cog-genomics.org/plink2
[3] LiftOver: https://genome.ucsc.edu/cgi-bin/hgLiftOver文末友情推薦
要想真正入門生物信息學建議務必購買全套書籍，一點一滴攻克計算機基礎知識，書單在：什麼，生信入門全套書籍僅需160 。如果大家沒有時間自行慢慢摸索著學習，可以考慮我們生信技能樹官方舉辦的學習班：
•數據挖掘學習班第5期（線上直播3周，馬拉松式陪伴，帶你入門），原價4800的數據挖掘全套課程， 疫情期間半價即可搶購。•生信爆款入門-第7期（線上直播4周，馬拉松式陪伴，帶你入門），原價9600的生信入門全套課程，疫情期間3.3折即可搶購。
如果你課題涉及到轉錄組，歡迎添加一對一客服：詳見：你還在花三五萬做一個單細胞轉錄組嗎？