1. raw vcf filter
mkdir vcf_filter && cd vcf_filter
mv ../raw.combine.vcf.gz ./
1.1 vcf snp過(guò)濾
#提取snp為單獨(dú)文件绽乔，
# -V 輸入文件
gatk SelectVariants -select-type SNP -V raw.combine.vcf.gz -O raw.snp.vcf.gz
#設(shè)置過(guò)濾參數(shù)進(jìn)行過(guò)濾
# -filter-expression "QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"`過(guò)濾參數(shù)
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
標(biāo)記過(guò)濾參數(shù)下要被過(guò)濾及合格的序列
#過(guò)濾
#`--exclude-filtered true`過(guò)濾掉非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
#查看過(guò)濾前和過(guò)濾后文件的內(nèi)容，第七列標(biāo)記了pass或是要被過(guò)濾的信息
# cut 文本切割符碳褒，默認(rèn)為制表符折砸，-f指定顯示第7列
#uniq -c 顯示每行連續(xù)出現(xiàn)的次數(shù)
過(guò)濾前后序列數(shù)量比較
1.2 vcf indel過(guò)濾
步驟與snp相同看疗，過(guò)濾參數(shù)不同
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
#過(guò)濾后的snp與indel合并（可以不合并）
gatk MergeVcfs -I filtered.snp.vcf.gz -I filtered.indel.vcf.gz -O combine.filtered.vcf.gz
2. ANNOVAR注釋
2.1 輸入文件格式轉(zhuǎn)換
# convert2annovar.pl 將多種格式轉(zhuǎn)為ANNOVAR使用.avinput格式;
# -format vcf4 表明輸入文件`filtered.indel.vcf.gz`是vcf格式的內(nèi)容
# -allsample 將輸入filtered.indel.vcf.gz中3個(gè)樣本按每個(gè)樣本輸出到單一的.aviput文件
dir=/home/qmcui/software/ANNOVAR/annovar/
$dir/convert2annovar.pl -format vcf4 -allsample filtered.indel.vcf.gz -outfile anno
2.2 ANNOVAR注釋功能
#根據(jù)基因注釋
cat >id.txt
anno.7E5239
anno.7E5240
anno.7E5241
# -geneanno 基于基因的注釋
# -buildver hg38 默認(rèn)使用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
基于基因注釋*為信息查看文件
分析參考代碼給學(xué)徒的WES數(shù)據(jù)分析流程
分析視頻全外顯子測(cè)序數(shù)據(jù)分析
