flanking sequence for candidate

1. 得到了cadidate SNPs list

Scaffold861: 64225
Scaffold44:1265163
Scaffold121: 842037
Scaffold225: 196363
Scaffold266: 463359
Scaffold68: 899972
...

2. 获取vcf文件

file=338cand.snps
sed -i 's/\t/:/g' $file
sed -i 's/ /:/g' $file
sed -i 's/:::/:/g' $file
sed -i 's/::/:/g' $file
sed -i 's/:/\t/g' $file
zgrep -F -f  $file gatk4.SNP_VF_SV.vcf.gz > 338cands_no_header.vcf
bcftools view -h gatk4.SNP_VF_SV.vcf.gz > header
cat header 338cands_no_header.vcf > 338cands.vcf

bcftools view -H 338cands.vcf|wc -l
grep -v "#" 338cands.vcf|grep Scaffold|wc -l

检查数目是否一致

SNP注释

DIR=/home/mmcui/poolseq/gatk_hard_filtering/howto_filter
VCF=338cands
java -Xmx4g -jar  ~/annotation/snpEff/snpEff.jar -v Agla_Btl03082013 $DIR/$VCF.vcf > $VCF.ann.vcf

Scaffold1:1235332-1237332
Scaffold1:4378505-4380505
Scaffold1:4078202-4080202
Scaffold1:1639666-1641666
Scaffold1:1555820-1557820

retrieve flanking region

while read line
do 
samtools faidx 210Scaffolds.fasta $line >>flanking_sequences.snps
done < 337snps.positions