群体遗传学亲缘关系分析
亲缘关系分析实操
前期准备
给标记加上ID
SNP data通常都是以VCF格式文件呈现，拿到VCF文件的第一件事情就是添加各个SNP位点的ID。
先看一下最开始生成的VCF文件：
原始VCF文件
可以看到，ID列都是"."，需要我们自己加上去。我用的是某不知名大神写好的perl脚本，可以去我的github上下载，用法：
perl path2file/VCF_add_id.pl YourDataName.vcf YourDataName-id.vcf`
当然也可以用excel手工添加。添加后的文件如下图所示（格式：CHROMID__POS）：
添加ID后VCF文件
SNP位点过滤（Missing rate and maf filtering）
SNP位点过滤前需要问自己一个问题，我的数据需要过滤吗？
一般要看后期是否做关联分析（GWAS）；如果只是单纯研究群体结构建议不过滤，因为过滤掉低频位点可能会改变某些样本之间的关系；如果需要和表型联系其来做关联分析，那么建议过滤，因为在后期分析中低频位点是不在考虑范围内的，需要保持前后一致。
如果过滤，此处用到强大的plink软件，用法：
plink --vcf YourDataName-id.vcf --maf 0.05 --geno 0.2 --recode vcf-iid -out YourDataName-id-maf0.05 --allow-extra-chr
参数解释：--maf 0.05：过滤掉次等位基因频率低于0.05的位点；--geno 0.2：过滤掉有2%的样品缺失的SNP位点；--allow-extra-chr：我的参考数据是Contig级别的，个数比常见分析所用的染色体多太多，所以需要加上此参数。
格式转换
将vcf文件转换为bed格式文件。
这里注意一点！！！！：应该是软件的问题，需要把染色体/contig名称变成连续的数字（1 to n），不然会报错无法算出结果！（坑）
plink --vcf YourDataName-id-maf0.05.vcf --make-bed --out snp --chr-set 29 no-xy
参数解释：--chr-set 给出染色体/contig的数目；no-xy 没有xy染色体。
用gcta做亲缘关系分析
gcta输出grm阵列（genetic relationship matrix）
gcta64 --make-grm-gz --out snp.gcta --bfile snp --autosome-num 29
参数解释：--autosome-num常染色体数目。
snp.gcta结果文件：
snp.gcta结果文件
解读：第一，第二列为样品编号；第三列为两样品间有多少个有效位点；第四列为两样品间的亲缘关系的值。
将上述阵列转化为矩阵形式
snp.gcta结果文件列举了两个样品间的关系，我们需要把它变成常见的矩阵形式，这里用R可以轻松完成，写好的R包我放在了Github中：GRM2normal_format.R，大家自行下载使用。
如果不想下载，可以复制如下代码：
library(reshape2)
tmp <- read.table(gzfile("snp.gcta.grm.gz"), header = F, stringsAsFactors = F)
ids <- read.table("snp.gcta.grm.id", header = F, stringsAsFactors = F)
tmp <- tmp[,c(1,2,4)]
result_matrix <- acast(tmp, V1~V2, value.var="V4", drop = F)
makeSymm <- function(m) {
m[upper.tri(m)] <- t(m)[upper.tri(m)]
return(m)
}
result_full <- makeSymm(result_matrix)
result_df <- as.data.frame(result_full)
row.names(result_df) <- ids$V2
colnames(result_df) <- ids$V2
write.table(result_df, file = "ldak.weight.kinship.txt", row.names = T, col.names = NA, sep = "\t", quote = F)
需要用到上诉步骤生成的2个结果文件：snp.gcta.grm.gz和snp.gcta.grm.id。
转换后的结果文件：
gcta.kinship.txt
解读：第一列和第一行都是对应的样品名称。
用LDAK做亲缘关系分析
相比gcta，能用LD对结果进行校正，具体来说，就是先用LD计算每个SNP位点的权重，根据权重再计算Kinship，这样的结果更接近真实情况。
LDAK输出grm阵列（genetic relationship matrix）
- 在不考虑权重的情况下，方法如下：
ldak5.linux --calc-kins-direct snp.ldak --bfile snp --ignore-weights YES --kinship-gz YES --power -0.25
- 用LD计算每个SNP位点的权重，根据权重再计算Kinship
#切割
ldak5.linux --cut-weights snp.sections --bfile snp
#查看有多少个section
cat snp.sections/section.number
#根据自己的section个数分别计算权重（我这里是31个）
for section in {1..31}; do ldak5.linux --calc-weights snp.sections --bfile snp --section $section; done
#weight文件整合，给SNP赋权重值
ldak5.linux --join-weights snp.sections --bfile snp
#输出grm阵列
ldak5.linux --calc-kins-direct snp.ldak.weight --bfile snp --weights snp.sections/weights.all --kinship-gz YES --power -0.25
结果文件和gcta类似，包含两个文件：snp.ldak.weight.grm.gz和snp.ldak.weight.grm.id，用上述同一个R包，同样的方法转化为矩阵形式。
同样的，如果不想下载，直接复制如下代码：
library(reshape2)
tmp <- read.table(gzfile("snp.ldak.weight.grm.gz"), header = F, stringsAsFactors = F)
ids <- read.table("snp.ldak.weight.grm.id", header = F, stringsAsFactors = F)
tmp <- tmp[,c(1,2,4)]
result_matrix <- acast(tmp, V1~V2, value.var="V4", drop = F)
makeSymm <- function(m) {
m[upper.tri(m)] <- t(m)[upper.tri(m)]
return(m)
}
result_full <- makeSymm(result_matrix)
result_df <- as.data.frame(result_full)
row.names(result_df) <- ids$V2
colnames(result_df) <- ids$V2
write.table(result_df, file = "ldak.weight.kinship.txt", row.names = T, col.names = NA, sep = "\t", quote = F)
数据可视化
用R画热图即可，各种热图的画法，我另外找个时间再详细说明，先直接分享一下我作图的代码：
library(pheatmap)
kinship <- read.table("ldak.weight.kinship.txt", header = F, row.names = 1, skip=1)
colnames(kinship) <- row.names(kinship)
diag(kinship) <- NA
hist_data <- hist(as.matrix(kinship), xlab = "Kinship", col = "red", main = "Histogram of Kinship")
pheatmap(kinship, fontsize_row = 0.3, fontsize_col = 0.3)
color = colorRampPalette(c("white", "red","red4"),bias=0.5)(500)
#调整cell大小
par(mar=c(10,10,10,10))
pheatmap(kinship,color=color,border_color = F,fontsize_row = 0.3, fontsize_col = 0.3,cellwidth = 2,cellheight = 2)
#调整聚类树高度
pheatmap(kinship,color=color,border_color = F,fontsize_row = 0.3, fontsize_col = 0.3,cellwidth = 2,cellheight = 2,treeheight_col= 40,treeheight_row = 40)
