1. raw vcf filter
mkdir vcf_filter && cd vcf_filter
mv ../raw.combine.vcf.gz ./
1.1 vcf snp过滤
#提取snp为单独文件，
# -V 输入文件
gatk SelectVariants -select-type SNP -V raw.combine.vcf.gz -O raw.snp.vcf.gz
#设置过滤参数进行过滤
# -filter-expression "QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"`过滤参数
gatk VariantFiltration -V raw.snp.vcf.gz --filter-expression "QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O filter.snp.vcf.gz
标记过滤参数下要被过滤及合格的序列
#过滤
#`--exclude-filtered true`过滤掉非pass的
gatk SelectVariants --exclude-filtered true -V filter.snp.vcf.gz -O filtered.snp.vcf.gz
zless -S filter.indel.vcf.gz|grep -v '^#'|cut -f 7 |sort|uniq -c
zless -S filtered.snp.vcf.gz|grep -v '^#'|cut -f 7 |sort|uniq -c
#查看过滤前和过滤后文件的内容，第七列标记了pass或是要被过滤的信息
# cut 文本切割符，默认为制表符，-f指定显示第7列
#uniq -c 显示每行连续出现的次数
过滤前后序列数量比较
1.2 vcf indel过滤
步骤与snp相同，过滤参数不同
gatk SelectVariants -select-type INDEL -V raw.combine.vcf.gz -O raw.indel.vcf.gz
gatk VariantFiltration -V raw.indel.vcf.gz --filter-expression "QUAL < 30.0 | QD < 2.0 || FS > 200.0 || SOR > 10.0 || Inbreeding Coeff < -0.8 || ReadPosRankSum < -20.0"
--filter-name "Filter" -O filter.indel.vcf.gz
gatk SelectVariants --exclude-filtered true -V filter.indel.vcf.gz \
-O filtered.indel.vcf.gz
#过滤后的snp与indel合并（可以不合并）
gatk MergeVcfs -I filtered.snp.vcf.gz -I filtered.indel.vcf.gz -O combine.filtered.vcf.gz
2. ANNOVAR注释
2.1 输入文件格式转换
# convert2annovar.pl 将多种格式转为ANNOVAR使用.avinput格式;
# -format vcf4 表明输入文件`filtered.indel.vcf.gz`是vcf格式的内容
# -allsample 将输入filtered.indel.vcf.gz中3个样本按每个样本输出到单一的.aviput文件
dir=/home/qmcui/software/ANNOVAR/annovar/
$dir/convert2annovar.pl -format vcf4 -allsample filtered.indel.vcf.gz -outfile anno
2.2 ANNOVAR注释功能
#根据基因注释
cat >id.txt
anno.7E5239
anno.7E5240
anno.7E5241
# -geneanno 基于基因的注释
# -buildver hg38 默认使用hg18
#`$id.avinput`输入文件
dir=/home/qmcui/software/ANNOVAR/annovar/
db=$dir/humandb/
ls $db
cat id.txt|while read id; do ($dir/annotate_variation.pl -geneanno -buildver hg38 $id.avinput $db);done
基于基因注释*为信息查看文件
分析参考代码给学徒的WES数据分析流程
分析视频全外显子测序数据分析
