搜索
查看: 2614|回复: 0

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-2 09:40:51 | 显示全部楼层 |阅读模式
【直播】我的基因组(十二):先粗略看看几个基因吧

昨天我们说到,测序得到的fastq文件map到基因组之后,我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format,而BAM就是SAM的二进制文件。通常sam文件太大,我们会生成bam文件来节省空间。sam文件和bam文件的转换用samtools这个软件就可以完成。


[Perl] 纯文本查看 复制代码
samtools view -h abc.bam > abc.sam
samtools view -b -S abc.sam > abc.bam

我们已经拿到了bam文件,我这里就先用公司给我的bam文件吧,根据我的帖子:[url=]仅仅对感兴趣的基因call variation[/url] ,可以先了解几个比较有趣的基因的变异情况。我自己呢,对以下几个位点和基因比较感兴趣,就用他们来讲一下今天的内容吧!




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基因内或周围。

那,我们就先关注这几个基因吧(不要问我为什么(-_-メ) )。
首先找到这些基因的坐标,看到如下:

b8cb935a7afbebe75aa053eb4e9af2aa.png

其中V1aR基因这个雄性标志性出轨基因,在标准的基因命名系统里面其实是AVPR1A:[url=]http://www.genecards.org/cgi-bin/carddisp.pl?gene=AVPR1A[/url] ,这里面涉及到HUGO symbol的概念,这个genecard数据库也非常赞,基因相关信息都可以在这里面查找的。

有了这些坐标信息,我们就进入我们的基因组工作目录:

cd data/project/myGenome/

然后把坐标文件做好


2105c245cb8251a07b4e82719d10b5af.png

因为公司给我的bam文件里面,用的参考基因组是GRCh37而不是hg19(两者区别在于chr是否标记),我们还是需要下载;

[Perl] 纯文本查看 复制代码
cd ~/reference
mkdir -p  genome/human_g1k_v37  && cd genome/human_g1k_v37
# [url=http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/]http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/[/url]
nohup wget [url=http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz]http://ftp.1000genomes.ebi.ac.uk ... an_g1k_v37.fasta.gz[/url]  &
gunzip human_g1k_v37.fasta.gz
wget [url=http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai]http://ftp.1000genomes.ebi.ac.uk ... n_g1k_v37.fasta.fai[/url]
wget [url=http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/README.human_g1k_v37.fasta.txt]http://ftp.1000genomes.ebi.ac.uk ... n_g1k_v37.fasta.txt[/url]

然后回到基因组工作目录,保证bam文件在上图中bamFiles那个目录,然后用下面这个脚本,批量提取我们感兴趣的基因的变异情况:


[Perl] 纯文本查看 复制代码
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

等三分钟就好了,结果如下:

03055eaf87a89126b1046f094e3860e6.png

前面我们说到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因,那么我们很容易去dbSNP数据库或者我最近强烈推荐 的snpedia数据库([url=]吐血推荐snpedia数据库,非常丰富的snp信息记录[/url])里面找到它的坐标。

6        32666295 :Rs9275319--HLA-DQ
2        191964633 :Rs7574865--STAT4
然后我检查了我刚才call到的variation文件,

zcat STAT4.vcf.gz |grep -w 191964633 显示为空。

zcat HLA-DQ* |grep 32666295  也是空。

哈哈,我完美的错过了这两个易感位点!!!!谢天谢地!!!
其余的我就不讲了,毕竟会涉及到隐私,我就讲这个方法吧!

文:Jimmy、吃瓜群众

图文编辑:吃瓜群众







上一篇:【直播】我的基因组(十一):测序数据的比对
下一篇:【直播】我的基因组(十三):了解sam格式比对结果
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-1 12:07 , Processed in 0.047078 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.