生信教程 | 基于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分析笔记丨检测不同种群之间由于选择引起的差异信息,群体遗传学经典方法
|
数据采集 存储 索引
转录组分析丨一套完整的操作流程简单案例(上)
转录组分析丨一套完整的操作流程简单案例
|
Perl
技术笔记:samtools统计重测序数据深度depth、depth
技术笔记:samtools统计重测序数据深度depth、depth
653 0
PSMC软件分析群体历史有效群体大小步骤(bcftools+PSMC))
PSMC软件分析群体历史有效群体大小步骤(bcftools+PSMC))
|
搜索推荐 Linux Perl
【专栏】Linux 中有趣的命令:`cowsay` 是 Linux 中一个趣味命令,可在终端创造“会说话的牛”效果
【4月更文挑战第28天】`cowsay` 是 Linux 中一个趣味命令,可在终端创造“会说话的牛”效果。基本用法是输入 `cowsay "text"`,展示带有文字的奶牛形象。使用 `-f` 可换不同牛的样式,`-e` 改变文字颜色。还有 `cowthink` 让牛思考,`cowbell` 添加铃铛声。可与其他命令结合,如 `grep` 或 `sed`,增加终端互动性与趣味性。不论新手还是老手,`cowsay` 都能为你的终端带来更多乐趣和个性化。
402 0
|
存储 算法 Linux
算法丨根据基因型VCF文件自动识别变异位点并生成序列fasta文件,基于R语言tidyverse
算法丨根据基因型VCF文件自动识别变异位点并生成序列fasta文件,基于R语言tidyverse
|
存储 数据可视化 数据挖掘
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
|
存储 NoSQL Linux
ArcGIS:如何创建地理数据库、创建要素类数据集、导入要素类、表?
ArcGIS:如何创建地理数据库、创建要素类数据集、导入要素类、表?
9421 1
|
机器学习/深度学习 编解码 人工智能
计算机视觉算法与应用汇总
计算机视觉算法与应用汇总
638 0
|
数据可视化
R可视化学习-- (8) 环形柱状图
柱状图可谓是论文最常见的图表了,但有时候这种笔直排列一排的柱状图用腻了, 换种呈现方式(比如拐个弯)也可换个心情嘛,今天就先来介绍俩种构建环形柱状图的方法,分别基于ggplot2和其他的R包来绘制。
569 0