要查找两个vcf文件的不同位点，使用vcftools的 --diff语句。
但是遇到一个问题，有两个文件不同时存在的chr，需要找出。
# 取第一列
awk '{print $1}' Ping72_filtered.vcf > ping72.txt
# 取所有不以#开头的行
grep -v "^#" ping48.txt > ping48-1.txt
# 排序，删除重复的
sort ping48-1.txt | uniq > ping48-2.txt
#find duplicate （将两个文件合并，只保留只出现过一次的，也就是差异的）
cat ping72-2.txt ping48-2.txt >ping72-48.txt
sort ping72-48.txt |uniq -u >ping72-48u.txt
本来想直接将txt文件作--not-chr 后的参数，但是vcftools不认识，只能把这个语句扩充为一个shell文件了。
# 将文件的换行全部变成相应的字符
awk '{print $1}' ping72-48u.txt |xargs |sed 's/ / --not-chr /g'>test72-48.sh
之后打开文件，插入缺失的语句（其实也可以写一个shell文件来运作，但是我只有几个需要操作的vcf，所以手动）
vcftools --vcf Ping72_filtered.vcf --diff Ping48_filtered.vcf --diff-site —-not-chr your file --out compare-72-48
2019-01-02
接着之前的继续做。
我们还需要找到差异的snp位点。
#找到含有这个基因注释的行
cat Ping_filtered.vcf| grep "BGIBMGA007793" >abc.vcf
#把需要的列提取出来
awk '{print $1, $2, $3, $4, $5}' 72.vcf > 72_abc.txt
如果我们批量化的话就写一个shell代码
for gene in $genes
do
geneid="BGIBMGA${gene}"
outdir="abc_${gene}"
# mkdir out
cd abc_gene_family
# 这个语句用来判断文件夹是否已经存在，不判断会报错
if [ ! -d $outdir ]
then mkdir $outdir
fi
cd ..
#find snps
ids="2 30 48 72"
for id in $ids
do
file="Ping${id}_filtered.vcf"
cat $file | grep $geneid > abc_gene_family/$outdir/${id}.vcf
awk '{print $1,$2,$3,$4,$5}' abc_gene_family/$outdir/${id}.vcf > abc_gene_family/$outdir/${id}.txt
done
这里我们可以看一下这个ids的赋值，采用" "
