十二 09

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

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

图文编辑:吃瓜群众

十一 23

【直播】我的基因组(七):从整体理解全基因组测序数据的变异位点

首先记住一个很重要的知识点,变异是相对的!

简单说一下什么是找变异,变异跟突变有什么区别呢?举个栗子:有国际组织规定了人类的参考基因组(如UCSC,ENSEMBL,NCBI等,前面帖子都有讲),就是 AAAAA(这里简化一下,就5个碱基,其实人类基因组多达30亿个)  。现在通过给自己测序得知,我与之对应的是AGCAA,那么我相比国际基因组来说,就是2个变异位点,位于基因组的坐标2和3,但是它们还不能说就是突变。 Continue reading

10

根据dbSNP的ID来转换成HGVS突变表示形式

dbSNP的ID直接在NCBI的dbSNP官网可以看到详细介绍,现在已经更新到146版本了,一般人看到一个ID肯定什么信息都获取不到,毕竟这只是人家NCBI规定的一个ID而已。但是HGVS突变形式就有非常详细的信息了。

人类基因组变异协会(HGVS)官方组织规定了mutation该如何记录:http://www.hgvs.org/mutnomen/recs.html  推荐大家都仔细阅读!!!

Continue reading

14

华盛顿大学把所有的变异数据都用自己的方法注释了一遍,然后提供下载

华盛顿大学把所有的变异数据都用自己的方法注释了一遍,然后提供下载:
文献是:Kircher M, Witten DM, Jain P, O'Roak BJ, Cooper GM, Shendure J. 

A general framework for estimating the relative pathogenicity of human genetic variants.
Nat Genet. 2014 Feb 2. doi: 10.1038/ng.2892.
PubMed PMID: 24487276.

文中的观点是:现在大多的变异数据注释方法都非常单一,通常是看看该位点是否保守,对蛋白功能的改变,在什么domain上面等等。
但这样是远远不够的,所以他们提出了一个新的注释方法,用他们自己的CADD方法把现存的一些公共数据库的变异位点(约86亿的位点)都注释了一下,并对每个位点进行了打分。
C scores correlate with allelic diversity, annotations of functionality, pathogenicity, disease severity, experimentally measured regulatory effects and complex trait associations, and they highly rank known pathogenic variants within individual genomes.
总之,他们的方法是无与伦比的!
所有他们已经注释好的数据下载地址是:http://cadd.gs.washington.edu/download
这些数据在很多时候非常有用,尤其是想跟自己得到的突变数据做交叉验证,或者做一下统计分析的时候!
clipboard
人的基因组才300亿个位点,他们就注释了86亿!!!
所以有三百多G的压缩包数据,我想,一般的公司或者单位都不会去用这个数据了!