通过VCF文件制作定制化的参考基因组

Posted 2023-04-13

参考技术A 最近做的项目要对 reference genome 基于突变进行一些modify，制作 personalized genome 或者说是 psuedo-genome 伪基因组。其实就是把某个测序样本call出来的 SNP && indel 替换掉参考基因组对应位置的碱基。
自己可以编写脚本修改，使用 Perl 中的 substr 来进行单个位点修改，把坐标按照从后往前的顺序，而不是从前到后。
其实有蛮多工具可以做这个的，在这里安利一下，大家可以使用体验一下。

GATK作为通用的变异检测的软件，其中有很多有用的工具，这里介绍一下 FastaAlternateReferenceMaker : Create an alternative reference by combining a fasta with a vcf.

准备好检测的 VCF 文件，参考基因组FASTA文件，使用如下命令：

这里的 -L 参数可以多单个基因或者基因组一段区域进行替换。
--snp-mask 参数，当构建时 psuedo-genome.fa 时，该VCF文件中的SNP用作掩码（在序列中插入N）。
运行结束，个人的参考基因组就构建好了，一般制作 psuedo-genome.fa 就是为了消除变异带来的影响，部分其他参数可以到gatk官网查阅。另外 FastaAlternateReferenceMaker 使用简单的 indel ，当VCF文件包含复杂的位点时（complex substitutions），会忽略。
PS:

perEditor 用于分析phased SNPs/indels，意味着VCF文件中的 GT 是 0|1/1|1/1|2 ，而不是我们常见的 0/1 这种，明确知道变异是在母本还是父本染色体上，第一个是母本，第二个是父本，使用方式如下：

这里对参数进行一下解释：
allele 两个选项 mother|father ，代表等位基因来自于哪一个亲本；
indivicual ，整数，这个是因为有一些VCF文件中不止包含一个样本的突变信息，参数代表选择第几个样本的突变来进行创建 personalized genome ， 1 代表第一个样本；

perEditor_ra ，允许用户把染色体重排信息考虑到创建个性化参考基因组中，一样只服务于phased的数据，另外还需要把涉及染色体重排的染色体放在一个目录下：

这里的 individual 和 allele 参数和 perEditor 是一样的含义；
rearrangements_t.vcf ，染色体重排注释结果，vcf format version 4.1；
chr_length_t.bed ，文件中包含基因组每一条染色体的长度，即核苷酸总数，格式如下：

centromeres_t.bed ，文件中包含每条染色体的着丝粒位置的文件，着丝粒坐标用于确定新染色体的身份，当两条染色体由于重排而合并时，新染色体将按照着丝粒来源的名称来命名，继承两个着丝粒，其名称则包含两个亲本染色体的名称，文件格式如下：

bed文件和vcf文件位于同一个工作目录下，最终生成的新的染色体命名中添加了 new ：

g2gtools 通过将SNP和indels整合到参考基因组中来创建自定义基因组，提取感兴趣的区域（例如外显子或转录本），并在两个基因组之间转换文件（bam，gtf，bed）的坐标。 与其他 liftover 工具不同， g2gtools 不会丢弃掉落在indel区域上的alignments。 版本0.2可以创建个人二倍体基因组。 新版本仍然将个人基因组坐标上的二倍体比对转换回参考基因组，因此我们可以比较种群样本之间的比对。

先激活virtual environment：

创建自定义的基因组，需要下列信息：

下面是操作流程：

vcf2diploid 通过将phased的变异整合到参考基因组中，从vcf文件创建二倍体基因组。

samtools 、 bcftools 和 vcfutils.pl 三个程序联合使用，从 BAM 文件开始操作，到获得新的参考基因组：

FastaAlternateReferenceMaker
perEditor A tool to create personalized genome sequences
g2gtools creates custom genomes by incorporating SNPs and indels into reference genome
Welcome to g2gtools’ documentation
Personal Genome Constructor (vcf2diploid tool)
construct DNA sequence based on variation and human reference

pysam - 多种格式基因组数据（sam/bam/vcf/bcf/cram/…）读写与处理模块（python）

在开发基因组相关流程或工具时，经常需要读取、处理和创建bam、vcf、bcf文件。目前已经有一些主流的处理此类格式文件的工具，如samtools、picard、vcftools、bcftools，但此类工具集成的大多是标准功能，在编程时如果直接调用的话往往显得不够灵活。

本文介绍的是一个处理基因组数据的python模块，它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能，能在编程时非常灵活的处理bam和bcf文件。

以下主要介绍pysam的安装和使用方法：

1. 安装

如果Linux上安装了pip，可以一键安装，在集群上的话，需要登录安装节点进行安装。

pip3 install pysam

检查是否安装成功

import pysam

 

2.读取bam文件(pysam.AlignmentFile)

bam是sam的二进制文件，因其占用空间少，所以都会使用bam进行存储和操作。

要读取bam文件,必须先创建一个AlignmentFile对象.

path_in = ‘./test.bam‘
samfile = pysam.AlignmentFile(path_in, "rb")

之后就可以逐行读取和处理bam文件了(顺序读取),以下打印出了bam的一行.

for line in samfile:
    print(line)
    break

 

但顺序读取还不够灵活,我们有时需要随机读取(提示:sam不能随机读取),pysam的fetch方法提供了随机读取功能.

直接使用fetch会报错

ValueError: fetch called on bamfile without index

提示我们需要建立(.bai)索引

samtools index corrected.bam

fetch返回的是一个迭代器(iterator),可以迭代读取内容.

for read in samfile.fetch(‘chr6‘, 28478220, 28478222):
...     print(read)

fetch方法的API如下,chr6为参考序列,后面数字分别为读取的起始和终止位置.

fetch(self, reference=None, start=None, end=None, region=None, tid=None, until_eof=False, multiple_iterators=False)

 

3.读取vcf/bcf文件(pysam.VariantFile)

读取方法同上,只是使用的是VariantFile方法:

gvcf = "./MHC.unified.g.vcf.gz"
vcf_in = pysam.VariantFile(gvcf)

若想随机读取,仍然需要建立索引:

首先使用bgzip压缩vcf

bgzip -c MHC.g.vcf > MHC.g.vcf.gz

然后用bcftools建立索引

bcftools index -c MHC.g.vcf.gz

使用fetch读取

for rec in vcf_in.fetch(‘chr6‘, 28577796, 28577896):
...     print(rec)
...     break

 

4.创建并写入到新的bam或vcf文件

pysam的核心功能是可以随心所欲的读取数据,处理之后,写入到一个新建的bam或bcf文件里.

我们可以完全自定义一些内容,然后写入到一个新的bam文件里,如下:

header = { ‘HD‘: {‘VN‘: ‘1.0‘},
            ‘SQ‘: [{‘LN‘: 1575, ‘SN‘: ‘chr1‘},
                   {‘LN‘: 1584, ‘SN‘: ‘chr2‘}] }

with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:
    a = pysam.AlignedSegment()
    a.query_name = "read_28833_29006_6945"
    a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
    a.flag = 99
    a.reference_id = 0
    a.reference_start = 32
    a.mapping_quality = 20
    a.cigar = ((0,10), (2,1), (0,25))
    a.next_reference_id = 0
    a.next_reference_start=199
    a.template_length=167
    a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")
    a.tags = (("NM", 1),
              ("RG", "L1"))
    outf.write(a)

同理,我们也可以读取一个已有的bam文件,逐个修改以上的属性,然后存储到一个新的bam文件里.这里不再举例.

上面设置header可能有点麻烦,容易出错,但我们可以复制一个已有bam文件的header到一个新的bam文件里.

outf = pysam.AlignmentFile(path_out, "wb", template=samfile)

以上template参数指定了模板bam文件.

 

5. 关闭文件

outf.close()

 

总结:

pysam模块非常实用,有了pysam模块,我们就可以非常灵活的操纵bam/bcf文件,而不必依赖于samtools或bcftools. pysam可以随机读取bam/bcf文件,也可以将处理后的内容自定义输出到bam/bcf文件.

 

以上只介绍了pysam最常见的功能,更多pysam功能请参照:http://pysam.readthedocs.io/en/latest/index.html