老规矩，还是资料整理者，边学习边做笔记，可能会改动原文章，按自己的理解将多篇优秀的推文自由整合，后面会附上所有引用推文的出处，感谢各位作者！
先说说PRS是用来干什么的？
The PRS is used to assess the cumulative genetic risk for a certain disorder (Purcell SM, et al. Common polygenic variation contributes to risk of schizophrenia and bipolar disorder. Nature 2009; 460: 748–52).
计算PRS的工具：PRSice software (http://prsice.info) (Euesden J. PRSice: polygenic risk score software.Bioinformatics 2015; 31: 1466–8.)
***注意：base和target样本不能是同一组人，会引起beta值的膨胀，如果要用同一组人，则必须用不同表型，而且这两个表型必须不相关，这样的话，其实做PRS也没什么意义了。所以，base和target要用不同的样本。
PRS常被计算为个体携带的风险等位基因数量的加权总和，其中风险等位基因及其权重由基因座定义及其测量的效应取决于全基因组关联研究检测。
image
n: 纳入的SNP的数量；
Xi：SNP I的基因型（0，1，2）
Wi：GWAS发现的SNP的权重，其中：
连续型变量的表型用Effect size表示，beta值
分类变量的表型用OR值表示
————————————————————————————
原来讲完啦，开始准备工具：
先得下载PRSice软件，官方网站https://choishingwan.github.io/PRSice/，现在是版本2.1.4，根据运行环境选择要下载的版本，我选的是linux64-bit 2.1.4，下载之后长这样，先拿他的示例试着跑一下。
image
————————————————————————————
工具已在手，开始准备“材料”：（这个过程是比较繁琐的，但也是很关键的一步，一定要核对好）
input文件包括4个：（根据自己的数据对应整理好）
PRSice.R file:上面的PRSice.R，画bar plot，high-resolution plot and quantile plot
PRSice executable file:上面的PRSice_linux，核心运算
Base data set:上面的TOY_BASE_GWAS.assoc
Target data set:上面的TOY_TARGET_DATA(3个)
得到两种“材料“后，就开始整理这些“材料”，看看都是否符合标准，就是我们说的质控啦，这步很关键哦。
质控：
Target Data质控：
– Sample size：我们建议用户仅对至少100个人的目标数据执行PRS分析。
– Genome build：基因组版本要和参考数据一致（这点很重要，两组数据需要用snp来匹配，必须确定相同版本的名称）
– Standard GWAS QC：目标数据的质量必须至少达到GWAS研究中实施的标准，例如 从Hardy-Weinberg平衡中去除具有低基因分型率，低等位基因频率的SNP，去除具有低基因分型率的个体。
plink \
--bfile EUR \ #文件前辍EUR，指定输入
--maf 0.05 \ # 删除所有等位基因频率小于0.05的SNP。基因分型错误通常对MAF低的SNP具有较大影响。样本量较大的研究可以应用较低的MAF阈值
--hwe 1e-6 \ #从Hardy-Weinberg平衡Fisher的精确检验或卡方检验中删除了低P值的SNP。HWE测试中具有显著P值的SNP更可能受到基因分型错误或自然选择的影响。应对对照样品进行过滤，以避免过滤因果关系的SNP（在某些情况下应选择）
--geno 0.01 \ #排除了在大部分受试者中缺失的SNP。通常执行两阶段过滤过程（请参阅Marees等）。
--mind 0.01 \ #排除基因型缺失率很高的人，因为这可能表明DNA样品或加工中存在问题。
--write-snplist \
--make-just-fam \
--out EUR.QC
个体中很高或很低的杂合率可能是由于DNA污染或近交水平高。因此，具有极高杂合性的样品通常在下游分析之前被去除。首先，我们执行修剪以删除高度相关的SNP：
plink \
--bfile EUR \
--keep EUR.QC.fam \
--extract EUR.QC.snplist \
--indep-pairwise 200 50 0.25 \ #EUR.QC.prune.in中的所有SNP都成对r2(平方)<0.25.
--out EUR.QC
然后可以使用plink计算杂合率,
plink \
--bfile EUR \
--extract EUR.QC.prune.in \
--keep EUR.QC.fam \
--het \
--out EUR.QC
这将生成EUR.QC.het文件，其中包含用于评估杂合性的F系数估计。我们将删除F系数与均值相比超过3个标准差（SD）单位的个人，可以使用以下R命令执行此操作。
library(data.table)
# Read in file
dat <- fread("EUR.QC.het")
# Get samples with F coefficient within 3 SD of the population mean
valid <- dat[F<=mean(F)+3*sd(F) & F>=mean(F)-3*sd(F)]
# print FID and IID for valid samples
fwrite(valid[,c("FID","IID")], "EUR.valid.sample", sep="\t")
– Ambiguous SNPs
在基础数据质量控制期间已将其删除
– Sex chromosomes
有时可能会出现样品贴错标签的情况，这可能导致无效的结果。样品标签错误的一个很好的指示是生物学性别和报告的性别不匹配。如果生物学性别与所报告的性别不符，则该样品可能贴错标签。
在进行性别检查之前，应进行修剪（请参见此处）。然后可以使用plink轻松进行性别检查：
plink \
--bfile EUR \
--extract EUR.QC.prune.in \
--keep EUR.valid.sample \
--check-sex \
--out EUR.QC
这将生成一个名为EUR.QC.sexcheck的文件，其中包含每个人的F统计信息。如果F统计量> 0.8，通常将个体称为生物学上的男性，如果F <0.2，则将个体称为生物学上的女性。
library(data.table)
#读入文件
valid <- 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")
— Mismatching genotypes
此外，当数据集之间的等位基因编码存在明确的不匹配时，例如基础中的A / C和目标数据中的G / T，则可以通过“链翻转”目标数据中的等位基因来解决互补的等位基因(大多数PRS软件将自动执行翻转，因此通常不需要此步骤)。这可以通过以下步骤实现：
a.将bim文件，GIANT摘要统计信息和QC SNP列表加载到R中:
library(data.table)
# Read in bim file
bim <- fread("EUR.bim")
setnames(bim, colnames(bim), c("CHR", "SNP", "CM", "BP", "B.A1", "B.A2"))
# Read in GIANT data (require data.table v1.12.0+)
height <- fread("Height.QC.gz")
# Change all alleles to upper case for easy comparison
height[,c("A1","A2"):=list(toupper(A1), toupper(A2))]
bim[,c("B.A1","B.A2"):=list(toupper(B.A1), toupper(B.A2))]
# Read in QCed SNPs
qc <- fread("EUR.QC.snplist", header=F)
b.识别需要链翻转的SNP
# Merge GIANT with target
info <- merge(bim, height, by=c("SNP", "CHR", "BP"))
info <- info[SNP %in% qc$V1]
# Function for calculating the complementary allele
complement <- function(x){
switch (x,
"A" = "T",
"C" = "G",
"T" = "A",
"G" = "C",
return(NA)
)
}
# Identify SNPs that are complementary between base and target
com.snps <- info[sapply(B.A1, complement) == A1 &
sapply(B.A2, complement) == A2, SNP]
# Now update the bim file
bim[SNP %in% com.snps, c("B.A1", "B.A2") :=
list(sapply(B.A1, complement),
sapply(B.A2, complement))]
# And update the info structure
info[SNP %in% com.snps, c("B.A1", "B.A2") :=
list(sapply(B.A1, complement),
sapply(B.A2, complement))]
c.识别需要在目标中重新编码的SNP（以确保目标数据中的编码等位基因是基本摘要统计中的有效等位基因）
# identify SNPs that need recoding & complement
com.recode <- info[sapply(B.A1, complement) == A2 &
sapply(B.A2, complement) == A1, SNP]
# Now update the bim file
bim[SNP %in% com.recode, c("B.A1", "B.A2") :=
list(sapply(B.A2, complement),
sapply(B.A1, complement))]
# And update the info structure
info[SNP %in% com.recode, c("B.A1", "B.A2") :=
list(sapply(B.A2, complement),
sapply(B.A1, complement))]
# Write the updated bim file
fwrite(bim, "EUR.QC.adj.bim", col.names=F, sep="\t")
d.识别在碱基和靶标中具有等位基因不同的SNP（通常是由于基因组构建或Indel的差异）
matched <- info[(A1 == B.A1 & A2 == B.A2) |
(A1 == B.A2 & A2 == B.A1)]
mismatch <- bim[!SNP%in%matched$SNP, SNP]
write.table(mismatch, "EUR.mismatch", quote=F, row.names=F, col.names=F)
e.将EUR.bim替换为EUR.QC.adj.bim：
# Make a back up
mv EUR.bim EUR.bim.bk
ln -s EUR.QC.adj.bim EUR.bim
--重复SNP
确保删除目标数据中的所有重复SNP（这些目标数据是模拟的，因此不包含重复的SNP）
--Sample overlap
由于目标数据是模拟的，因此此处的基础数据和目标数据之间没有重叠的样本（有关避免样本重叠的重要性的讨论，请参见本文的相关部分）。
--Relatedness
目标数据中关系密切的个人可能会导致结果过拟合，从而限制了结果的通用性。在计算相关性之前，应进行修剪（请参见此处）。样品中具有一级或二级亲属的个人（π>0.125）可以使用以下命令删除：
plink \
--bfile EUR \
--extract EUR.QC.prune.in \
--keep EUR.QC.valid \
--rel-cutoff 0.125 \
--out EUR.QC
贪心算法用于以最优化保留样本大小的方式删除紧密相关的个体。但是，该算法取决于所使用的随机种子，这可能会产生不同的结果。因此，要重现相同的结果，您将需要指定相同的随机种子。
PLINK的去除相关个体的算法不能解释所研究的表型。为了最大程度地减少疾病的清除，可以使用以下算法代替：GreedyRelated。
Generate final QC'ed target data file
执行完全面分析后，您可以使用以下命令生成质量控制数据集：
plink \
--bfile EUR \
--make-bed \
--keep EUR.QC.rel.id \
--out EUR.QC \
--extract EUR.QC.snplist
以上质控这段参考（这个写的很棒）：https://mp.weixin.qq.com/s/xDUrZvJwvBXJk7TDj_Ke0Q
需要注意的问题：
1.若target中fam文件最后一列表型缺失（-9或NA），则需要另外加表型文件，加命令： --pheno-file
2.bp、chr、snp：软件会自动审核base和target文件中这些信息是否匹配，不匹配则去除，一定要注意文件名，大小写也需要注意，column名称必须完全一致，顺序不对应没关系，会自动匹配。
3.亲测了下，不一定要 .assoc.linear文件，我用的txt文件，里面包含了所有它规定的内容，也可以跑通。
———————————————————————————————
工具、材料都在手啦~开始启动啦:就那么一条命令就能搞定。
直接在linux下输入以下命令：
Rscript PRSice.R --dir . --prsice ./PRSice_linux --base TOY_BASE_GWAS.assoc --target TOY_TARGET_DATA --thread 1 --stat OR --binary-target T
这里涉及一个是连续变量还是二进制变量的问题，改换不同的BEAT（连续）或OR（二进制），后面的F（连续）或T（二进制）
image
结果报错了，原因不明，后来输入命令
chmod 775 PRSice_linux
再返回输入刚才那条命令，就好了。结果文件就生成了。也是神奇了，上面那条命令啥意思啊没弄明白，如果有高人看到了帮忙指点下。
image
附：如果是下载V1版本的，需要在Linux下启动R，把需要的R包都安装好，再退出R，再linux下输入命令即可。需要的文件一定要准备好。
自己的数据上一遍：
input：PRSice.R PRSice_linux TAU_age_gwas.assoc.linear ADNI1_Genotypes_Filt_CEU_final.bed ADNI1_Genotypes_Filt_CEU_final.bim ADNI1_Genotypes_Filt_CEU_final.fam
命令：
chmod 775 PRSice_linux
Rscript PRSice.R --dir . --prsice ./PRSice_linux --base TAU_age_gwas.assoc.linear --target ADNI1_Genotypes_Filt_CEU_final --thread 1 --stat BETA --binary-target F
输出的文件
image
展示一下各结果文件
image
image
image
image
image
PRSice will automatically calculate the PRS for different p-value thresholds and perform a regression to test the level of association of the PRS with the target phenotype. This allow the identification of PRS that "best" predicts the phenotype and can be used for downstream analysis.
It is vital that the human genome build is the same for the GTF file, bed files, target file and the base file. Otherwise the coordinates of the SNPs can be wrong and PRSice will not be able to correctly assign the gene membership, leading to invalid results.
Nagelkerke’s pseudo R2 was calculated to measure the proportion of variance explained by the PRS.
参考推文：
1.https://www.jianshu.com/p/656573127d11
2.https://www.cnblogs.com/chenwenyan/p/10686136.html
3.英文github官网https://choishingwan.github.io/PRSice/step_by_step/#input-data
网友评论