基于变异和人类参考构建DNA序列答案

【问题标题】：construct DNA sequence based on variation and human reference基于变异和人类参考构建DNA序列
【发布时间】：2013-09-17 14:15:52
【问题描述】：

千人基因组计划为我们提供了有关数千人 DNA 序列相对于人类参考 DNA 序列“变异”的信息。变体存储在VCF 文件中
格式。基本上，对于该项目中的每个人，我们都可以从 VCF 文件中获取他/她的 DNA 变异信息，例如，变异的类型（例如插入/删除和 SNP）以及变异相对于参考的位置。参考采用 FASTA 格式。通过结合 VCF 文件中一个人的变异信息和 FASTA 文件中的人类参考，我想为那个人构建 DNA 序列。

我的问题是：是否已经存在一些工具可以很好地执行任务，或者我必须自己编写脚本。

【问题讨论】：

标签： bioinformatics vcf-variant-call-format

【解决方案1】：

来自 VCFtools 的 perl 脚本 vcf-consensus 似乎与您正在寻找的内容接近：

vcf-consensus  
Apply VCF variants to a fasta file to create consensus sequence.

Usage: cat ref.fa | vcf-consensus [OPTIONS] in.vcf.gz > out.fa
Options:
   -h, -?, --help         This help message.
   -H, --haplotype <int>  Apply only variants for the given haplotype (1,2)
   -s, --sample <name>    If not given, all variants are applied
Examples:
   samtools faidx ref.fa 8:11870-11890 | vcf-consensus in.vcf.gz > out.fa

发布在 Biostar 上的问题 New fasta sequence from reference fasta and variant calls file? 的答案可能也会有所帮助。

【讨论】：

【解决方案2】：

您可以使用 bcftools (https://github.com/samtools/bcftools) 来执行此任务：

bcftools consensus <file.vcf> \
  --fasta-ref <file> \
  --iupac-codes \
  --output <file> \
  --sample <name>

安装 bcftools：

git clone --branch=develop git://github.com/samtools/bcftools.git
git clone --branch=develop git://github.com/samtools/htslib.git
cd htslib && make && cd ..
cd bcftools && make && cd ..
sudo cp bcftools/bcftools /usr/local/bin/

您还可以将 bcftools 共识与 samtools faidx (http://www.htslib.org/) 结合使用，以从 fasta 文件中提取特定间隔。更多信息参见 bcftools 共识：

About:   Create consensus sequence by applying VCF variants to a reference
         fasta file.
Usage:   bcftools consensus [OPTIONS] <file.vcf>
Options:
    -f, --fasta-ref <file>     reference sequence in fasta format
    -H, --haplotype <1|2>      apply variants for the given haplotype
    -i, --iupac-codes          output variants in the form of IUPAC ambiguity codes
    -m, --mask <file>          replace regions with N
    -o, --output <file>        write output to a file [standard output]
    -c, --chain <file>         write a chain file for liftover
    -s, --sample <name>        apply variants of the given sample
Examples:
   # Get the consensus for one region. The fasta header lines are then expected
   # in the form ">chr:from-to".
   samtools faidx ref.fa 8:11870-11890 | bcftools consensus in.vcf.gz > out.fa

【讨论】：

【解决方案3】：

如果您有一个 fasta 参考基因组和一个 bam 文件，您想通过更改 SNP 和 N 将其转换为参考文件，那么任何仍然访问此页面的人，您可以使用 samtools、bcftools 和 vcfutils 尝试这个单行。 pl（初学者ps：samtools和bcftools都可以在计算集群或Linux中编译，如果可以，只需在软件名称前添加每个位置即可；vcfutils已经是bcftools的perl脚本）

samtools mpileup -d8000 -q 20 -Q 10 -uf  REFERENCE.fasta Your_File.bam | bcftools call -c | vcfutils.pl vcf2fq  > OUTPUT.fastq

d, --max-depth == -q, -min-MQ 要使用的比对的最低映射质量 == -Q, --min-BQ 要考虑的碱基的最低碱基质量 == (你当然可以使用不同的值，见http://www.htslib.org/doc/samtools.html)

这会生成一种看起来像 fastq 但不是的奇怪格式，因此您无法使用转换器对其进行转换，但您可以使用以下 sed 命令，我为此输出专门编写了该命令：

sed -i -e '/^+/,/^\@/{/^+/!{/^\@/!d}}; /^+/ d; s/@/>/g' OUTPUT.fastq

最后，确保将您的新 fasta 文件与您的参考文件进行比较，以确保一切正常。

请谨慎使用 SED 命令，它可能会在质量评分不同的情况下删除您的一些阅读内容

【讨论】：