样本间亲缘关系矩阵(kinship matrix)和同源性矩阵(IBS matrix)构建的方式
1. 可以使用plink的–make-rel计算个体之间的亲缘关系(强调个体之间的遗传相似性)
/opt/software/plink --bfile vcf_bfile--make-rel --out relatedness_matrix # 得到亲缘关系距离矩阵:
# relatedness_matrix.rel
2. kinship
# 利用tassel计算
run_pipeline.pl -Xmx1536m-Xms512m -SortGenotypeFilePlugin -inputFile 你的vcf文件 -outputFile outvcf -fileType VCF
run_pipeline.pl-Xmx1536m -Xms512m -importGuess outvcf -KinshipPlugin -methodCentered_IBS -endPlugin -export tassel_kinship.txt -exportType SqrMatrix
# 利用gcta计算
使用 --make-grm-alg 1 或 --make-grm 0
gcta --make-grm --make-grm-alg 1 --out snp.gcta --bfile vcf_bfile snp --autosome-num 90
3. IBS
/opt/software/plink --bfile vcf_bfile --make-bed --out IBS_matrix --maf 0.05 --recode --double-id --allow-extra-chr --chr-set 90 --distance square ibs
要计算遗传距离,使用1-ibs
群体关系矩阵如何构建?
转换方法:平均IBS(1-IBS)(个体对间均值)
计算所有个体两两之间的IBS(1-IBS)均值,反映群体内遗传相似性、遗传距离。
python3 /home/huguang/script/GWAS/group_dis.py distance_matrix group_labels # 具体group_dis.py代码如下
import sys
import numpy as np
import pandas as pd
from itertools import combinationsdef compute_intergroup_distances(distance_matrix, group_labels):"""计算组间遗传距离矩阵(组间平均距离)参数:distance_matrix: 全样本的遗传距离矩阵 (n x n)group_labels: 样本分组标签列表 (length = n),例如 [0,0,1,1,2,2]返回:group_dist_matrix: 组间距离矩阵 (k x k),k为组数"""# 参数检查assert distance_matrix.shape[0] == distance_matrix.shape[1], "距离矩阵必须是方阵"assert len(group_labels) == distance_matrix.shape[0], "标签数与矩阵维度不匹配"groups = np.unique(group_labels)k = len(groups)group_dist_matrix = np.zeros((k, k))# 计算每对组间的平均距离for i, j in combinations(range(k), 2):# 获取组i和组j的样本索引idx_i = np.where(group_labels == groups[i])[0]idx_j = np.where(group_labels == groups[j])[0]# 提取子矩阵submatrix = distance_matrix[np.ix_(idx_i, idx_j)]# 计算平均距离(可替换为中位数等其他统计量)avg_dist = np.nanmean(submatrix)# 对称赋值group_dist_matrix[i, j] = avg_distgroup_dist_matrix[j, i] = avg_dist# 对角线(组内距离)特殊处理for i in range(k):idx = np.where(group_labels == groups[i])[0]submatrix = distance_matrix[np.ix_(idx, idx)]group_dist_matrix[i, i] = np.nanmean(submatrix[np.triu_indices(len(idx), k=1)])return group_dist_matrix# 示例用法
if __name__ == "__main__":full_matrix = pd.read_csv(sys.argv[1], sep='\t', header=None, index_col=None, engine='c').to_numpy()groups = pd.read_csv(sys.argv[2], sep="\t", header=None, index_col=None).iloc[:, 0].to_numpy()# 计算组间距离矩阵result = compute_intergroup_distances(full_matrix, groups)print("组间遗传距离矩阵:")print(np.round(result, 4))