文章
标题：Genotyping structural variants in pangenome graphs using the vg toolkit
中文：使用vg对泛基因组结构变异进行基因分型
杂志：Genome Biol.
时间：2020
被引：66（谷歌学术2021.12.17）
Github: https://github.com/vgteam/vg
Bioconda: https://anaconda.org/bioconda/vg
conda安装
conda create -n SV
conda activate SV
conda install -c bioconda vg
# version v1.36.0
conda安装中没找到测试数据，因此重新下载Linux同版vg v1.36.0
地址：https://github.com/vgteam/vg/releases/download/v1.36.0/vg-v1.36.0.tar.gz
主要piplines
main mapping and calling pipeline:
-- autoindex mapping tool-oriented index construction from interchange formats
-- construct graph construction
-- rna construct splicing graphs and pantranscriptomes
-- index index graphs or alignments for random access or mapping
-- map MEM-based read alignment
-- giraffe fast haplotype-aware short read alignment
-- mpmap splice-aware multipath alignment of short reads
-- augment augment a graph from an alignment
-- pack convert alignments to a compact coverage index
-- call call or genotype VCF variants
-- help show all subcommands
1 Variation graph construction
# 工作地址
#/hwfssz1/ST_HEALTH/P18Z10200N0423/hutongyuan/softwares/vg-v1.36.0/test
cd /hwfssz1/ST_HEALTH/P18Z10200N0423/hutongyuan/softwares/vg-v1.36.0/test/
mkdir small_out
vg construct \
-r small/x.fa \
-v small/x.vcf.gz \
> small_out/x.vg
# x.vg是一个二进制文件
参数
x.fa 基因组
x.vcf.gz 基因组VCF
-r, --reference FILE
-v, --vcf FILE
2 Viewing, conversion
# GFA output
vg view small_out/x.vg > small_out/x.gfa
# dot output suitable for graphviz
vg view -d small_out/x.vg > small_out/x.dot
# And if you have a GAM file
cp small/x-s1337-n1.gam small_out/x.gam
# json version of binary alignments
vg view -a small_out/x.gam > small_out/x.json
3 Mapping
# construct the graph (paths below assume running from `vg/test` directory)
vg construct -r small/x.fa -v small/x.vcf.gz > small_out/x.vg
# store the graph in the xg/gcsa index pair
vg index \
-x small_out/x.xg \
-g small_out/x.gcsa \
-k 16 small_out/x.vg
# align a read to the indexed version of the graph
# note that the graph file is not opened, but x.vg.index is assumed
vg map -s CTACTGACAGCAGAAGTTTGCTGTGAAGATTAAATTAGGTGATGCTTG \
-x small_out/x.xg \
-g small_out/x.gcsa \
> small_out/read.gam
# simulate a bunch of 150bp reads from the graph, one per line
vg sim \
-n 1000 -l 150 \
-x small_out/x.xg \
> small_out/x.sim.txt
# now map these reads against the graph to get a GAM
vg map \
-T small_out/x.sim.txt \
-x small_out/x.xg \
-g small_out/x.gcsa \
> small_out/aln.gam
# surject the alignments back into the reference space of sequence "x", yielding a BAM file
vg surject \
-x small_out/x.xg \
-b small_out/aln.gam \
> small_out/aln.bam
# or alternatively, surject them to BAM in the call to map
vg sim -n 1000 -l 150 \
-x small_out/x.xg \
> small_out/x.sim.txt
vg map \
-T small_out/x.sim.txt \
-x small_out/x.xg \
-g small_out/x.gcsa --surject-to bam \
> small_out/aln.bam
4 Augmentation
vg augment \
small_out/x.vg \
small_out/aln.gam \
-A small_out/aug.gam \
> small_out/aug.vg
vg augment \
small_out/x.vg \
small_out/aln.gam -i -S \
> small_out/aug_with_paths.vg
5 Variant Calling
Calling variants using read support
# Compute the read support from the gam
# -Q 5: ignore mapping and base qualitiy < 5
# -s 5: ignore first and last 5bp from each read
vg pack \
-x small_out/x.xg \
-g small_out/aln.gam \
-Q 5 -s 5 \
-o small_out/aln.pack
# Generate a VCF from the support.
vg call small_out/x.xg \
-k small_out/aln.pack \
> small_out/graph_calls.vcf
############################################
vg call \
small_out/x.xg \
-k small_out/aln.pack -a \
> small_out/snarl_genotypes.vcf
############################################
# Index our augmented graph
vg index \
small_out/aug.vg \
-x small_out/aug.xg
# Compute the read support from the augmented gam (ignoring qualitiy < 5, and 1st and last 5bp of each read)
vg pack \
-x small_out/aug.xg \
-g small_out/aug.gam \
-Q 5 -s 5 \
-o small_out/aln_aug.pack
# Generate a VCF from the support
vg call small_out/aug.xg \
-k small_out/aln_aug.pack \
> small_out/calls.vcf
#########################################
# Re-construct the same graph as before but with `-a`
vg construct \
-r small/x.fa \
-v small/x.vcf.gz -a \
> small_out/xa.vg
# Index the graph with `-L' to preserve alt paths in the xg
vg index small_out/xa.vg -x small_out/xa.xg -L
# Compute the support (we could also reuse aln.pack from above)
vg pack \
-x small_out/xa.xg \
-g small_out/aln.gam \
-o small_out/aln.pack
# Genotype the VCF (use -v)
vg \
call small_out/xa.xg \
-k small_out/aln.pack \
-v small/x.vcf.gz \
> small_out/genotypes.vcf
##################################################
# filter secondary and ambiguous read mappings out of the gam
vg filter \
small_out/aln.gam \
-r 0.90 -fu -m 1 -q 15 -D 999 \
-x small_out/x.xg \
> small_out/aln.filtered.gam
# then compute the support from aln.filtered.gam instead of aln.gam in above etc.
vg pack \
-x small_out/xa.xg \
-g small_out/aln.filtered.gam \
-o small_out/aln.pack
vg call \
small_out/xa.xg \
-k small_out/aln.pack \
-v small/x.vcf.gz \
> small_out/genotypes.vcf
vg snarls small_out/x.xg \
> small_out/x.snarls
# load snarls from a file instead of computing on the fly
vg call \
small_out/x.xg \
-k small_out/aln.pack \
-r small_out/x.snarls \
> small_out/calls.vcf
Calling variants from paths in the graph
mkdir small_hsp
# create a graph from a multiple alignment of HLA haplotypes (from vg/test directory)
vg msga \
-f GRCh38_alts/FASTA/HLA/V-352962.fa \
-t 1 -k 16 | vg mod -U 10 - | vg mod -c - \
> small_hsp/hla.vg
# index it
vg index \
small_hsp/hla.vg \
-x small_hsp/hla.xg
# generate a VCF using gi|568815592:29791752-29792749 as the reference contig. The other paths will be considered as haploid samples
vg deconstruct \
small_hsp/hla.xg \
-e -p "gi|568815592:29791752-29792749" \
> small_hsp/hla_variants.vcf
网友评论