昨天我们说到，测序得到的fastq文件map到基因组之后，我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format，而BAM就是SAM的二进制文件。通常sam文件太大，我们会生成bam文件来节省空间。sam文件和bam文件的转换用samtools这个软件就可以完成。
samtools view -h abc.bam > abc.sam samtools view -b -S abc.sam > abc.bam
我们已经拿到了bam文件，我这里就先用公司给我的bam文件吧，根据我的帖子：仅仅对感兴趣的基因call variation ，可以先了解几个比较有趣的基因的变异情况。我自己呢，对以下几个位点和基因比较感兴趣，就用他们来讲一下今天的内容吧！
1.STAT4上的rs7574865和HLA-DQ的 rs9275319是中国人群中乙型肝炎病毒（HBV）相关肝细胞癌（HCC）遗传易感基因
2.V1aR基因是雄性标志性出轨基因。
3.GLI3和PAX1基因控制鼻孔的大小，而RUNX2基因控制鼻梁的宽度。DCHS2基因调控鼻子的突起程度，即决定鼻尖是否朝上和鼻尖的角度，或者说它决定了你的鼻子是否迷人挺拔。
4.肥胖有关的基因FTO（Fat Mass and Obesity Associated），最近发现了调控肥胖（主要是脂肪燃烧）的基因是IRX3 和IRX5。大约100个基因位点与BMI(身体质量指数)相关，600个基因位点与身高相关，160个基因位点与肥胖特征如腰臀比相关。6个新基因位点，这些位点位于LEMD2、CD47、GANAB、RPS6KA5/C14orf159、ANP32和ARL15基因内或周围。
那，我们就先关注这几个基因吧（不要问我为什么(-_-メ) ）。
首先找到这些基因的坐标，看到如下：
其中V1aR基因这个雄性标志性出轨基因，在标准的基因命名系统里面其实是AVPR1A：http://www.genecards.org/cgi-bin/carddisp.pl?gene=AVPR1A ，这里面涉及到HUGO symbol的概念，这个genecard数据库也非常赞，基因相关信息都可以在这里面查找的。
有了这些坐标信息，我们就进入我们的基因组工作目录：
cd data/project/myGenome/
然后把坐标文件做好
因为公司给我的bam文件里面，用的参考基因组是GRCh37而不是hg19(两者区别在于chr是否标记)，我们还是需要下载；
cd ~/reference
mkdir -p genome/human_g1k_v37 && cd genome/human_g1k_v37
# http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/
nohup wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz &
gunzip human_g1k_v37.fasta.gz
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/README.human_g1k_v37.fasta.txt
然后回到基因组工作目录，保证bam文件在上图中bamFiles那个目录，然后用下面这个脚本，批量提取我们感兴趣的基因的变异情况：
cat key_gene.list |while read id;
do
chr=$(echo $id |cut -d" " -f 1|sed 's/chr//' )
start=$(echo $id |cut -d" " -f 2 )
end=$(echo $id |cut -d" " -f 3 )
gene=$(echo $id |cut -d" " -f 4 )
echo $chr:$start-$end $gene
samtools mpileup -r $chr:$start-$end -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta bamFiles/P_jmzeng.final.bam | \
bcftools call -vmO z -o $gene.vcf.gz
done
等三分钟就好了，结果如下：
前面我们说到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒（HBV）相关肝细胞癌（HCC）遗传易感基因，那么我们很容易去dbSNP数据库或者我最近强烈推荐 的snpedia数据库（吐血推荐snpedia数据库，非常丰富的snp信息记录）里面找到它的坐标。
6 32666295 :Rs9275319--HLA-DQ
2 191964633 :Rs7574865--STAT4
然后我检查了我刚才call到的variation文件，
zcat STAT4.vcf.gz |grep -w 191964633 显示为空。
zcat HLA-DQ* |grep 32666295 也是空。
哈哈，我完美的错过了这两个易感位点！！！！谢天谢地！！！
其余的我就不讲了，毕竟会涉及到隐私，我就讲这个方法吧！
文：Jimmy、吃瓜群众
图文编辑：吃瓜群众