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 文件。一些参数解释：

samtools：
- mpileup中的-Q和-q分别确定baseQ和mapQ的截止值
- -v 告诉 mpileup 生成 vcf 输出，-u 表示应该解压缩
- -f 是使用的参考fasta（需要建立索引）
- -r 是调用 mpileup 的区域（在本例中，是基于数组任务 id 的特定染色体）
- P964.bam是要使用的bam文件
bcftools：
- call -c 使用原始调用方法从 mpileup call consensus 序列
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

微信扫一扫分享文章