【直播】我的基因组(十二):先粗略看看几个基因吧

昨天我们说到,测序得到的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、吃瓜群众

图文编辑:吃瓜群众

Comments are closed.