GWAS - PRS多基因风险评分计算学习笔记

2020-05-31 本文已影响0人 SnowPye

一、安装PRSice（mac版）

经试验我觉得直接从git hub中下载对应的安装包是最快的：https://github.com/choishingwan/PRSice，下载之后解压，解压文件如图所示

PRSice_mac为软件的运行脚本
PRSice.R是执行脚本的封装
TOY开头的是数据集，分为BASE和TARGET两部分
*接下来可以先用这个数据集测试一下

跟plink的安装一样，打开terminal，cd到解压的文件夹，./PRSice_mac如果显示下面的画面表现为安装成功：

之后下载R包

Rscript PRSice.R --dir .

Rscript 表示用R脚本来调用该软件
dir参数为指定R包ggplot2安装的路径
因为结果展示会调用ggplot2画图进行可视化，如果你的R中已经安装了这个包，这个参数可以不要!!!!(我的R中有，所以后文中的此条命令我均不需要运行)
**因为直接运行该命令总容易报错，因此也可以通过R studio来安装这个包：点击面板中的package，选择ggplot2点击install，如果太慢的话选择清华的镜像就会快一点

二、文件准备

上面安装成功的图中可以看到PRSice的注释信息，需要准备base和target两个文件

Base Dataset：以空格分割的SNP关联分析文件，为.assoc结尾的。关联分析的教程见前面。
必须要有的几列是：SNP, A1, OR or BETA（对于分类型表型提供OR值）, P，其他几列不是必须要的，如下图：

BASE

如果表头没有以上述名称命名，必须用参数--chr, --A1, --A2, --stat, --snp, --bp, --se, --pvalue指明表头。
参数可以这样设置：--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产生的二进制文件.bed, .bim, 和 .fam file，要经过质控QC，过程翻看之前的教程。

三、PRSice运行

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

#数量性状用的是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参数指定可执行文件的路径；（安装过ggplot2的可以直接忽略掉这一条命令）
base参数指定base data（验证样本）的关联分析结果，包含每个遗传变异的效应值和 p 值。
target参数指定target data（测试样本）的分型结果，指 target 样本二进制 plink 格式的基因型数据前缀，也支持 BGEN 格式
thread参数指定线程数
stat指定base data关联分析结果中的效应值，有OR和BETA两种取值
binary-target参数指定target data中分型结果是否为二分类的表型，T=ture
·

如果下载下来的数据第六行并不是按我们的需求写上是表型的，为了后面的关联分析，要手工定义第六行是实验组，即需要把.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中没有表型文件，需要另外再创建一个表型文件：
表型文件一般表型文件为txt格式，表型文件有三列，分别为：
FID（Family ID）、IID（Individual ID）、Pheno（Phenotype）
在命令语句中加入【--pheno xxx.txt】
用bar-level设定你想计算的阈值，或者PRSice会默认计算upper和lower间一定间隔的p值下的所有PRS，命令语句中加入【--all-score】命令可以得到这所有的PRS结果

四、运行结果解释

运行完毕后，PRSice文件夹中应该有如下的文件，两个图和几个文本文件

PRSice_BARPLOT_.png 显示了不同阈值得到的多基因评分的对应 R2 分布。柱状图最高的点表示模型最优，如该图显示的是P值阈值为0.4463时，模型最优，R方=0.0520082（P=4.7*10-18）这些信息也可以在.summary文件中查到。

BARPLOT
PRSice_HIGH-RES_PLOT_ .png条形图显示了许多不同的 p值阈值，以及对应 R2 的 p 值（黑点），绿线为趋势线。y 轴最高点代表最优解best-fit PRS
这里最佳的P值阈值为绿色最高点处，此时P值的阈值为0.4463

high-res

这两个图均表明，许多影响 base 样本性状的 SNP 可用于预测 target 样品中的性状。注意，这两个性状可以相同也可以不同。如果使用相同的性状，则预测值与该性状的遗传性（以及 base 样本的样本量）有关。如果分析不同的性状，则预测值与两个性状之间的遗传重叠有关。无论哪种方式，多基因风险评分分析通常都表明，宽松 p 值阈值的模型通常比更严格阈值的模型有更好的预测能力，这表明许多统计学上无关紧要的 SNP 在多基因性状上具有预测价值。

PRSice.prsice文件包含：在给定不同阈值的P值以后，符合要求的SNP数量（Num_SNP）,SNP的解释度（R2），回归系数
PRSice.best文件包含FID,IID,是否回归，PRS值。这个文件计算的是每个个体最优的PRS值。
PRSice.summary文件，包含表型，P的阈值，PRS的解释方差，所有变量的解释方差，协变量的解释方差，回归系数，P值，该阈值下的SNP数量。这个文件给出的是该表型下最优的模型。
Clumping ：识别并选择每个LD block中最显著的 SNP（即 p 值最低）以进行进一步分析。这样可以减少 SNP 之间的相关性，同时保留具有最强统计证据的 SNP。

附1：OR值与log OR值的转换

下载PGC数据库关于精神分裂症的数据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值（计算PRS需要用OR），可以用R进行转换：

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可以打开任何一个文件查看，什么后缀的都可以但是！！！！！！排列不整齐！！！！进行后续筛选非常不方便！！！！！！
这个时候sort -o可以帮你

sort -o happy.txt sad.assoc

上面的命令就可以把关联分析得到的sad.assoc文件另存为 happy.txt，然后就可以用EXCEL打开了，尽情查看吧。

参考资料

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