标签归档:Genome

思考:是否升级参考基因组版本

Should I switch to a newer reference?

GRCh38 consists of several components: chromosomal assembly, unlocalized contigs (chromosome known but location unknown), unplaced contigs (chromosome unknown) and ALT contigs (long clustered variations). The combination of the first three components is called the primary assembly. It is recommended to use the complete primary assembly for all analyses.
参考基因组包括chromosomal assembly, unlocalized contigs (chromosome known but location unknown), unplaced contigs (chromosome unknown)和ALT contigs (long clustered variations),前三个是primary assembly,alt contigs代表的是部分区域单体型的多样性,这些区域过于复杂不能用一条序列表示。

In addition to adding many alternate contigs, GRCh38 corrects thousands of SNPs and indels in the GRCh37 assembly that are absent in the population and are likely sequencing artifacts. It also includes synthetic centromeric sequence and updates non-nuclear genomic sequence.
除了添加了很多alternative contigs,GRCh38更正了数以千计的GRCh37版本中的SNPs和indels,这些SNPs和indels在人群中没有出现过、可能因为测序错误导致。GRCh38版本也包含了人工的着丝粒序列和更新了非核基因组序列。

The ability to recognize alternate haplotypes for loci is a drastic improvement that GRCh38 makes possible. Going forward, expanding genomics data will help identify variants for alternate haplotypes, improve existing and add additional alternate haplotypes and give us a better accounting of alternate haplotypes within populations. We are already seeing improvements and additions in the patch releases to reference genomes, e.g. the seven minor releases of GRCh38 available at the time of this writing.
GRCh38大幅提高了识别alternate haplotype的能力,进一步提高了识别alternate haplotype的突变的能力。

BWA的作者Li Heng推荐GRCh37 primary assembly+ALT+decoy组成的参考基因组hs38DH,可以通过bwa下载,见https://github.com/lh3/bwa/blob/master/README-alt.md

bwa.kit/run-gen-ref hs38DH

作者的用NA12878做测试的比对结果如下,

Assembly	hs37	hs38	hs38DH
FP	255706	168068	142516
TP	2142260	2163113	2150844

PS:
从改进来看,这些改进都将提高比对和突变检测的质量,
从结果看,38版本提高了TP的灵敏度,并降低了FP,
从趋势来看,37升级成38是大趋势,so请及时升级。

参考
https://software.broadinstitute.org/gatk/documentation/article.php?id=7857
https://gatkforums.broadinstitute.org/gatk/discussion/8180/9-takeaways-to-help-you-get-started-with-grch38https://software.broadinstitute.org/gatk/documentation/article.php?id=8017
https://github.com/lh3/bwa/blob/master/README-alt.md

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#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/

继续阅读

hg19、GRCH37、b37、hs37d5介绍和区别

大家经常用UCSC的hg19和NCBI的GRCh37版本的,但还有其他的版本,比如b37,hg37d5,比如在分析NIST的genome in a bottle(GIAB)提供的bam数据时,就遇到了hg37d5的版本,在用GATK的时候会遇到b37版本。

GRCh37

Genome Reference Consortium(基因组参照序列联盟),由英国Wellcome Trust Sanger研究中心(the Wellcome Trust Sanger Center)、华盛顿大学基因组中心(The Washington University Genome Center)、欧洲生物信息研究所(the European Bioinformatics Institute)和美国国家生物技术信息中心(NCBI)联合组成。

GRCH37版本发布之后,也会有小的更新,比如GRCh37.p2,大的更新比如由GRCh37升级到GRCh38,填补gap,修改部分序列,其目的是提供一个完整的基因组序列assemble。GRCh38已经在2013年发布,多数基因组数据库正在兼容或者更新到该版本。

该版本包含人类chr1到chr22,chrX,chrY,MT染色体以及

  • “unlocalized sequences”:知道来自哪条染色体但不知道具体位置的序列
  • “unplaced sequences”:知道来自人类基因组序列,但不知道与染色体的关系
  • “alternate loci”:来自基因组特定区域,代表该区域序列的多样性

“1” to “22”, “X”, “Y” and “MT”命名比较规范,ENSEMBL, genome browser, the NCBI dbSNP (in VCF files), the Sanger COSMIC (in VCF files),都依照该规范。

下载地址:ftp://ftp.ncbi.nih.gov/genomes/Homo_sapiens

GRCh37lite

只包含GRCh37版本中的chr1到chr22,chrX,chrY,MT染色体

下载地址:ftp://ftp.ncbi.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Homo_sapiens/GRCh37/special_requests/

hg19

UCSC提供,容易下载,因为UCSC方便下载各种坐标文件(bed,gtf等),该版本可以与这些坐标对应。与GRCh38对应的是hg38版本。

该版本序列包括chr1到chr22,chrX,chrY序列与GRCh37完全一致(完全一致,完全一致),线粒体序列稍微不一样,以及

  • “chr*_random sequences” 知道来自哪条染色体但不知道具体位置的序列
  • “chrUn_* sequences” 知道来自人类基因组序列,但不知道与染色体的关系
  • UCSC与GRCh不同的地方有:

  • 在重复区域repeat region有小写来表示,这点和GRCh不同
  • 此外染色体有chr前缀,而GRCh没有chr前缀。
  • 线粒体序列版本不一样
  • 下载地址:ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes

    b37

    来自1000 Genome第一阶段,被称为b37(GATK和IGV社区小组中经常使用),包含了GRCh37, the rCRS mitochondrial sequence(MT序列)
    unlocalized sequences和unplaced sequences以他们的检索号命名,比如”GL000191.1″, “GL000194.1”, etc.
    但是不(不不不)包含 alternate loci

    下载地址:ftp://ftp.broadinstitute.org/pub/seq/references

    hs37d5

    可以理解成b37的升级版,在1000 Genome第二阶段使用。该版本包含了b37数据,以及

  • human herpesvirus 4 type 1 sequence 人类疱疹病毒序列(“NC_007605”).
  • “decoy” sequence诱饵序列(名为hs37d5) 来自HuRef、BAC或者质粒克隆和NA12878,可以提高序列比对的正确率。
  • 此外在Y染色体上的性染色体同源区域PAR标为N碱基,这样对应的X染色体的区域当作二倍体
  • 这些更改有利于序列比对和突变检测,降低假阳性与b37兼容。

    PS:正如以前讲到的(文章《思考在比对时,关于是否将chr*_random和chrUn_*序列放在参考基因组中的思考》),序列比对不要只放chr1到chr22,chrX,chrY和MT(我见过国内知名测序公司有这么干的),其他序列和诱饵序列非常重要,可以提高比对的准确率,降低假阳性。

    详情见:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.slides.pdf
    下载地址:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence

    参考

    http://googlegenomics.readthedocs.io/en/latest/use_cases/discover_public_data/reference_genomes.html
    https://wiki.dnanexus.com/Scientific-Notes/human-genome

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

    思考–在比对时,关于是否将chr*_random和chrUn_*序列放在参考基因组中的思考

    通常认为chr1-22,chrY,chrX和chrM为参考基因组序列,于是包括我在内的很多人,分别下载了25条染色体序列,合并成一个fasta文件,用bowtie2或者BWA构建index,用于下一步的read比对,然后是各种分析(包括突变、转录表达等)。

    UCSC下载的HG19版本的整个参考基因组文件http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz中,除还包括chr*_random和chrUn_*序列(暂时理解为补丁序列,真实的补丁序列称呼常见assemble过程,见http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/info/patches.shtml,有fix 和novel patch,这里我们现在只讨论chr*_random和chrUn_*)。

    The chr*_random sequences are unplaced sequence on those reference
    chromosomes.

    The chrUn_* sequences are unlocalized sequences where the corresponding
    reference chromosome has not been determined.

    当然如果DNA或RNA测序的read比对到chr*_random 和 chrUn_* 序列上,显示大多数人都不关注这些序列上信息,我想这也是很多人不把chr*_random 和 chrUn_* 放到参考基因组fasta文件中的原因。但是,chr*_random 和 chrUn_* 显然是存在的,只不过现在暂时没有确定位置或者后续用于基因组序列更正。

    如果参考基因组序列中不包含chr*_random 和 chrUn_*序列,那么原来属于chr*_random 和 chrUn_*的read则有可能比对到(不是一定)chr1-22,chrX,chrY上的相似区域(这些区域与chr*_random 和 chrUn_*中的部分区域相似),造成假阳性比对,后续这些reads提供的信息都是不可靠的。

    如果参考基因组序列中包含chr*_random 和 chrUn_*序列,那么来自这些区域的reads则会正确的比对到这个地方,没有假阳性比对,只不过后续分析不需要考虑chr*_random 和 chrUn_*即可。

    假设有一条read来自chr1_random,
    条件                                     比对结果             分析结果
    基因组fa文件包含chr1_random序列          比对到chr1_random    后续不考虑
    基因组fa文件中不包含chr1_random序列      比对到chr1           造成假阳性

    举个例子,以前,我们确认一个突变是否存在,看覆盖这个点的read上有多少突变的碱基,如果覆盖这个点的read本来属于chr*_random 和 chrUn_*序列的,但比对到这个地方,即使这个位点的突变碱基再多,也是个假阳性突变,影响后续分析。

    综上,参考基因组需要放chr*_random和chrUn_*序列,降低reads比对时的假阳性。

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