文件格式 | 对VCF文件的常见操作¶
发表于: 2024-02-20
压缩¶
bcftools view sativa332.indel.Lsat_1_v5_gn_7_15020.hard.vcf -Oz -o sativa332.indel.Lsat_1_v5_gn_7_15020.hard.vcf.gz
-Oz: 指定输出文件为压缩格式 gz.
-o: 指定输出文件名。
建立索引¶
只能对压缩后的vcf文件建立索引
bcftools index -t sativa332.indel.Lsat_1_v5_gn_7_15020.hard.vcf.gz
-t: 填压缩后的 vcf.gz 文件.
tabix -p vcf view.vcf.gz
-p: 填压缩后的 vcf.gz 文件.
上述命令任选其一, 均可完成建索引.
提取变异¶
gatk SelectVariants -R Lactuca_sativa_Salinas_V8.fa -V Lactuca.snp.vcf.gz -L chr7:18846094-18857694 -L chr8:1884094-188009694 -O sativa332.snp.Lsat_1_v5_gn_7_15020.hard.vcf.gz
-R: 参考基因组fasta文件
-V: 输入vcf文件
-L: 需要提取的变异区间,可以为 chrx:xx-xx 或直接提取整条染色体 chrx, 但格式必须与输入vcf文件相同,如: chr1, Chr1, Chr01, 可以同时提取多个区间 (重复使用该参数)
-O: 输出vcf文件名
bcftools view -S keep.list test.vcf > sub_indv.vcf
-S: 需要提取的变异的位置信息,可以是 染色体\t具体位置 两列，也可以是 染色体\t起始\t终止 三列.
上述命令任选其一即可
提取样品¶
gatk SelectVariants -R Lactuca_sativa_Salinas_V8.fa -V Lactuca.snp.vcf.gz -sn s001 -sn s002 -sn s003 -O sativa332.snp.Lsat_1_v5_gn_7_15020.hard.vcf
-R: 参考基因组fasta文件
-V: 输入vcf文件
-O: 输出vcf文件名
-sn: 需要提取的样品名,可以同时提取多个样品 (重复使用该参数)
bcftools view view.vcf.gz -s NA00001,NA00002 -o subset.vcf
-s: 需要提取的样品名,,分割
-o: 输出文件名
bcftools view input.SNP.revhet.recode.vcf.gz -S input_samples.list -Oz -o output_fixed.vcf.gz
-S: 需要提取的样品名称列表
vcftools --gzvcf in.vcf.gz --recode --recode-INFO-all --stdout --keep id.txt > out.vcf.gz
--gzvcf: 输入文件名
--recode: 告诉程序使用过滤器写入一个新的vcf文件
--recode-INFO-all: 保留旧vcf文件中的所有INFO标志
--stdout: 配合后面的管道函数，如果输出文件是压缩文件 --stdout > out.vcf.gz (也可以不压缩),不压缩时，可以直接 --out out 此处第一个out是命令，第二个out是文件名 (也可以 --stdout > out.vcf)
--keep: 所需要的样品名列表
选择一个即可完成
去除样品¶
适用于去除少数样品
bcftools view example.vcf.gz -s ^NA00001,NA00002 -o subset.vcf
-s ^ 后接需要去除的样品名, , 分割
-o 输出文件名
适用于较多样品
vcftools --gzvcf in.vcf.gz --recode --recode-INFO-all --stdout --remove id.txt > out.vcf.gz
--gzvcf: 输入文件名
--recode: 告诉程序使用过滤器写入一个新的vcf文件
--recode-INFO-all: 保留旧vcf文件中的所有INFO标志
--stdout: 配合后面的管道函数，如果输出文件是压缩文件 --stdout > out.vcf.gz (也可以不压缩), 不压缩时，可以直接 --out out 此处第一个 out 是命令，第二个 out 是文件名 (也可以 --stout > out.vcf)
--remove: 所需要的样品名列表
修改样品名¶
bcftools的常用命令
bcftools reheader --samples sub.list -o output.vcf.gz input.vcf.gz
--samples: 内容为oldname\tnewname, 可有多行 (修改多个样品名)
MAF 过滤¶
参考: VCFtools的使用(参数说明)
/zfssz5/BC_PUB/Software/03.Soft_ALL/vcftools-0.1.13/bin/vcftools --gzvcf sativa332.hard.Lsat_1_v5_gn_7_15020.snp_indel.vcf.gz --max-missing 0.9 --maf 0.05 --recode --recode-INFO-all --stdout | /ldfssz1/ST_AGRIC/P18Z10200N0148_LETTUCE/USER/bin/01.software/smcpp1.13.1/bin/bgzip -c > sativa332.maf.Lsat_1_v5_gn_7_15020.snp_indel.vcf.gz
--gzvcf: 输入文件名
--max-missing: 缺失率，0为接受完全缺失，1为接受数据全都存在
--maf: Minor Allele Frequency二等位基因频率进行过滤，通常保留大于0.05的。
--recode : 告诉程序使用过滤器写入一个新的vcf文件
--recode-INFO-all: 保留旧vcf文件中的所有INFO标志
--stdout: 配合后面的管道函数，如果输出文件是压缩文件 --stdout > out.vcf.gz (也可以不压缩), 不压缩时，可以直接 --out out 此处第一个 out 是命令，第二个 out 是文件名 (也可以 --stout > out.vcf)
bgzip -c >: 压缩文件
提取染色体名称¶
bioawk -c vcf ' { print $chrom } ' old.vcf.gz | sort -n | uniq > ChrName.txt
sort -n :按数值排序
uniq: 去重
更改染色体名称¶
bcftools annotate --rename-chrs NewChrName.txt old.xxx.vcf.gz -Oz -o new.xxx.gz
rename-chr: 接新老名称的对应关系,其中NewChrName.txt文件中第一列为旧ID，第二列为对应的新ID
-Oz: 指定输出文件为压缩格式 gz
-o: 指定输出文件名
按染色体拆分vcf¶
for i in {chr1..chr9};do echo -e "gatk SelectVariants -R Lactuca_sativa_Salinas_V8.fa -V Lactuca.snp.vcf.gz -L chr7:18846094-18857694 -L $i -O sativa332.snp.$i.hard.vcf.gz";done
循环使用提取变异命令
按染色体合并vcf (推荐转换成bed文件合并，更快速)¶
参考: 多款软件进行vcf合并-gatk、vcftools、bcftools
gatk GatherVcfs -I concat-a.vcf -I concat-b.vcf -O combine_a_b_samesample_diffsites.vcf
-I: 输入文件名(可重复使用该参数)
-O: 输出文件名
gatk MergeVcfs -I concat-a.vcf -I concat-b.vcf -O combine_a_b_diffsample_allsites_gatk.vcf
-I: 输入文件名(可重复使用该参数)
-O: 输出文件名
bcftools concat concat-a.vcf concat-b.vcf -o combine_a_b_samesample_diffsites_bcftools.vcf
concat:后接需要合并的vcf名
-O: 输出文件名
合并SNP和indel¶
bcftools concat -R Lsat_1_v5_gn_7_15020.region -a sativa332.indel.hard.vcf.gz sativa332.snp.vcf.gz -Oz -o sativa332.hard.snp_indel.vcf.gz
-R: 只有包含在这个文件中指定的区域的记录才会被合并,格式为 chr1\t1389\t111333,可有多行(多个区间)
-a: 需要合并的文件,示例中的两个文件将按顺序合并
-Oz: 指定输出文件为压缩格式 gz.
-o: 指定输出文件名。
vcf的基本信息¶
bcftools stats example.vcf.gz > example.vcf.stats
提取样品ID信息¶
bcftools query -l example.vcf.gz >example.id.list
-l: 后接需要查询的vcf文件
提取变异信息¶
参考: bcftools常用命令详解
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' view.vcf.gz
%CHROM: 代表VCF文件中染色体那一列，其他的列，比如POS, ID, REF, ALT, QUAL, FILTER也是类似的写法
[]: 对于FORMAT字段的信息，必须要中括号括起来
%SAMPLE: 代表样本名称
%GT: 代表FORMAT字段中genotype的信息