家系WES的vcf文件探索

应该是走的GATK流程，得到每个三口之家的每个单独个体的gvcf然后再合并成一个vcf文件，如下：

355M Jul 25  2016 H13_H14_H15.snp.vcf
441M Jul 25  2016 H16_H17_H18.snp.vcf
350M Jul 25  2016 H1_H2_H3.snp.vcf
297M Jul 25  2016 H22_H23_H24.snp.vcf
293M Jul 25  2016 H25_H26_H27.snp.vcf
304M Jul 25  2016 H29_H30_H31.snp.vcf
308M Jul 25  2016 H33_H34_H35.snp.vcf
304M Jul 25  2016 H39_H40-H41.snp.vcf
310M Jul 25  2016 H42_H43_H44.snp.vcf
317M Jul 25  2016 H45_H46_H47.snp.vcf
288M Jul 25  2016 H53_H54_H55.snp.vcf
445M Jul 25  2016 H5_H6_H7.snp.vcf
294M Jul 25  2016 H60_H61_H62.snp.vcf
309M Jul 25  2016 H63_H64_H65.snp.vcf
324M Jul 25  2016 H66_H67_H68.snp.vcf
401M Jul 25  2016 H70_H71_H72.snp.vcf
409M Jul 25  2016 H74_H75_H76.snp.vcf
390M Jul 25  2016 H9_H10_H11.snp.vcf

根据性染色体的杂合纯合SNP数量及比例判断性别

首先，使用脚本进行统计性染色体的基因型数量：

ls /home/yyfu/snp.vcf/*.vcf|while read id;
do
echo $id
cat $id|grep -v "^#" |grep -w "PASS" |grep -w chrX|cut -f 10|cut -d":" -f 1|sort |uniq -c
cat $id|grep -v "^#" |grep -w "PASS" |grep -w chrX|cut -f 11|cut -d":" -f 1|sort |uniq -c
cat $id|grep -v "^#" |grep -w "PASS" |grep -w chrX|cut -f 12|cut -d":" -f 1|sort |uniq -c
cat $id|grep -v "^#" |grep -w "PASS" |grep -w chrY|cut -f 10|cut -d":" -f 1|sort |uniq -c
cat $id|grep -v "^#" |grep -w "PASS" |grep -w chrY|cut -f 11|cut -d":" -f 1|sort |uniq -c
cat $id|grep -v "^#" |grep -w "PASS" |grep -w chrY|cut -f 12|cut -d":" -f 1|sort |uniq -c
done 

然后对统计后的结果进行汇总：

grep  -v '2/' log2.txt |grep -v '/2'|grep vcf|cut -d"/" -f 5|sed 's/.snp.vcf//g' |tr '\n' '\t' |perl -alne '{print "$_\t$_" foreach @F}'|sed 's/[_-]/\t/g' |tr '\n' '\t'
grep  -v '2/' log2.txt |grep -v '/2'|grep -v "vcf"|perl -alne '{print if $.%4==1}'|awk '{print $1}'
grep  -v '2/' log2.txt |grep -v '/2'|grep -v "vcf"|perl -alne '{print if $.%4==1}'|awk '{print $1}'>1
grep  -v '2/' log2.txt |grep -v '/2'|grep -v "vcf"|perl -alne '{print if $.%4==2}'|awk '{print $1}'>2
grep  -v '2/' log2.txt |grep -v '/2'|grep -v "vcf"|perl -alne '{print if $.%4==3}'|awk '{print $1}'>3
grep  -v '2/' log2.txt |grep -v '/2'|grep -v "vcf"|perl -alne '{print if $.%4==0}'|awk '{print $1}'>4
paste 1 2 3 4

还可以提取这些VCF的性染色体信息拿去做生物信息学工程师练习题目：

ls /home/yyfu/snp.vcf/*.vcf|while read id;
do
echo $id
file=$(basename $id)
sample=${file%%.*} 
echo calling $sample 
cat $id|grep -v "^#" |grep -w "PASS" |grep -w "chr[XY]" >$sample.sex.snp.vcf
done 

统计全部的基因型比例

因为性染色体上面的遗传不准确，可以略去

ls /home/yyfu/snp.vcf/*.vcf|while read id;
do
echo $id
file=$(basename $id)
sample=${file%%.*} 
echo calling $sample 
cat $id|grep -v "^#" |grep -w "PASS"|grep -v "chr[XY]" |cut -f 10-12|grep -v  "\.\/\."|grep  -v '2/'|grep -v '/2'|perl -ane '{print substr($_,0,3)."\t" foreach @F;print "\n"}' |sort |uniq -c |sort -k1,1nr|awk '{print $1"\t"$2"\t"$3"\t"$4}' >$sample.stat
done 

这样得到的文件，比如 H13_H14_H15.stat 就可以判断哪些符合孟德尔遗传定律，哪些不符合，结合上面的性别信息：

H13是女性，H14是男性，H15是女性！！！ 因为男性只有一条X染色体，所以其X染色体的突变都应该是纯合的！

chromosome	sample	no_call	wt	het	hom	het/hom
chrX	H13	1709	551	1293	4432	0.291741877
chrX	H14	3197	1034	233	3525	0.066099291
chrX	H15	1875	597	1261	4253	0.296496591
chrY	H13	665	194	663	1553	0.426915647
chrY	H14	494	345	939	1295	0.725096525
chrY	H15	665	192	674	1541	0.437378326

可以得到 左边是女儿，中间是父亲，右边是母亲

54378   1/1 1/1 1/1
18852   0/1 0/1 0/1
14724   0/1 0/1 0/0
14333   0/1 0/0 0/1
12700   0/0 0/0 0/1
12345   0/0 0/1 0/0
7816    1/1 1/1 0/1
7653    1/1 0/1 1/1
6610    0/1 1/1 0/0
6567    0/1 0/0 1/1
6512    0/0 0/1 0/1
6499    0/1 0/1 1/1
6162    1/1 0/1 0/1
6155    0/1 1/1 0/1
4830    1/1 1/1 0/0
4807    1/1 0/0 1/1
1856    0/0 1/1 1/1
1778    0/1 1/1 1/1
1456    0/0 0/0 1/1
1423    0/0 1/1 0/0
1353    1/1 0/0 0/1
1297    1/1 0/1 0/0
927 0/0 0/1 1/1
873 0/0 1/1 0/1
809 0/1 0/0 0/0
565 1/1 0/0 0/0