由于OMIM上的位置是参考基因组GRch38,所以在进行hg19版本的annovar注释时,需要转化为hg19的,根据OMIM数据库上的提示,我们可以从gff文件中获取对应的位置。(说一下个人想法:我们看到要进行位置转化,首先想到是用UCSC上的liftover进行转化,或者是NCBI上的remap进行转化,但是据说这个在2023年11月就要停止了。如果对整个数据的查看,就会发现位置不是唯一性的,里面很多同一个位置对应好几个MIM number或者好几个Entrez Gene ID,可能是这个位置太长了,里面对应的基因太多,所以如果我们简单的进行位置转化,在注释的时候,根据位置匹配,就会出现多个基因,最好的是根据OMIM数据库给的提示操作,把对应的位置范围缩小,才能更精准)
从gff文件根据MIM number或者Entrez Gene ID来获取位置,尝试使用R,python(biopython、bcbio-gff)。R是因为安装包容易报错,bcbio-gff可以解析gff文件,但是因为gff文件太大,python在处理时速度慢,均放弃。也尝试过在线网站bioDBnet转化,结果为GRch38,也放弃。最后用shell处理。
cut -f1 id.txt |sed 1d | sed '/^$/d'| while read i;do zgrep $i GRCh37_latest_genomic.gff.gz;done | perl -alne '{next unless $F[2] eq "gene"; ;/ID=gene:(.*?);/; print "$F[0]\t$F[3]\t$F[4]\t$F[8]"}'
id.txt文件的列名:MIM number Gene Symbols Entrez Gene ID,先提取第一列(MIM number ),然后去掉第一行,再去掉空行,之后对每行MIM number循环读取,在GFF文件中寻找(而且只要GFF中gene的信息),并将相关信息匹配提取出来,因为还需要核对,所以输出内容是整行完整的。