【基于R语言群体遗传学】-5-扩展到两个以上等位基因及多基因位点

我们现在继续对于群体遗传学进行统计建模,书接上回,我们讨论了孤雌生殖的物种违反哈代温伯格遗传比例的例子,那我们现在来看多于两个等位基因的情况的计算。

如果没有看过之前文章的同学,可以先去看一下之前的文章:

群体遗传学_tRNA做科研的博客-CSDN博客

多等位基因情况

到目前为止,我们一直专注于一个双等位基因(bi-allelic)系统,其中一个等位基因的频率表示为p,另一个等位基因的频率表示为(1-p)。然而,基因型预测可以很容易地扩展到超过两个等位基因。由于我们总是假设完全随机交配,我们理论上预期的纯合子数量仍然将等于p²,无论我们考虑多少个等位基因。一组j个等位基因的预期纯合率会是每个等位基因频率平方的和

p1 <- 0.2
p2 <- 0.3
p3 <- 0.5
sum(sapply(c(p1,p2,p3),function(x) x^2))

假设我们有三个等位基因频率,p1、p2和p3,并希望计算总体预期的纯合子基因型频率

这给了我们一个总体的纯合子基因型频率为38%。然后我们可以相当容易地得到一个预期的杂合子频率:

多位点情况

接下来让我们读取一个来自东地中海地区阿勒颇松(Pinus halepensis)的真实多采样等位基因数据集,整理这些数据,并计算我们预期的以及观察到的杂合子频率(改编自Gershberg等人,2016)。使用popgenr包,安装请看前面的博客:

基因型中的数据都来自微卫星。微卫星是由短重复的核苷酸组成特征的DNA片段。这类遗传标记因其高度变异(即具有高突变率)而常用于研究。让我们用str()函数来了解一下数据:

library("popgenr")
data(genotypes)
str(genotypes)

我们可以看到这个数据框有181个观测值(行)和20列。其中十八列代表不同的位点,每个独特的数字是一个不同的等位基因,前两列是每个样本的个体ID($ID)和种群分配($Pop)。我们想要编写一个迭代代码来遍历数据集进行计算;为了简化这个过程,我们先对数据进行一些简化处理。

rownames(genotypes) <- genotypes$ID
genotypes <- genotypes[,-c(1,2)]

根据这个数据集的设置,每个个体有两列,代表在一棵采样的二倍体松树中一个位点的两个拷贝。让我们计算数据集中位点的总数。对数据框使用length()函数应该只返回列的数量。由于每个位点由两列表示,我们可以将其除以2得到采样的位点总数: 

(num.loci <- (length(genotypes))/2)

现在我们知道我们正在处理九个不同的位点。接下来我们需要弄清楚每个位点实际上有多少个等位基因。我们将使用for循环来为每个位点进行等位基因计数。、

Hom_exp <- NULL
Het_exp <- NULL
Hom_obs <- NULL
Het_obs <- NULL
for(n in 1:(num.loci)){ # 对于每一个基因座current <- n*2-1 # 计算当前基因座的起始位置locus <- c(genotypes[,current],genotypes[,current+1]) # 获取当前基因座的两个等位基因alleles <- unique(locus) # 获取该基因座的所有独特等位基因alleles <- alleles[alleles!=-1] # 移除非等位基因标记(例如缺失数据标记为-1)p_allele <- NULL # 初始化等位基因频率向量for(a in 1:length(alleles)){ # 对于每个独特的等位基因p_allele <- c(p_allele,sum(alleles[a]==locus)/sum(locus!=-1)) # 计算等位基因频率}Hom_exp <- c(Hom_exp, sum(sapply(p_allele,function(x) x^2))) # 期望纯合子频率是每个等位基因频率的平方和obs <- 0 # 初始化观察到的纯合子计数for(i in 1:length(genotypes[,current])){ # 对于当前基因座的每个个体if(genotypes[i, current]!=-1){ # 如果等位基因不是缺失数据if(genotypes[i, current]==genotypes[i,current+1]){ # 如果两个等位基因相同obs <- obs+1 # 增加纯合子计数}}}Hom_obs <- c(Hom_obs,obs/(sum(locus!=-1)/2)) # 观察到的纯合子频率是纯合子计数除以有效等位基因总数的一半
}

在R语言中,极少数会使用繁琐的循环,但是为了处理这里的数据我们不得不这样做,所以这也是为什么我们课题组会使用Java进行计算,我简述一下这个代码的意思:我们写一个循环遍历每个基因座(locus)--内部循环计算每个等位基因的频率--计算期望纯合子频率--计算观察到的纯合子频率:

现在要找到预期和观察到的杂合子频率,我们只需从频率的总和中减去我们的纯合子频率:

Het_exp <- 1- Hom_exp
Het_obs <- 1- Hom_obs

让我们绘制这九个位点的观察到的杂合度频率与预期杂合度频率的对比图,看看它们之间的关系如何。我们将从简单地使用plot绘制Het_obs和Het_exp开始,然后绘制一条回归线。我们将使用lm(线性模型)函数来估计数据集之间的线性关系(最小二乘线性回归)

# 绘制观测值与期望值的散点图
plot(Het_obs, Het_exp)# 添加线性回归线
abline(lm(Het_exp ~ Het_obs))# 进行线性回归分析并打印摘要
reg <- summary(lm(Het_exp ~ Het_obs))
print(reg)# 提取并打印决定系数(r-squared)
rr <- reg$r.squared
rrlabel <- paste("r-squared =", round(rr, digits = 3))
text(0.6, 0.2, rrlabel)# 提取并打印P值
pv <- reg$coefficients[2, 4]
pvlabel <- paste("P-value =", pv)
text(0.6, 0.15, pvlabel)

查看我们的图,我们应该看到一个很好的验证,即哈代-温伯格预测可以扩展到多个位点,在这种情况下,它在预测广泛的杂合性值方面表现得相当好。 

下一篇博客我们将讲述血液型及血液等位基因频率的内容,欢迎大家点赞关注!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/diannao/40171.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

开源租房项目

项目名称项目地址描述体验地址后端代码前端代码小程序端代码gitHubstart租房或房屋交易项目https://github.com/saysky/manland?tabreadme-ov-filePC端 管理端http://manland.liuyanzhao.com/有有无房适–房屋租赁管理平台https://github.com/LiuXIn011/rightHouse开源房屋管理…

非对称加密算法原理与应用1——秘钥的生成

作者:私语茶馆 1.前言 非对称算法有非常多的用途,实现license管控,数字签名,加密内容等等,由于涉及场景和标准非常多,因此实际使用过程中还是存在一定门槛,这里记录一下利用非对称算法RSA的应用关键点,并提供实现license管理的案例。预计拆分为以下几个章节: (1)秘…

Apipost接口测试工具的原理及应用详解(三)

本系列文章简介: 随着软件行业的快速发展,API(应用程序编程接口)作为不同软件组件之间通信的桥梁,其重要性日益凸显。API的质量直接关系到软件系统的稳定性、性能和用户体验。因此,对API进行严格的测试成为软件开发过程中不可或缺的一环。在众多API测试工具中,Apipost凭…

【分布式数据仓库Hive】HivQL的使用

目录 一、Hive的基本操作 1. 使用Hive创建数据库test 2. 检索数据库&#xff08;模糊查看&#xff09;&#xff0c;检索形如’te*’的数据库 3. 查看数据库test详情 4. 删除数据库test 5. 创建一个学生数据库Stus&#xff0c;在其中创建一个内部表Student&#xff0c;表格…

ubuntu20.04在anaconda环境下不能使用catkin_make

ubuntu20.04在anaconda环境下不能直接使用catkin_make编译&#xff0c;报错显示需要安装python3-empy 这时候查询会发现该软件包已经安装了&#xff0c;但是是在ROS环境中&#xff0c;安装anaconda环境后python解释器的指向变了&#xff0c;所以需要在anaconda环境中再装pytho…

Unity udp通信详解

在Unity中实现UDP通信&#xff0c;需要使用C#的System.Net和System.Net.Sockets命名空间。UDP&#xff08;用户数据报协议&#xff09;是一种无连接的网络协议&#xff0c;它允许数据包在网络上发送和接收&#xff0c;但不保证数据包的到达顺序、完整性或可靠性。这使得UDP非常…

沃德校园助手丨校园跑腿-校园外卖-校园论坛三合一系统

校园跑腿项目其实由来已久&#xff0c;由于大学校园生活的特殊性&#xff0c;除了日常的课程学习之外&#xff0c;大学生的所有生活基本长期处于全模式形态下的校园封闭环境中&#xff0c;再加之当前大学生一部分学业繁忙&#xff0c;办事不便。另一部分自理能力较差&#xff0…

【kafka】可视化工具cmak(原kafka-manager)安装问题解决

众所周知&#xff08;反正不管你知不知道&#xff09;&#xff0c;kafka-maneger更名了&#xff0c;现在叫cmak&#xff01;原因是什么呢&#xff1f;据不可靠小道信息说&#xff0c;原kafka-manager这个名字涉及到kafka商标使用问题&#xff0c;应该是被律师函警告了&#xff…

如何批量创建、提取和重命名文件夹!!!

你是否还在一个一个手动创建文件名&#xff01; 你是否还在一个一个手动提取文件名&#xff01; 你是否还在一个一个手动修改文件名&#xff01; 请随小生一起批量自动创建、提取、重命名&#xff01; 1、批量创建文件夹 【案例】创建1日-31日共31个文件夹 【第一步】在A列…

Android SurfaceFlinger——动画播放流程(十六)

前两篇文章介绍了系统启动动画服务的启动和准备阶段,并且我们选择了自定义动画的分支,该分支的动画播放流程主要包含一下几个阶段: loadAnimation:解析 zip 包的动画数据。playAnimation:播放解析好的纹理数据。releaseAnimation:播放完毕释放资源。一、动画播放流程 1、…

Gradle学习-5 发布二进制插件

注&#xff1a;以下示例基于Gradle8.0 1、发布插件 复制一分 buildSrc&#xff0c;执行命令行&#xff0c;生成一个新目录 leon-gradle-plugin cp -rf buildSrc leon-gradle-plugin在 leon-gradle-plugin 目录下的 build.gradle 中引入maven plugins{// 引用 Groovy 插件&…

(五十二)第 8 章 动态存储管理(边界标识法)

1. 背景说明 2. 示例代码 1) errorRecord.h // 记录错误宏定义头文件#ifndef ERROR_RECORD_H #define ERROR_RECORD_H#include <stdio.h> #include <string.h> #include <stdint.h>// 从文件路径中提取文件名 #define FILE_NAME(X) strrchr(X, \\) ? strr…

Linux环境下的字节对齐现象

在Linux环境下&#xff0c;字节对齐是指数据在内存中的存储方式。字节对齐是为了提高内存访问的效率和性能。 在Linux中&#xff0c;默认情况下&#xff0c;结构体和数组的成员会进行字节对齐。具体的对齐方式可以通过编译器选项来控制。 在使用C语言编写程序时&#xff0c;可…

【Linux】线程——线程的概念、线程的特点、线程的优点和缺点、线程和进程、线程函数的使用

文章目录 Linux线程1. 线程的概念1.1 什么是线程 2. 线程的特点2.1 线程的优点2.2 线程的缺点2.4 线程和进程 3. 线程函数的使用pthread_create() 创建线程pthread_self() 获取线程IDpthread_exit() 线程终止pthread_cancel() 线程取消pthread_join() 线程等待pthread_detach()…

【第14章】MyBatis-Plus批量操作

文章目录 前言一、功能概览二、类结构说明1.MybatisBatch<?>2.MybatisBatch.Method<?>3. BatchMethod<?>4.使用步骤5.返回值说明 三、使用示例1. execute方法2. 示例一&#xff1a;实体类型数据3. 示例二&#xff1a;非实体类型数据4. 示例三&#xff1a;…

茗鹤 | 如何借助APS高级计划排程系统提高汽车整车制造的效率

在我们做了详尽的市场调研及头部汽车制造企业排程需求沟通后&#xff0c;我们发现尽管企业有很多的业务系统做支撑&#xff0c;在计划排程领域&#xff0c;所有的汽车制造总装厂仍旧使用人工“Excel”做排产规划&#xff0c;其中少部分也会借助MRP、第三方辅助排产工具。鉴于我…

ChatGPT:Java中的对象引用实现方式

ChatGPT&#xff1a;Java中的对象引用实现方式 如果使用句柄的话&#xff0c;那么 Java 堆中将会划分出一块内存来作为句柄池&#xff0c;reference 中存储的就是对象的句柄地址&#xff0c;而句柄中包含了对象实例数据与对象类型数据各自的具体地址信息。 你提到的句柄机制是…

JVM原理(十一):JVM虚拟机六种必需对类进行初始化的情况

Java虚拟机把描述类的数据从Class文件加载到内存&#xff0c;并对数据进行校验、转换解析和初始化&#xff0c;最终形成可以被虚拟机直接使用的Java类型&#xff0c;这个过程被称作虚拟机的类加载机制。Java天生可以动态扩展的语言特性就是依赖运行期间动态加载和动态链接这个特…

104.二叉树的最大深度

给定一个二叉树 root &#xff0c;返回其最大深度。 二叉树的 最大深度 是指从根节点到最远叶子节点的最长路径上的节点数。 示例 1&#xff1a; 输入&#xff1a;root [3,9,20,null,null,15,7] 输出&#xff1a;3 示例 2&#xff1a; 输入&#xff1a;root [1,null,2] 输出…

相机参数与图像处理技术解析

01. 相机内参和外参的含义&#xff1f;如果将图像放大两倍&#xff0c;内外参如何变化&#xff1f; 相机有两个最基础的数据&#xff1a;内参(Instrinsics)和外参(Extrinsics)&#xff0c;内参主要描述的是相机的CCD/CMOS感光片尺寸/分辨率以及光学镜头的系数&#xff0c;外参主…