标签归档:liftover

liftover,crossmap进行坐标转换时用到的chain文件介绍

chain description
在做基因组坐标转换的时候,用crossmap和liftover的时候,会用到chain file。大家都在讲坐标转换会用到这个文件,却没有讲过chain文件的具体内容(百度和谷歌搜索结果都没有中文介绍,本文应该是第一个)。本文的内容都翻译自UCSC网站,原文https://genome.ucsc.edu/goldenpath/help/chain.html,希望能帮到大家了解这个文件。

UCSC chain文件 http://hgdownload.soe.ucsc.edu/goldenPath/hg38/liftOver/
Ensembl chain文件 https://sourceforge.net/projects/crossmap/files/Ensembl_chain_files/

chain file里面包含许多块alignment的信息(个人觉得可以理解为同源的地方,chain),其中每一块有一个header,记录alignment在两个版本中坐标,以及许多行alignment data line记录具体比对情况。

即每一块有 Header Line 和 Alignment Data Lines组成。形式如下,有两个chain

    chain 4900 chrY 58368225 + 25985403 25985638 chr5 151006098 - 43257292 43257528 1
     9       1       0
     10      0       5
     61      4       0
     16      0       4
     42      3       0
     16      0       8
     14      1       0
     3       7       0
     48

     chain 4900 chrY 58368225 + 25985406 25985566 chr5 151006098 - 43549808 43549970 2
     16      0       2
     60      4       0
     10      0       4
     70

Header Lines

记录两个版本的序列对应区域

chain score tName tSize tStrand tStart tEnd qName qSize qStrand qStart qEnd id

header line以chain开头
score — chain 值
tName — 染色体 (reference sequence)
tSize — 染色体大小 (reference sequence)
tStrand — 链strand (reference sequence)
tStart — 比对开始 (reference sequence)
tEnd — 比对结束 (reference sequence)
qName — 检索染色体 (query sequence,对应的另外一个基因组版本 )
qSize — 检索染色体大小 (query sequence)
qStrand — 链strand (query sequence)
qStart — 比对开始 (query sequence)
qEnd — 比对结束 (query sequence)
id — chain ID

位置从0开始,前100bp的表示方法为0到100,接下的100bp为100到200,strand的值为-号时,对应的是反向互补序列。

Alignment Data Lines

Alignment Data Line有三列,记录区域内的具体对应信息,最后一行只有ungapped的序列长度。

    size dt dq

size — ungapped(一致的区域)的长度
dt — 该block块(一致的区域)的结束位置,到下一个block快(下一个一致的区域)起始位置的长度,其实就是gapped block不一致的区域的长度 (reference sequence)
dq — 该block块(一致的区域)的结束位置,到下一个block快(下一个一致的区域)起始位置的长度,其实就是gapped block不一致的区域的长度(query sequence)

以chain id为2的chain为例

     chain 4900 chrY 58368225 + 25985406 25985566 chr5 151006098 - 43549808 43549970 2
     16      0       2
     60      4       0
     10      0       4
     70

我个人理解的,不一定正确,欢迎指正

chain description

坐标转换在reference和query的位置为:
reference: start 25985406 end 25985566
query: start 43549808 end 43549970

区域长度分别为:
reference region length = 25985566 – 25985406 = 160 bp
query region length = 43549970 – 43549808 = 162 bp

一致的区域(ungapped)长度为
reference = query = 16 + 60 +10 + 70 = 156 bp

不一致的区域(gapped)长度为
reference = 0 + 4 +0 =4
query = 2 + 0 + 4 =6

区域长度也可以通过ungapped和gapped区域计算,分别为:
reference = ungapped + gapped = 156 + 4 =160 bp
query = ungapped + gapped = 156 + 6 =162 bp

更多信息详见:

http://genomewiki.ucsc.edu/index.php/LiftOver_Howto

https://genome.ucsc.edu/goldenpath/help/chain.html

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

基因组坐标转换工具-以BED文件为例,从hg19转换到hg38坐标

分析时使用的基因组版本,可能会与其他来源数据所使用的基因组版本不一致,需要统一成同一个版本的坐标,才能方便下一步的分析。

常用的有NCBI的Remap在线服务和UCSC的liftover,其实还有很多,本文暂时总结部分工具的用法。以将APOA1的编码区坐标(利用UCSC的genome browser下载,或者下载该文件APOA1.bed)转换为例,从hg19转到hg38版本坐标上。需要注意的是,在使用的时候,需要注意是否支持对应的格式。

类型 支持格式 地址 推荐指数
Liftover 在线 bed http://genome.ucsc.edu/cgi-bin/hgLiftOver 一般
Liftover 本地 bed和gff http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver 推荐
Remap 在线 hgvs,bed,gvf,gff,gtf,Text ASN.1,Binary ASN.1,UCSC Region和VCF https://www.ncbi.nlm.nih.gov/genome/tools/remap 推荐
CrossMap 本地 SAM/BAM,,Wiggle/BigWig, bed, gff/gtf,VCF http://crossmap.sourceforge.net/ 推荐
picard 本地 interval和VCF http://broadinstitute.github.io/picard/ 推荐VCF转换

在线版

Liftover

hg-liftover-ucsc
服务地址: http://genome.ucsc.edu/cgi-bin/hgLiftOver

    • 1,首先选择检索和目标基因组
Original Genome:	Human
Original Assembly:	Feb. 2009 (GRCh37/hg19)
New Genome:  	Human
New Assembly:	Dec. 2013 (GRCh38/hg38)   
  • 2,将bed格式内容复制或者上传,如果是复制内容的话,点击submit;如果是上传文件的 话,点击submit
  • 3,在Result部分下载结果,Result部分会显示结果说明,点击View Conversions,即能得到转换后的bed文件。

Remap

NCBI-remap
支持hgvs,bed,gvf,gff,gtf,Text ASN.1,Binary ASN.1,UCSC Region和VCF,如果数据量较少,可以考虑该方法。

服务地址:https://www.ncbi.nlm.nih.gov/genome/tools/remap

    • 1,在Assembly-Assembly中选择
Source Organism *:Homo sapiens
Source Assembly *:   GRCh37(hg19)
Target Assembly *:    GRCh38(hg38)
  • 2,在Data处,可以选择复制文件内容还是上传文件,还可以指定文件格式。
  • 3,点击submit,页面会跳转到结果页面,Summary Data告诉你有多少条匹配,Mapping Report告诉你对应关系。这些都可以下载,NCBI remap的结果文件前几列是原来的内容,后几列是转换后的坐标。

Ensembl Assembly Converter

http://asia.ensembl.org/Homo_sapiens/Tools/AssemblyConverter?db=core

本地版

本地工具需要利用chain file,才能知道两个版本的坐标对应关系。chain文件可以从UCSC下载。我们需要的是hg19Thg38.over.chain.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz
chain文件后续会有介绍。

CrossMap

CrossMap支持SAM/BAM, Wiggle/BigWig, BED, GFF/GTF, VCF格式。支持的格式较多,用起来方便。

软件网址:http://crossmap.sourceforge.net/

安装:pip install CrossMap

当然也可以从源码编译安装,官网有说明(http://crossmap.sourceforge.net/#install-crossmap-from-source-code)。

python CrossMap.py bed    hg19ToHg38.over.chain      APOA1.bed      APOA1.hg38.bed

Liftover

支持bed和gff格式

wget -c http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver
chmod 755 liftOver
# liftOver oldFile map.chain newFile unMapped
liftOver   APOA1.bed    hg19ToHg38.over.chain    APOA1.hg38.bed   unMapped.txt

Picard

地址: http://broadinstitute.github.io/picard/

继续阅读