生信教程 | 基于PSMC估计有效群体大小

简介: 生信教程 | 基于PSMC估计有效群体大小

简介

PSMC 模型使用单个个体的完整二倍体序列中的信息来推断种群规模变化的历史。它最初于 2011 年发布,现已成为基因组学领域非常流行的工具。在本教程中,我们将逐步完成为 PSMC 生成必要的输入数据的步骤,并在发布的猛犸象数据上运行它。

数据

Genome: https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001905.1/

Bam: https://www.ebi.ac.uk/ena/browser/view/ERX935618

这些数据最初是从 Broad 研究所(大象参考基因组)和 ENA( bam 文件)下载的。如果您自己下载数据,则需要在开始之前使用 samtools 索引 fasta 文件和 BAM 文件。

请注意,对于此分析,我们从 BAM 文件开始,其中包含已映射到参考基因组(在本例中为大象)的读数。要在您自己的数据上运行 PSMC,您需要首先将您的读数映射到参考基因组,然后再调整这些脚本。

Install


conda create -n psmc  -c bioconda psmc samtools bcftools

conda activate psmc

索引数据

# genome
samtools faidx loxAfr4.fa 

# bam
samtools index P964.bam

Call consensus 序列

从映射读数开始,第一步是生成 FASTQ 格式的一致序列。为此,我们将使用 samtools/bcftools 工具,遵循论文中描述的方法。

生成consensus序列背后的基本思想是首先使用 samtools mpileup 获取映射读取并生成 VCF 文件。然后,bcftools 使用原始共识调用模型生成consensus序列,并通过 vcfutils.pl 转换为 fastq(带有一些额外的过滤)。

  • 由于 Palkopoulou 等人仅分析了常染色体,因此我们将做同样的事情,依赖于参考文献中 27 个常染色体被命名为 chr1 - chr27 。
samtools mpileup -Q 30 -q 30 -u -v -f loxAfr4.fa -r $CHR P964.bam | bcftools call -c |  \
vcfutils.pl vcf2fq -d 5 -D 34 -Q 30 > P964.$CHR.fq

# $CHR: chr1 - chr27

这将对齐的 bam 文件和参考基因组作为输入,使用 samtools 生成 mpileup,使用 bcftools call consensus序列,然后过滤并将共有序列转换为 fastq 格式,将每个染色体的结果写入单独的 fastq 文件。一些参数解释:

  1. samtools:

    • mpileup中的-Q和-q分别确定baseQ和mapQ的截止值
    • -v 告诉 mpileup 生成 vcf 输出,-u 表示应该解压缩
    • -f 是使用的参考fasta(需要建立索引)
    • -r 是调用 mpileup 的区域(在本例中,是基于数组任务 id 的特定染色体)
    • P964.bam是要使用的bam文件
  2. bcftools:

    • call -c 使用原始调用方法从 mpileup call consensus 序列
  3. vcfutils.pl:

    • -d 5 和 -d 34 确定允许 vcf2fq 的最小和最大覆盖范围,该范围之外的任何内容都会被过滤
    • -Q 30 将均方根映射质量最小值设置为 30

PSMC

PSMC 使用 consensus fastq 文件,并推断种群规模的历史。尽管需要多种参数来控制模型拟合的细节,但我们将遵循 Palkopoulou 等人的做法并使用默认值。

我们需要做的第一件事是将所有单染色体 fastq 文件合并到一个consensus序列中,我们将使用 unix 工具 cat 来完成此操作。

cat P964.chr*.fq > P964.consensus.fq

现在我们需要将此 fastq 文件转换为 PSMC 的输入格式:

$PSMC_HOME/utils/fq2psmcfa P964.consensus.fq > P964.psmcfa

然后我们可以使用默认选项运行 PSMC——但请注意,我们指定 -p 参数,因为论文中报告的默认值与当前默认值不同。

psmc -p "4+25*2+4+6" -o P964.psmc P964.psmcfa

最后,我们使用论文中报告的每代突变率 -u 和以年为单位的世代时间 -g 绘制 PSMC 图。因为论文没有给出他们如何绘制绘图的确切参数,所以这可能看起来与图有点不同,但它会非常接近。

$PSMC_HOME/utils/psmc_plot.pl -u 3.83e-08 -g 31 -p P964_plot P964.psmc
相关文章
XP-CLR分析笔记丨检测不同种群之间由于选择引起的差异信息,群体遗传学经典方法
XP-CLR分析笔记丨检测不同种群之间由于选择引起的差异信息,群体遗传学经典方法
|
1月前
|
存储 人工智能 自然语言处理
大模型备案攻略—2025全网最新最详细解读版
随着AI技术的发展,大模型备案成为行业热点。本文详解备案所需具体条件与注意事项,涵盖模型功能、适用场景、研制情况、安全评估及备案材料等核心内容,帮助企业全面了解备案流程,规避合规风险,顺利推进产品上线。
|
9月前
|
监控 安全 网络安全
社会工程学:概念、技术与防范
社会工程学是一种利用人性弱点而非技术漏洞来获取敏感信息或进行攻击的策略。常见技术包括钓鱼攻击、预设信任、尾随、垃圾箱搜寻、电话欺诈和社交媒体工程。防范措施包括提高意识、双重验证、物理安全、信息管理和技术防护等。通过综合措施,可以有效降低社会工程学攻击的风险,保护信息安全。
489 10
|
机器学习/深度学习 算法 数据挖掘
算法金 | K-均值、层次、DBSCAN聚类方法解析
**摘要:** 这篇文章介绍了聚类分析的基本概念和几种主要的聚类算法。聚类是无监督学习中用于发现数据内在结构的技术,常用于市场分析、图像分割等场景。K-均值是一种基于划分的算法,简单高效但易受初始值影响;层次聚类包括凝聚和分裂方式,形成层次结构但计算复杂;DBSCAN基于密度,能处理任意形状的簇,但参数选择敏感。文章还讨论了这些算法的优缺点和适用场景,并提供了相关资源链接和Python实现。
389 9
算法金 | K-均值、层次、DBSCAN聚类方法解析
screenfull全屏、退出全屏、指定元素全屏的使用步骤
screenfull全屏、退出全屏、指定元素全屏的使用步骤
824 5
|
Python Windows
Python中4种方法实现 xls 文件转 xlsx
【8月更文挑战第6天】以下是Python中将`xls`文件转换为`xlsx`格式的四种方法:1) 使用`pandas`库,通过读取和重新保存文件实现转换;2) 利用`openpyxl`库加载并复制工作簿内容;3) 结合`xlrd`与`xlwt`读取旧格式并写入新格式;4) 在Windows系统下,采用`win32com`自动化Excel应用完成转换。例如,可将`example.xls`文件转换为`converted.xlsx`。
1797 0
PSMC软件分析群体历史有效群体大小步骤(bcftools+PSMC))
PSMC软件分析群体历史有效群体大小步骤(bcftools+PSMC))
|
存储 数据可视化 数据挖掘
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
|
数据采集 搜索推荐 安全
谷歌最快多久能收录我的网站?
答案是:最快是24小时内。 为什么谷歌的收录时间很重要 在网站建设和优化过程中,站长们都非常关心自己的站点何时能被谷歌搜索引擎收录。 这主要是因为,被谷歌收录意味着您的网站内容可以出现在用户的搜索结果中,从而带来流量和潜在的转化机会。 影响谷歌收录时间的因素 站点质量 谷歌对于网站的质量要求颇高。 拥有高质量原创内容、良好的用户体验和高速的加载速度的网站更容易受到搜索引擎的欢迎,并且可能被更快地收录。
350 0
谷歌最快多久能收录我的网站?
|
Linux Python
SGAT丨利用GAPIT进行GWAS分析的方法
SGAT丨利用GAPIT进行GWAS分析的方法