生信log40|只有GRCh37或GRCh38的人基因组数据怎么
2023-10-05 本文已影响0人
小周的万用胶囊
最近在查看一些数据集，发现有些非常老的数据库里面只有GRCh37的数据，一些经常更新的数据库里面只有GRCh38的数据，如果恰好手头上有两个版本的基因组数据可就非常头疼了。为了找到数据，我掘地三尺，把网络数据都翻了一下，发现有些官方网站如VEP会提供两个版本的数据，并且上面标注了"remapped with CrossMap"，说明这些数据是转换而来的。
VEP数据库
1、CrossMap的实现原理
- 根据文件中记载的序列之间的位点关系进行基因坐标的转换（但软件的不进行相似性的比较，而是从外部直接获取）
2、CrossMap支持什么文件
CrossMap支持什么文件
目前官方支持BED、GFF/GTF、VCF文件，这些文件里面都包含了染色体、染色体的起始位点信息。
哪些信息是这些文件里面共有的？
这些文件里面都包含了染色体、染色体的起始位点信息
3、CrossMap的用法
准备文件： 此处以转换VCF为例子
- 待转换的VCF文件
- 转换链文件： .chain或.lz格式
- 转换版本的参考基因组：如要把数据转换成GRCh38的数据，那么
PS：转换链文件是通过UCSC的LiftOver工具得到的，官方有提供可直接下载。
- 1、软件安装
# 法一
conda install -c bioconda crossmap
#法二
pip3 install CrossMap
- 2、转换数据
基础命令：CrossMap.py vcf input.chain input.vcf refgenome.fa out_vcf
以转换ExAC数据为例
- 由于资源有限，此处仅抽取一万行的数据做实验。
CrossMap.py vcf \
references_grch37_to_grch38.over.chain.gz \
ExAC.hg19.test.vcf \ #you
references_Homo_sapiens_assembly38.fasta \
ExAC.hg38.test.out.nofilter.vcf --no-comp-allele
- 3、数据验证
因为已经有现成的GRCh38版本的ExAC数据了，因此直接拿转换好的数据进行对比，由于转换时间和版本的问题，后续的一些VCF补充信息可能与官方的不一致，此处仅截取头五列五行（#CHROM,POS,ID,REF,ALT）的数据作对比。
本地转换好的数据
本地转换好的数据
网上的数据
网上的数据
有上两图可见这5X5的数据是一致的！
其他
- 下载参考基因组和链文件需要注意版本，否则会发现转换之后的文件没有对应位点的数据。
PS：未完待续 ...
参考
CrossMap文章
CrossMap官网
文中链文件和参考基因组下载地址
