bioinformatics - 如何将所有染色体组合在一个文件中

Question

我下载了 VCF 格式的 1000 个基因组数据（染色体 1 -22）。如何将所有染色体组合在一个文件中？我应该先将所有染色体转换为 plink 二进制文件，然后再执行--bmerge mmerge-list？或者有没有其他方法可以将它们结合起来？请问有什么建议吗？

Answer 1

您可以按照您的建议使用 PLINK。您还可以使用BCFtools：

https://samtools.github.io/bcftools/bcftools.html

具体来说，concat命令：

bcftools concat ALL.chr{1..22}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz -Oz -o  ALL.autosomes.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

如果您使用 PLINK，您可能会遇到 1000 Genomes 的问题，因为它包含与 PLINK 不兼容的多等位基因 SNP。此外，有些 SNP 具有相同的 RS 标识符，这也与 PLINK 不兼容。

您需要通过将多等位基因 SNP 拆分为多个记录并删除具有重复 RS 标识符的记录（或制作新的唯一标识符）来解决这些问题，以使 PLINK 正常工作。

此外，PLINK 二进制 PED 不支持基因型概率。我不记得 1000 Genomes 是否包含此类信息。如果确实如此并且您想保留它，则无法将其转换为二进制 PED，因为基因型概率将被硬调用，请参阅：

https://www.cog-genomics.org/plink2/input

具体来说：

由于 PLINK 1 二进制格式不能表示基因型概率，因此不确定性大于 0.1 的调用通常被视为缺失，其余则被视为硬调用。

因此，如果您打算为输出保留 VCF 格式，我建议您不要使用 PLINK。

编辑

以下是将 VCF 转换为 PLINK 的方法：

要从 VCF 文件构建 PLINK 兼容文件，需要合并或删除重复的位置和 SNP id。在这里，我选择删除所有重复的条目。目录重复 SNP id：

grep -v '^#' <(zcat ALL.chr${chrom}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz) | cut -f 3 | sort | uniq -d > ${chrom}.dups

使用 BCFTools，拆分多等位基因 SNP，并使用 plink 删除在上一步中找到的重复 SNP id：

bcftools norm -d both -m +any -Ob ALL.chr${chrom}.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | plink --bcf /dev/stdin --make-bed --out ${chrom} --allow-extra-chr 0 --memory 6000 --exclude ${chrom}.dups

重要的是，这不是解决将 VCF 转换为 PLINK 问题的唯一方法。例如，您可以为重复的 RS id 唯一地分配标识符。

Answer 2

picard GatherVcfs https://broadinstitute.github.io/picard/command-line-overview.html

将分散操作中的多个 VCF 文件收集到单个 VCF 文件中。输入文件必须按基因组顺序提供，并且不得在重叠位置包含事件。