将汇报所用到的分析脚本和绘图脚本进行整合
1.个体识别基因组相似度热图
# 读取数据文件
data <- read.table("A.txt", header = TRUE, row.names = 1)
# 安装并加载必要的包
library(pheatmap)
library(RColorBrewer)
library(grid) # 用于图形调整
# 使用YlGnBu调色板
color_palette <- colorRampPalette(brewer.pal(9, "YlGnBu"))(100)
# 绘制热图
heatmap <- pheatmap(data,
treeheight_row = 0, # 不绘制行树状图
treeheight_col = 0, # 不绘制列树状图
border_color = NA, # 去除边框
scale = "none", # 不对数据进行标准化
cluster_rows = FALSE, # 行不进行聚类
cluster_cols = FALSE, # 列不进行聚类
fontsize_row = 7, # 行标签字体大小
fontsize_col = 7, # 列标签字体大小
width = 8, # 图的宽度
height = 8, # 图的高度
color = color_palette, # 使用自定义颜色调色板
display_numbers = TRUE, # 显示热图数值
number_color = "black", # 数值颜色为黑色
angle_col = 45) # 将列标签倾斜45度
image.png
2.个体识别深度降低造成准确度下降箱线图和散点图
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# 读取txt文件，假设是以制表符分隔
df = pd.read_csv('1.txt', sep='\t')
# 检查数据是否正确加载
print(df.head())
# 定义测序深度的顺序
depth_order = ['15X', '10X', '5X', '3.5X', '2.5X', '1X']
# 将深度列转换为有序类别，以确保按顺序显示
df['deepth'] = pd.Categorical(df['deepth'], categories=depth_order, ordered=True)
# ------------------ 绘制箱线图 ------------------
plt.figure(figsize=(6, 6))
# 使用Seaborn绘制箱线图，区分pair和unpair
sns.boxplot(data=df, x='deepth', y='score', hue='pair_or_not')
# 设置标题和标签
plt.title('Sequencing Depth vs Similarity Score (Boxplot)', fontsize=16)
plt.xlabel('Sequencing Depth', fontsize=12)
plt.ylabel('Similarity Score', fontsize=12)
# 旋转x轴标签，避免重叠
plt.xticks(rotation=45)
# 将图表保存为pdf格式
plt.savefig('plot1.png')
plt.savefig('plot1.pdf')
# 显示图表
plt.tight_layout()
plt.show()
# ------------------ 绘制散点图 ------------------
plt.figure(figsize=(6, 6))
# 使用Seaborn的scatterplot绘制数据，区分匹配与不匹配个体
sns.scatterplot(data=df, x='deepth', y='score', hue='pair_or_not', style='pair_or_not', s=100)
# 设置标题和标签
plt.title('Sequencing Depth vs Similarity Score (Scatter)', fontsize=16)
plt.xlabel('Sequencing Depth', fontsize=12)
plt.ylabel('Similarity Score', fontsize=12)
# 旋转x轴标签，避免重叠
plt.xticks(rotation=45)
# 将图表保存为pdf格式
plt.savefig('plot2.png')
plt.savefig('plot2.pdf')
# 显示图表
plt.tight_layout()
plt.show()
image.png
3.利用少量位点进行个体识别
3.1 绘制MAF频率分布图(Ry语言）
maf_freq <- read.table("MAF_check.frq", header =TRUE, as.is=T)
# 设置图形比例
par(pin = c(5, 4)) # 第一个值是宽度，第二个值是高度，单位为英寸
hist(maf_freq[,5],breaks=seq(0, 0.5, by = 0.01),main = "MAF distribution", xlab = "MAF")
image.png
3.2 绘制染色体上筛选出来的位点数量散点图
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# 读取CSV文件
df_markers = pd.read_csv('selected_markers.csv')
# 检查数据是否正确加载
print(df_markers.head())
# 将染色体编号作为分类变量，以确保图中只显示存在的编号
df_markers['Chromosome'] = df_markers['Chromosome'].astype(str) # 将染色体编号转换为字符串类型，避免自动填充
# 设置图形大小
plt.figure(figsize=(8, 6))
# 绘制带有线条的点图
sns.lineplot(data=df_markers, x='Chromosome', y='Selected Markers', marker='o', color='black')
# 设置标题和标签
plt.title('Number of Selected Markers per Chromosome', fontsize=16)
plt.xlabel('Chromosome', fontsize=12)
plt.ylabel('Number of SNPs markers', fontsize=12)
# 手动设置x轴刻度，确保所有实际存在的染色体编号都显示
plt.xticks(rotation=45) # 将x轴的标签旋转45度，避免重叠
# 设置y轴的范围，确保视觉效果更好
plt.ylim(0, max(df_markers['Selected Markers']) + 5)
# 将图表保存为pdf格式
plt.savefig('plot1.png')
plt.savefig('plot1.pdf')
# 显示图表
plt.tight_layout()
plt.show()
3.3 绘制MAF和PIC小提琴图
# 导入必要的库
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# 读取数据
file_path = 'allele_freq_with_PIC.txt'
data = pd.read_csv(file_path, delim_whitespace=True)
# 提取MAF和PIC列数据
maf = data['MAF']
pic = data['PIC']
# 设置图形大小
plt.figure(figsize=(6, 6))
# 使用小提琴图绘制MAF和PIC列，调整缩放比例使图形更美观
sns.violinplot(data=[maf, pic], scale='width', inner='quartile')
# 设置X轴标签
plt.xticks([0, 1], ['MAF', 'PIC'])
# 设置标题
plt.title('Violin Plot of MAF and PIC (Adjusted for Aesthetics)')
# 调整布局以确保图形不重叠
plt.tight_layout()
# 保存图片
plt.savefig('plot_violin.pdf')
# 显示图形
plt.show()
image.png
3.4 计算位点的平均间距
# 导入所需的库
import pandas as pd
# 读取CSV文件
file_path = 'selected_snps.csv'
data = pd.read_csv(file_path)
# 将数据按照染色体编号和位点位置进行排序
data = data.sort_values(by=['Chromosome', 'Position'])
# 初始化一个列表用于存储计算的位点距离
distances = []
# 循环处理每条染色体
for chrom in data['Chromosome'].unique():
# 提取当前染色体的所有SNP位点
chrom_data = data[data['Chromosome'] == chrom]
# 计算当前染色体上每个相邻位点之间的距离
for i in range(1, len(chrom_data)):
pos_current = chrom_data.iloc[i]['Position'] # 当前位点的坐标
pos_previous = chrom_data.iloc[i - 1]['Position'] # 前一个位点的坐标
distance = pos_current - pos_previous # 计算两个位点之间的距离
distances.append(distance) # 将距离添加到列表中
# 计算所有位点之间的平均距离（单位为bp）
mean_distance_bp = sum(distances) / len(distances)
# 将平均距离从bp转换为kb
mean_distance_kb = mean_distance_bp / 1000
# 输出平均距离
print(f"位点之间的平均距离为：{mean_distance_kb:.2f} kb")
3.5 比较每个个体之间314个位点之间的相似度
生成一个相似度矩阵
# 导入所需的库
import pandas as pd
# 读取 VCF 文件
file_path = 'filtered315.vcf'
# 跳过 VCF 文件的头部元数据，只读取数据部分
vcf_data = pd.read_csv(file_path, sep='\t', comment='#', header=None)
# 提取样本信息的列，即从第10列开始为每个样本的基因型信息
# 具体情况要看VCF格式，假设前几列为常见VCF格式字段如 CHROM, POS, ID 等
genotype_data = vcf_data.iloc[:, 9:]
# 初始化一个相似度矩阵（用 DataFrame 存储）
num_individuals = genotype_data.shape[1]
similarity_matrix = pd.DataFrame(index=genotype_data.columns, columns=genotype_data.columns)
# 定义一个函数来计算两个个体之间的基因型相似度
def calculate_similarity(genotype1, genotype2):
# 计算两个个体基因型的相似度，简单版本就是直接比较每个位点是否相同
matches = sum(genotype1 == genotype2)
total_snps = len(genotype1)
return matches / total_snps
# 计算相似度矩阵
for i in range(num_individuals):
for j in range(i, num_individuals):
# 提取第i个和第j个个体的基因型数据
genotype1 = genotype_data.iloc[:, i]
genotype2 = genotype_data.iloc[:, j]
# 计算相似度
similarity = calculate_similarity(genotype1, genotype2)
# 将相似度填入矩阵
similarity_matrix.iloc[i, j] = similarity
similarity_matrix.iloc[j, i] = similarity # 对称矩阵
# 输出相似度矩阵
print(similarity_matrix)
# 你可以将矩阵保存到文件
similarity_matrix.to_csv('similarity_matrix.csv')
3.6 根据生成的csv文件，生成热图
# 导入必要的库
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
# 读取相似度矩阵CSV文件
file_path = 'similarity_matrix.csv'
similarity_matrix = pd.read_csv(file_path, index_col=0)
# 设置图形大小
plt.figure(figsize=(10, 8))
# 使用Seaborn的heatmap函数绘制热图
sns.heatmap(similarity_matrix, annot=False, cmap='coolwarm', cbar=True)
#修改标签文字大小和方向
plt.xticks(rotation=45, ha='right', fontsize=6, color='#515151') # 进一步缩小字体到6
plt.yticks(fontsize=6, color='#515151')
# 调整布局以确保图形不重叠
plt.tight_layout()
# 保存为pdf
plt.savefig('hotmap_individuals_similiar.pdf')
image.png
4.溯源工作
4.1 绘制Rped与Rgenome相关性曲线图
# 导入必要的库
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from scipy.stats import linregress
# 加载数据集
data = pd.read_excel('Rped_Rgenome.xlsx')
# 使用linregress计算线性回归模型，包含斜率、截距、相关系数r值、p值和标准误差
regression_result = linregress(data['Rped'], data['Rgenome'])
# 计算R²值，它是相关系数r值的平方
r_squared = regression_result.rvalue ** 2
# 设置整体图表的美观风格为白色背景，以保持图表的清洁
sns.set_style("white")
# 创建图表，设置图表的大小为6x6英寸，使图表在水平方向上更加紧凑
plt.figure(figsize=(6, 6))
# 生成散点图和回归线
# 使用scatter_kws设置散点图点的外观，如透明度alpha=0.6和颜色#6fc7b2
# 使用line_kws设置回归线的颜色为#14a885
# ci=95设置95%的置信区间，颜色为灰色，用于绘制回归线周围的阴影区域
ax = sns.regplot(x='Rped', y='Rgenome', data=data, scatter_kws={'alpha':0.6, 'color': '#6fc7b2'},
line_kws={'color': '#14a885'}, ci=95, color='gray')
# 移除顶部和右侧的边框，以清理图表的外观
sns.despine()
# 加强左侧和底部边框的外观，通过设置边框宽度为2，使其更加醒目
ax.spines['left'].set_linewidth(2)
ax.spines['bottom'].set_linewidth(2)
# 调整刻度参数，使其更大更粗，提高可读性
# axis='both'表示同时调整x轴和y轴，which='major'表示调整主刻度，labelsize=12设置标签大小，width=2设置刻度线宽度
ax.tick_params(axis='both', which='major', labelsize=12, width=2)
# 在图表上添加R²值和p值作为参考
# 使用plt.text在指定位置添加文本，文本内容包括R²值和p值，fontsize=12设置字体大小
plt.text(0.00045, 0.4, f'R²={r_squared:.2f}, p={regression_result.pvalue:.2e}', fontsize=12)
# 设置图表的标题和轴标签，并增强字体大小和重量，以便更好的可视化
# plt.title设置标题，plt.xlabel和plt.ylabel分别设置x轴和y轴的标签，fontsize和fontweight设置字体大小和加粗
plt.title('Rped vs Rgenome with Regression Line', fontsize=14, fontweight='bold')
plt.xlabel('Rped', fontsize=12, fontweight='bold')
plt.ylabel('Rgenome', fontsize=12, fontweight='bold')
# 将更加紧凑的图形保存为PDF文件
output_file = 'Rped_Rgenome.pdf'
plt.savefig(output_file, format='pdf') # 指定保存格式为PDF
# 显示图表
plt.show()
image.png
4.2 绘制随着snp位点数量的变化，R2（Rped与Rgenome）变化的折线图
# 导入必要的库
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# 读取Excel文件
file_path = 'relation.xlsx'
data = pd.read_excel(file_path)
# 设置绘图风格和图形大小
sns.set_style("whitegrid") # 使用白色网格背景
plt.figure(figsize=(8, 6)) # 设置图像大小为8x6英寸
# 使用lineplot绘制线条图，loci作为x轴，R²作为y轴，并在每个数据点显示标记
sns.lineplot(x='loci', y='R2', data=data, marker='o', color='#6fc7b2')
# 移除网格线
plt.grid(False)
# 移除上方和右侧的轴线
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
# 增加坐标轴线的粗细
plt.gca().spines['bottom'].set_linewidth(2) # x轴线
plt.gca().spines['left'].set_linewidth(2) # y轴线
# 设置x轴的刻度标签为实际的SNP位点数量，并旋转以防止重叠
plt.xticks(ticks=data['loci'], labels=[f'{int(l)}' for l in data['loci']], rotation=45)
# 设置图表标题和轴标签
#plt.title('R²值随SNP位点数量的变化趋势', fontsize=16, fontweight='bold') # 图表标题
plt.xlabel('SNP Loci Number', fontsize=14) # x轴标签
plt.ylabel('R² (Correlation Coefficient)', fontsize=14) # y轴标签
# 将x轴反转，使SNP位点数量从高到低显示
plt.gca().invert_xaxis()
# 将更加紧凑的图形保存为PDF文件
output_file = 'SNP_number_Rped_Rgenome.pdf'
plt.savefig(output_file, format='pdf') # 指定保存格式为PDF
# 调整布局以防止标签和图形重叠
plt.tight_layout()
# 显示图表
plt.show()
image.png
5.判断各种亲缘关系的亲缘系数的取值范围
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
# 加载数据
file_path = 'SNP0.02fig.xlsx' # 替换为您的文件路径
data = pd.read_excel(file_path)
# 将亲缘系数值保留两位小数
data['correlation_index'] = data['correlation_index'].round(2)
# 检查并转换数据类型
data['correlation_index'] = pd.to_numeric(data['correlation_index'], errors='coerce')
data = data.dropna(subset=['correlation_index'])
# 提取独特的关系类型
relationships = data['relationship'].unique()
# 准备绘制正态分布曲线
plt.figure(figsize=(10, 6))
colors = ['blue', 'green', 'red', 'purple']
for i, relationship in enumerate(relationships):
subset = data[data['relationship'] == relationship]
mean = subset['correlation_index'].mean()
std_dev = subset['correlation_index'].std()
# 生成正态分布数据
x = np.linspace(mean - 3*std_dev, mean + 3*std_dev, 1000)
y = norm.pdf(x, mean, std_dev)
# 绘制正态分布曲线
plt.plot(x, y, color=colors[i], label=f'{relationship} (µ={mean:.2f}, σ²={std_dev**2:.2f})')
# 美化图表
plt.title('Normal Distribution of Correlation Indexes by Relationship', fontsize=14, fontweight='bold')
plt.xlabel('Correlation Index', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.legend(title='Relationship', title_fontsize='13', fontsize='11')
plt.grid(True)
# 保存图片
plt.savefig('relationship_kinship.png')
plt.savefig('relationship_kinship.pdf')
# 显示图表
plt.show()
image.png
