如何對(duì)比vcf文件

如果我們要比較的兩個(gè)vcf文件的參考基因組版本不一致，就需要使用CrossMap等軟件進(jìn)行參考基因組版本轉(zhuǎn)換，然后里使用  SnpSift 軟件的 Concordance 命令比較它們。其中CrossMap軟件依賴pyBigWig，使用conda進(jìn)行安裝，代碼如下：

conda create -n py3 python=3.6
conda activate py3
conda install -c bioconda pyBigWig
pip3 install CrossMap

進(jìn)行參考基因組版本轉(zhuǎn)換的命令如下：

# 需要自行下載 hg19ToHg38.over.chain.gz  文件，以及參考基因組 Homo_sapiens_assembly38.fasta 
python ~/miniconda3/envs/py3/bin/CrossMap.py \
vcf ~/data/liftover/hg19ToHg38.over.chain.gz test.snp.hg19.vcf \
~/data/Homo_sapiens_assembly38.fasta test.snp.hg38.vcf

可以把snp和indel的vcf文件都轉(zhuǎn)換一下，然后拿到的轉(zhuǎn)換好的文件如下：

1.3M Jul  8 05:16 test.indel.hg38.vcf
 23K Jul  8 05:16 test.indel.hg38.vcf.unmap
1003K Jun 19 11:10 test.indel.vcf
 13M Jul  8 05:18 test.snp.hg38.vcf
245K Jul  8 05:18 test.snp.hg38.vcf.unmap
 13M Jun 19 18:29 test.snp.vcf

可以看到轉(zhuǎn)換的成功率是非常高的！unmap的文件很小，因?yàn)榇_實(shí)參考基因組有變化，總有一下基因組片段被修改了。

但是，有意思的是，之前我們的vcf文件是嚴(yán)格按照基因組坐標(biāo)排好序的，但是轉(zhuǎn)換過(guò)后，出現(xiàn)了部分坐標(biāo)亂序情況，如下：

這個(gè)很容易理解，因?yàn)橥粋€(gè)物種的不同版本參考基因組肯定是有

chr1    119955031       .       G       A
chr1    148483282       rs7513869       C       T
chr1    144995248       rs6600697       A       G
chr1    144995236       rs6600696       A       C
chr1    144995050       rs1884147       C       T
chr1    144995033       rs1884146       A       G

也就是說(shuō)，人類的參考基因組在由hg19進(jìn)化到hg38的時(shí)候，不僅僅是片段的自然擴(kuò)充，還包括一些以前組裝順序弄錯(cuò)了的片段的糾正。

這樣坐標(biāo)亂序的vcf文件，在很多下游分析都是不友好的，所以可以使用下面的代碼進(jìn)行簡(jiǎn)單過(guò)濾。

input=test.snps.VQSR.vcf
cat $input | java -jar ~/biosoft/snpEff/SnpSift.jar  filter "( DP > 20 & FILTER = 'PASS' )" | \
perl -alne '{print unless $F[0]  =~  /_/}' | \
awk '$1 ~ /^#/ {print $0;next} {print $0 | "sort -k1,1 -k2,2n"}' | \
grep -v '1/2' >  test.filter.sort.vcf
# 檢查不同染色體分布情況： 
cat new.filter.sort.vcf |grep -v '^#' |cut -f 1 |sort |uniq

# 接下來(lái)就可以對(duì)干凈的VCF文件進(jìn)行注釋啦 
java -jar ~/biosoft/snpEff/snpEff.jar GRCh48.86 \
test.filter.sort.vcf  > test.filter.sort.eff.vcf