12.19 vcftools找两组snp的不同位点

要查找两个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的赋值，采用" "