我对H3F3A这个基因做了两个突变的cellline,分别是G34V和K27M,现在知道这个基因在hg38上面的坐标是:
Genomic Location for H3F3A Gene
Chromosome: 1
Start:226,061,851 bp from pter End:226,072,002 bp from pter
Size:10,152 bases Orientation:Plus strand
然后我用samtools结合bcftools把该基因区域的snp位点call出来:
samtools mpileup -r chr1:226061851-226072001 -t "DP4" -ugf ~/reference/genome/hg38/hg38.fa *sorted.bam | bcftools call -vmO z -o H3F3A.vcf.gz
但是得到的vcf只有DP4和染色体起始终止坐标坐标信息,我并不知道该坐标是蛋白质的第几个位点,所以需要注释,我首先想到的就是ANNOVAR啦,毕竟用了它很久。
~/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old H3F3A.vcf >tmp.annovar
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg38 --geneanno --outfile tmp.anno tmp.annovar ~/biosoft/ANNOVAR/annovar/humandb/
但是注释过后,很诡异的事情发生了!只有一个位点被认为是exon什么的,而且造成的蛋白质改变是G35R,很明显不是我所设计的突变位点,我设计的是G34V,它们这么近,我怀疑还是基因坐标表现形式的问题,而且该位点测序深度高达6000,应该是没有问题 的
line4 nonsynonymous SNV H3F3A:NM_002107:exon2:c.G103A:p.G35R, chr1 226064454 226064454 G A hom 219 6592 60
然后我查看了那些不在exon区域的位点,发现了更奇怪的事情,居然全部在H3F3AP4上面,这个时候我就傻眼了,这个假基因命名定位在
/home/jianmingzeng/reference/gtf/gencode/allGene.hg19.position:chr2 175584636 175585046 H3F3AP4
/home/jianmingzeng/reference/gtf/gencode/allGene.hg38.position:chr2 174719908 174720318 H3F3AP4
怎么也不可能跑到chr1来呀!!!!ANNOVAR到底是如何给我注释的!!!!
我只好去查ANNOVAR的database,发现它居然真的有如此无厘头的记录:
grep H3F3AP4 humandb/hg38_refGene.txt
2309 NR_002315 chr1 + 226062726 226072002 226072002 226072002 4 226062726,226064328,226065655,226071350, 226062811,226064479,226065809,226072002, 0 H3F3AP4 unk unk -1,-1,-1,-1,
1918 NR_002315 chr2 + 174719799 174720841 174720841 174720841 1 174719799, 174720841, 0 H3F3AP4 unk unk -1,
一个基因被记录两个位置,让我好生郁闷!!!而且H3F3AP4很明显是与H3F3A重合了的,我敢打包票,肯定是某人写脚本的时候,没有考虑周全,跟我上一个文章提到的原因一模一样,搞这些数据库维护的单位太多了,总会有不一致的地方。
2309 NM_002107 chr1 + 226062706 226072002 226064351 226071479 4 226062706,226064328,226065655,226071350, 226062811,226064479,226065809,226072002, 0 H3F3A cmpl cmpl -1,0,2,0,
所以,当我们尤其是想确认某一个问题的事情,请务必再三检查!!!