vcftools是一種可以對VCF文件和BCF文件進行格式轉(zhuǎn)換及過濾的工具,功能非常強大胡陪,而且運算速度也很快沥寥。
1.下載及安裝
1.1 下載地址
http://vcftools.sourceforge.net/downloads.html
1.2 安裝
進入壓縮包目錄碍舍,進行解壓。
$ tar xvf vcftools_0.1.13.tar.gz
$ cd ~/vcftools_0.1.13/bin
關(guān)于安裝营曼,有一些小伙伴留言遇到了問題乒验,我再詳細寫一下安裝的問題:
$ cd ~/vcftools_0.1.1
$ ./configure
$ make
$ make install
檢查安裝是否成功
$ vcftools
VCFtools (v0.1.13)
? Adam Auton and Anthony Marcketta 2009
Process Variant Call Format files
For a list of options, please go to:
https://vcftools.github.io/examples.html
Questions, comments, and suggestions should be emailed to:
vcftools-help@lists.sourceforge.net
但是vcftools安裝確實容易出現(xiàn)各種報錯,所以建議用conda來安裝:
$ conda install -c bioconda vcftools
2. 基礎(chǔ)用法
2.1 vcf文件加ID
add id.pl拷貝至vcftools/bin目錄下
add id.pl是一個老師寫的腳本蒂阱,這里不好直接放上來锻全,所以需要添加id的話,請大家再去查找其他的教程录煤,這里只是我自己做個備份鳄厌。
perl add_id.pl root.hic.vcf root.hic.id.vcf
2.2 分開indel和SNP
只輸出indel
vcftools --vcf root.hic.id.vcf --keep-only-indels --recode --recode-INFO-all --out root.hic.id.indel
只保留SNP
vcftools --vcf root.hic.id.vcf --remove-indels --recode --recode-INFO-all --out root.hic.id.snp
2.3 vcf文件過濾
vcftools --vcf root.hic.id.vcf --max-missing 0.8 --maf 0.05 --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out root.hic.id.int0.8maf0.05.allele2
max-missing:分型完整度
maf:第二等位基因頻率
min-alleles:最小等位基因個數(shù)
max-alleles:最大等位基因個數(shù)
哈迪溫伯格平衡 hwe filtering(通常在人類中使用)
vcftools --vcf root.hic.id.vcf --remove-indels --max-missing 0.8 --maf 0.05 --min-alleles 2 --max-alleles 2 --hwe 0.01 --recode --recode-INFO-all --out root.hic.id.snp.hwe0.01
2.4 文件格式轉(zhuǎn)換
轉(zhuǎn)換vcf格式為ped map格式
vcftools --vcf root.hic.id.vcf --plink --out root.hic
2.5 vcf文件拆分
準備一個samplelist文件,即需要的樣本的ID
vcftools --vcf hic.sort.ref.vcf --recode --recode-INFO-all --keep samplelist.txt --out root
詳細說明書見官網(wǎng):
http://vcftools.sourceforge.net/
2.6 統(tǒng)計等位基因頻率
第一條染色體上的等位基因頻率妈踊,注意vcf文件中chr的大小寫了嚎。
結(jié)果中第一列是染色體;第二列SNP位置廊营;第三列是這個位置一共出現(xiàn)了幾種堿基歪泳,這里是兩個;第四列是該位置出現(xiàn)的堿基總數(shù)露筒,這里一個樣本貢獻了兩個堿基位點呐伞;后面是該位置出現(xiàn)的堿基對應(yīng)的頻率。
$ vcftools --vcf root.id.vcf --freq --chr Chr1 --out Chr1_analysis
$ head -n 5 Chr1_analysis.frq
CHROM POS N_ALLELES N_CHR {ALLELE:FREQ}
Chr1 5154 2 286 C:0.713287 T:0.286713
Chr1 5187 2 296 A:0.614865 G:0.385135
Chr1 5220 2 282 A:0.241135 T:0.758865
Chr1 50889 2 294 A:0.870748 G:0.129252
計算整體文件中的等位基因頻率
$ vcftools --vcf root.id.vcf --freq --out allel_analysis
$ tail -n 5 allel_analysis.frq
scaffold995 3028 2 304 G:0.927632 A:0.0723684
scaffold995 3082 2 304 G:0.842105 A:0.157895
scaffold995 3168 2 298 T:0.946309 C:0.0536913
scaffold995 3185 2 292 C:0.89726 T:0.10274
scaffold995 3228 2 302 A:0.407285 C:0.592715
2.7 比較兩個vcf文件的變異位點
結(jié)果文件中的第一列是染色體慎式;第二列和第三列是文件1和文件2中SNP位置伶氢;第四列中B表示兩個文件中都有這個堿基,如果是1則表示只有文件1中有這個堿基瘪吏,如果是2同理癣防。
$ vcftools --vcf root.id.vcf --diff fei.id.vcf --diff-site --out in1_v_in2
$ head -n 5 in1_v_in2.diff.sites_in_files
CHROM POS1 POS2 IN_FILE REF1 REF2 ALT1 ALT2
Chr1 5154 5154 B C C T T
Chr1 5187 5187 B A A G G
Chr1 5220 5220 B A A T T
Chr1 50889 50889 B A A G G