要查找兩個vcf文件的不同位點井仰，使用vcftools的 --diff語句峰尝。
但是遇到一個問題儒旬，有兩個文件不同時存在的chr温兼，需要找出秸滴。
# 取第一列
awk '{print $1}' Ping72_filtered.vcf > ping72.txt
# 取所有不以#開頭的行
grep -v "^#" ping48.txt > ping48-1.txt
# 排序，刪除重復(fù)的
sort ping48-1.txt | uniq > ping48-2.txt
#find duplicate （將兩個文件合并募判，只保留只出現(xiàn)過一次的荡含，也就是差異的）
cat ping72-2.txt ping48-2.txt >ping72-48.txt
sort ping72-48.txt |uniq -u >ping72-48u.txt
本來想直接將txt文件作--not-chr 后的參數(shù)，但是vcftools不認識届垫，只能把這個語句擴充為一個shell文件了释液。
# 將文件的換行全部變成相應(yīng)的字符
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
接著之前的繼續(xù)做。
我們還需要找到差異的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
# 這個語句用來判斷文件夾是否已經(jīng)存在寝蹈，不判斷會報錯
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的賦值，采用" "
