PRS多基因评分教程学习笔记(二）
之前学习了Base Data质控过程，下面继续，最近一直没有开启博客写作，十月将过，加紧补点。
2、Target Data质控
教程的Target数据来自于基于1000个基因组计划欧洲样本的模拟数据。 目标数据包括通常在您的实验室/部门/协作中生成的个人级别的基因型-表型数据（看来还是需要原始数据的，由于多数数据不公开，涉及遗传资源保护，有的数据用就不错了）。对于本教程，我们使用1000个Genomes Project欧洲样本对一些基因型-表型数据进行了模拟。您可以在此处下载数据，也可以使用以下命令下载数据：
目标数据的质控清单：
质量检查清单的质量检查步骤。 – Sample size 我们建议用户仅对至少100个人的目标数据执行PRS分析。我们目标数据的样本量为503个人。 – File transfer 这条是使用md5值确保数据完整性的，如果不知道md5可以搜索下。 – Genome build 基因组版本要和参考数据一致。 – Standard GWAS QC（标准的GWAS质控） 目标数据的质量必须至少达到GWAS研究中实施的标准，例如 从Hardy-Weinberg平衡中去除具有低基因分型率，低等位基因频率的SNP，去除具有低基因分型率的个体（参见Marees等）。
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
通常，我们可以使用新的样本列表生成新的基因型文件。但是，这将占用大量存储空间。使用plink的–extract，–exclude，–keep，–remove，–make-just-fam和–write-snplist函数，我们可以仅处理样品和SNP列表，而无需复制基因型文件， 减少存储空间的使用。 个体中很高或很低的杂合率可能是由于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)- 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