亲缘关系分析实操

前期准备

给标记加上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

ids

tmp

result_matrix

makeSymm

m[upper.tri(m)]

return(m)

}

result_full

result_df

row.names(result_df)

colnames(result_df)

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

ids

tmp

result_matrix

makeSymm

m[upper.tri(m)]

return(m)

}

result_full

result_df

row.names(result_df)

colnames(result_df)

write.table(result_df, file = "ldak.weight.kinship.txt", row.names = T, col.names = NA, sep = "\t", quote = F)

数据可视化

用R画热图即可，各种热图的画法，我另外找个时间再详细说明，先直接分享一下我作图的代码：

library(pheatmap)

kinship

colnames(kinship)

diag(kinship)

hist_data

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)