人基因组分析环境搭建

生信喵

一、软件安装

1. #1 下载biconda
   
2. wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh  
   
3. #2 安装
   
4. sh Miniconda3-latest-Linux-x86_64.sh  
   
5. source ~/.bashrc
   
6. #3 添加软件源
   
7. conda config --add channels bioconda
   
8. conda config --add channels conda-forge
   
9. #4 安装mamba
   
10. conda install -y mamba
    
11. 
    
12. #创建虚拟环境
    
13. conda create -n human -y
    
14. conda activate human
    
15. conda install -y fastqc
    
16. conda install -y multiqc
    
17. conda install -y fastp
    
18. conda install -y gatk4
    
19. conda install -y bwa
    
20. conda install -y bwa-mem2
    
21. conda install -y samtools
    
22. conda install -y bcftools
    
23. conda install -y minimap2
    
24. conda install -y longshot
    
25. conda install -y cutesv 没装上
    
26. conda install -y ngmlr
    
27. conda install -y sniffles
    
28. conda install -y seqkit
    
29. conda install -y snpeff
    
30. conda install -y snpsift 没装上
    
31. conda install -y freebayes

复制代码

二、参考序列下载
       变异检测需要将测序数据与参考序列进行序列比对，比较二者基因组上的差异。因此在做变异检测的过程中，除了样品本身测序对分析结果有影响，参考序列同样有重要的影响。一些重要物种往往有多个版本的参考序列可供选择，因此选择合适的参考序列是非常重要的准备工作。
2.1 人全基因组版本介绍
       做人全基因的变异检测，首先需要下载人全基因组参考序列。参考序列主要是一些已发表物种基因组的一些信息，可以包括基因组序列，基因的核酸序列，基因的氨基酸序列，坐标文件，注释信息等，需要针对不同的分析目的，选择合适的数据进行下载。由于分析目的的不同，人全基因组参考序列包含很多的版本，同一版本有衍生出很多特殊版本，例如重复序列是否屏蔽，是否带端粒，中心粒，孢疹病毒，ALT 序列，rCRS mitochondrial，是否包含没有定位到基因组上的序列，以及参考序列 ID 是否与后续注释数据库相同等。需要考虑的因素比较多。
       常见的参考序列主要存储在国际三大核酸数据库 NCBI，EMBL，DDBJ 等，还包括 UCSC 以及一些物种单独的网站上。由于参考序列一般比较大，这里推荐使用 ftp 工具进行下载，ftp工具是专门的针对 ftp 文件传输协议的工具，下载速度更快，并且支持断点续传，可以使用xftp 或者 filezilla 等访问 ftp 进行下载，也可以直接通过命令行 ftp 工具进行下载，例如 lftp命令等。

1. NCBI: ftp://ftp.ncbi.nlm.nih.gov/
   
2. EMBL：ftp://ftp.ensembl.org/pub/
   
3. UCSC: http://genome.ucsc.edu/
   
4. JGI：https://jgi.doe.gov/
   
5. ENSEMBL：http://asia.ensembl.org/info/about/species.html
   

复制代码

2.2 不同版本人参考序列之间的关系
       众所周知人类基因组主要构成包括 22 条常染色体（1-22），2 条性染色体（X，Y）和线粒体 DNA（mtDNA），但是由于拼接结果不完整等原因，人全基因组参考序列存在多个版本。
       而且即使是同样的序列，存储在不同数据库的命名也有很大差别，例如来自 UCSC 数据库中，一般命名为 hg18,19，hg38。而在 NCBI 数据库中命名为 GRch36,37,38，GRC 是 Genome Reference Consortium 的简称，而 embl 数据库中则是 release 然后跟一个版本号，有非常多的版本，千人基因组命名为 b36，b37。
       目前，广泛使用的版本是GRCH37和GRCH38。2009年，the Genome Reference Consortium (GRC)发布了第 19 版人类基因组 GRCH37，也常被称为 hg19。GRCH37 被广泛应用于数据分析。2013年，GRC 发布了 GRCH38。但由于注释工具、数据库的不健全及升级基因组工作繁杂，时至今日，GRCH37 仍被相当程度地使用。
       不同数据库采用不同的命名方式，UCSC 命名最简单，就是 hg 加数字，比如 hg18，hg19 和hg38 这些版本比较常用；NCBI 的命名更复杂一些，是 GRCH 开头，然后数字，36,37,38 等。而 ENSEMBL 命名就更复杂了，只给出版本比如 52,59,61,80,81,82 等。这三种命名方式有一个对应关系。如下表所示：
       不同版本基因组比较

发布时间	UCSC	NCBI	ENSEMBL	1000G
-	hg18	GRCh36	release_52	b36
2009	hg19	GRCh37	release_59/61/64/68/69/75	b37
2013	hg38	GRCh38	76/77/78/80/81/82	---

2.3 GRCH37 与 GRCH38 版本详细差别
       根据 GRC 的官方文件，GRCH38 是最精确的人类基因组。GRCH38 基于金标准 Sanger 测序组装，读长约为 1000bp，精确度是高通量测序的 10 倍。与 GRCH37 相比，GRCH38 替换了 8000 个等位基因位点，校正了数个组装错误的基因组区域，补全了 gap，添加了着丝粒序列，在 178个区域组装了 261 条 alternate loci，丰富了基因组的多样性。
       已发表的论文认为 GRCH38 是 GRCH37 的重大升级，可提供更精确的生物信息学和基因组学分析。我们设计了实验量化基于 GRCH38 和 GRCH37 的数据分析差异。

1. 参考材料 1：https://doi.org/10.1016/j.ygeno.2017.01.005
   
2. 参考材料 2：
   
3. https://bitesizebio.com/38335/get-to-know-your-reference-genome-grch37-vs-grch38/

复制代码

GRCH37 与 GRCH38 详细比较
       不算线粒体 DNA，GRCH37 和 GRCH38 分别有 3095677412 和 3088269832 个核苷酸。最常用的线粒体基因组是 1999 年剑桥发布的 rCRS，因此两者线粒体基因组是一样的。在基因组 fasta文件中，’N’表示 gap 或者未注释区域，GRCH37 共有 234350281 个‘N’，而 GRCH38 中有150630719 个，减少了 83719562 个，占比 35.7%。从表 1 中看出，每条染色体上的‘N’数量都有减少。有文献研究表明 GC 含量影响 Illumina 测序深度及测序均一性，这与后续的CNV 检测密切相关。GC 位点的总数从 GRCH37 的 1170371008 增加到 GRCH38 的 1200551672，共计增加了 30180664 个核苷酸。
       外显子可以编码氨基酸，是人类基因组最重要的组成部分。从 Ensembl (GRCh37 v37.75, GRCh38 v38.82)下载最新的 Gene Feature Format (GTF)文件统计外显子区域。外显子区域由 GRCH37 的 75231228 个核苷酸增加到 GRCH38 的 95505476 个，约有 26.9%的增幅。从全基因组水平看，外显子占比由 2.43%增至 3.09%。外显子区域扩大的主要原因有 3 个：i.在GRCH38 中，外显子的总数从 327058 个增加到 457748 个；ii.每个基因的外显子数从 13 个增加到 19 个；iii.每个外显子核苷酸的中位数从 140 增加到 146。

2.4 选择合适的基因组版本
       那么到底该选择哪个版本的基因组呢？首先要清楚不同版本之间的详细差别。不同版本的参考基因组，主要有几点差别：
       第一、序列不同，早前测序的基因组版本与后面测序的基因组版本，例如 hg19 和 hg38 之间，序列会有所差别，导致基因组坐标也不同。
       第二、里面包含的内容不同，由于人的全基因组并不是非常完美的 24 条，会有很多序列指导染色体号但是无法定位到具体位置，例如 chr2_KI270716v1_random；还有一些无法定位到具体染色体上，例如 chrUn 开头的，还有一些包含 EBA，HLA 序列等。
       第三、同一个序列，在不同库之间命名会有所有不同，例如有的版本基因组命名包括 chr，例如 chr1，chr2，而有些则只包含染色体号。
       第四、重复序列屏蔽方式，对于重复区域，有些基因组中会将重复序列替换为 N，成为“hardmask”，或者替换为小写字母，称为“soft mask”，有些则不作任何处理。
       除此之外，还与具体的分析目的有关系，例如要做 SNP 检测，就需要屏蔽重复序列，如果做RNAseq 则不需要。这里有有一篇博客，详细阐述不同版本基因组选择的问题。

1. https://lh3.github.io/2017/11/13/which-human-reference-genome-to-use

复制代码

三、不同版本人全基因组序列下载
       下面具体来介绍一下各个版本之间的差别，以及如何下载具体的序列。
3.1 ncbi
       Genome Reference Consortium(基因组参照序列联盟)，由英国 Wellcome Trust Sanger 研究中心（the Wellcome Trust Sanger Center）、华盛顿大学基因组中心（The WashingtonUniversity Genome Center）、欧洲生物信息研究所（the European Bioinformatics Institute）和美国国家生物技术信息中心（NCBI）联合组成。
       该版本包含人类 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),都依照该规范。
       下载地址：

1. https://www.ncbi.nlm.nih.gov/projects/genome/guide/human/

复制代码

生信喵

1. #重复序列是否屏蔽，是否带端粒，中心粒，孢疹病毒，ALT序列， rCRS mitochondrial，没有定位到基因组上的序列，序列ID不同。
   
2. #1 NCBI
   
3. #NCBI hs37-1kg  与千人基因组计划一致
   
4. ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz
   
5. #NCBI hs37d5 deocy版本
   
6. ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
   
7. #NCBI hg38
   
8. ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

复制代码

3.2 ensembl
       Ensembl 来自于欧洲分子生物学中心。Ensembl 中主要是参考序列的基因组，按照动物，植物，真菌，原核生物的方式划分。ENSEMBL 的优势是具有较全的注释信息。包括 GTF 类型的文件格式。Ensembl 采用按照文件格式组织的组织方式。current 表示最新版本，然后最下面有很多 release，最下面也有一个 README 文件。因为 embl 定期会对不同物种参考序列注释信息进行完善，每完善一次就会发布出来，所以就有了很多 release，其实每个 release里中包含所有参考序列物种信息，现在最新的是 release95。
       如果要要下载人全基因组的 fasta 文件和 gtf 文件，分别进到对应目录查找。首先找到home_sapiens。目录中的 READM 给出详细信息。每个染色体单独给出，也包含都在一个文件内的。那么文件中包含 primary_assembly，toplevel，rm，sm，这些有什么差别呢。其实主要是因为重复序列的影响，因为有些时候做变异检测的需要 Mask 屏蔽掉重复区，处理重复序列有多种方式，可以直接替换为 N，这个称为 hard masker，也就是 rm，也可以将重复序列替换为小写字母，这个称为 soft masker，sm，primary 是原始没有处理过的。而 toplevel表示包含位点多态性信息，在做基因组变异检测时可以使用。

1. http://ensemblgenomes.org/

复制代码

1. #EnsEMBL primary
   
2. http://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
   
3. #EnsEMBL toplevel
   
4. http://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.toplevel.fa.gz

复制代码

3.3 ucsc
       UCSC 是加利福尼亚大学圣克鲁兹分校（University of California, Santa Cruz，简称 UCSC），该网站提供常见参考序列基因组的信息。除了提供常规基因组序列下载之外，还有一个基因组浏览器，可以查看基因组上的不同区域。常见的 hg18，hg19 以及 hg38 等命名方式就来自于 UCSC。找到下载地址，直接就可以进行下载。除此之外，UCSC 还可以下载各种坐标文件（bed，gtf 等），还提供了不同版本基因组坐标转换的工具 LiftOver。需要注意的是，UCSC 中的基因组序列与 NCBI 的内容有所不同，最大的一点就是 UCSC 的染色体有 chr 前缀，而 GRCh 没有 chr 前缀。此外 chr1 到 chr22，chrX，chrY 序列与 GRCh37 完全一致，但是线粒体序列稍微不一样，以及在重复区域 repeat region 有小写来表示，这点和 GRCh 不同。
       详细文档：

1. http://hgdownload.soe.ucsc.edu/downloads.html#hg38sequence

复制代码

1. #UCSC hg19
   
2. wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
   
3. #UCSC hg38
   
4. wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

复制代码

3.4 千人基因组
       国际千人基因组计划依托中国深圳华大基因研究院、英国桑格研究所、美国国立人类基因组研究所。“国际千人基因组计划”自 2008 年 1 月 22 日启动，测序的总任务为 1200 个人（故称为千人基因组计划） ，旨在绘制迄今为止最详尽、最有医学应用价值的人类基因组遗传多态性图谱。深圳华大基因研究院作为发起单位之一，不仅承担了 400 个黄种人全基因组样本的测序和分析工作，而且还帮助完成了非洲人群的全部测序和分析任务。
       “千人基因组计划”将测序的人群包括：尼日利亚伊巴丹区域的约鲁巴人；居住于东京的日本人；居住于北京的中国人；美国犹他州的北欧和西欧人后裔；肯尼亚 Webuye 的 Luhya 人和 Kinyawa 的 Maasai 人；意大利的 Toscani 居民；居住于休斯顿的 Gujarati 印第安人；居住于丹佛的中国人；居住于洛杉矶的墨西哥人后裔；居住于美国西南部的非洲人后裔。

1. https://www.internationalgenome.org/

复制代码

      数据下载页：

1. https://www.internationalgenome.org/data#download

复制代码

3.4.1 b37 版本
       千人基因组计划和 NCBI 基因组通用，其中来自千人基因组计划第一阶段的参考基因组，命名为 b37，其实就是 NCBI 的 GRCh37 版本,但是包含了线粒体 MT 序列（the rCRS mitochondrial sequence），unlocalized sequences 和 unplaced sequences 以他们的检索号命名，但是不包含 alternate loci。
       数据可以从 NCBI 站点也可以从 EBI 站点进行下载。

       NCBI千人：

1. ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp

复制代码

      EBI千人：

1. ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/

复制代码

1. #1000 genome   https://www.internationalgenome.org/data
   
2. #ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/
   
3. ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz
   
4. #decoy版本，带EBV病毒
   
5. http://www.cureffi.org/2013/02/01/the-decoy-genome/

复制代码

3.4.2 hs37d5 版本
       hs37d5 来自于千人基因组计划第二阶段的参考基因组，可以理解是 b37 的升级版，在 1000 Genome 第二阶段使用。hs37d5 包含了 b37 的数据，以及 b37 在 GRCH37 的基础上进行命名和坐标系统规范，并增加了一条病毒序列（疱疹病毒），一条 decoy 序列（并且在 Y 染色体上把 X,Y 染色体的同源区 mark 成了 N。
       hs37d5=b37+人类疱疹病毒序列+ “decoy” sequence+ HuRef、BAC 或者质粒克隆和NA12878。

1. #ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/
   
2. ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz

复制代码

3.5 GATK
       gatk 在官网提供了一个 resource bundle，里面包含了所需要的很多数据，如果使用 gatk 软件，最好把这些数据下载下来。
这些数据提供 FTP 和 Google Cloud bucket 两种下载方式。

1. https://gatk.broadinstitute.org/hc/en-us

复制代码

1. lftp ftp://[email protected]/bundle/
   
2. ls
   
3. cd buldle
   
4. #mget或者mirrors进行下载
   
5. mirror hg38

复制代码

3.6 不同版本基因组之间坐标转换
       由于人全基因组存在多个版本，不同版本之间又有差别，主要是坐标位置的不同。每次基因组版本的升级，比如从 hg18 到 hg19，再到 hg38，坐标系统已经不一样，如果选择一个版本，需要与下游的注释数据库 ID 相匹配，例如 hg19 的版本需要选择对应的数据库版本，否则就会得到错误的结果。例如去公共数据库查询频率，位置等信息时，都要对应到使用的参考基因组查询相关信息，才能保持信息的一致性。虽然 2013 年发布了 GRCh38 基因组版本，要比 GRCh37 的基因组版本有更高的质量，但目前 GRCh37 版本使用时间更长，注释信息更全。并且 GRCh38 每年还会在不改变序列和坐标的情况下发布一些更新补丁 Patches。
       详情请见下面网址：

1. https://www.ncbi.nlm.nih.gov/grc/help/patches/

复制代码

      对于不同版本基因组坐标不同的问题，可以使用一些工具进行转换，主要包括网页版工具LiftOver，remap，以及命令行工具 crossmap。

       进行 liftover 需要一个 chain data，用于描述新旧 build 之间的差异。

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

复制代码

1. #
   
2. conda install crossmap -y
   
3. CrossMap.py
   
4. # Get the chain file that maps from hg19 to hg38.
   
5. # 下载 chain data
   
6. wget
   
7. http://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.
   
8. gz
   
9. # Get a test data file that will be remapped.
   
10. # bed 文件
    
11. wget http://data.biostarhandbook.com/data/ucsc/test.hg19.bed
    
12. # Run the remapping process.
    
13. # 进行 remap
    
14. CrossMap.py bed hg19ToHg38.over.chain.gz test.hg19.bed test.hg38.bed

复制代码

生信喵

四、注释数据库下载
       通过将测序样品与参考基因组比对得到的突变定位到基因组上，例如定位到具体的基因区，基因间区，进行 genotype 分型，计算等位频率等，这些都需要与已知数据库进行比对。目前已经累计了很多关于人基因组突变的信息。例如 dbSNP，dbVAR，cgi69ExAC.vcf.gz（broadinstitute 提供的外显子联盟），Cosmic_v73.ann.vcf.gz （癌症突变信息集），finalTCGA.vcf.gz （TCGA 计划癌症相关），1000g-ph3v5.gff.gz（千人基因组计划），ESP6500（Variants from the Exome Sequencing Project (ESP)），瓶中基因组计划，deCODE 计划，ENCODE 计划，还有各个国家级的基因组计划等。
4.1 dbSNP
       dbSNP: database of SNP，单核苷酸多态性数据库。是由 NCBI 与人类基因组研究所（National Human Genome Research Institute）合作建立的，它是关于单碱基替换以及短插入、删除多态性的资源库。因为开发 dbSNP 是为了补充和辅助 GenBank, 所以它包含了来自任何生物体的核苷酸序列。
SNP 数据库的数据内容分为两类：一类是提交数据，即观察所得的原始序列变异；另一类是计算内容，即通过对原始提交数据的计算在每个“build”周期中产生的内容。

1. https://www.ncbi.nlm.nih.gov/SNP/

复制代码

1. #dbSNP hg38
   
2. wget ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_dbSNP_all.vcf.gz

复制代码

4.2 clinvar
       ClinVar 是 NCBI 主办的与疾病相关的人类基因组变异数据库。它的强大在于整合了 dbSNP、dbVar、Pubmed、OMIM 等多个数据库在遗传变异和临床表型方面的数据信息，形成一个标准的、可信的遗传变异-临床相关的数据库。当前数据库可支持 XML、VCF 以及制表符定界文件格式的文件。ClinVar 同时支持在线和下载到本地两种形式。
       ClinVar 是一个公开的数据库，其中收集了与疾病相关的遗传变异。这一数据库由美国国立卫生研究院 2013 年为了生物技术信息开发而构建，来自美国联盟医疗体系（Partners Healthcare）的临床遗传学家 Heidi Rehm 表示，到目前为止，已经从研究人员和其它数据库中获得了包含超过 125,000 份独特突变的临床注释。
       ClinVar 将基因突变对健康影响的多方面性质都考虑在内了，比如对于一个突变，这个研究组说它是良性，但另外一个研究组又认为它其实性质更加严重。而且 ClinVar 也有自己的分类，“可能致病性（likely pathogenic）”就是一种更清楚的定义和标准化。
工作机理：
       ClinVar 采用的是星标系统（star-based system），可以评估某个特定突变在疾病中的本来或者注释作用。四星级是最高级，也就是说这个突变的功能已经经过了多位专家的测评。这样详细审查过程的好处在于用户能信任三、四星突变的功能注释，Rhem 说，但是在 ClinVar 数据库中只有少量此类数据（3800 个）。
       更多时候是只有一星的突变，这通常只基于单个提出注释功能的研究成果，还有没有星的，也就是提交者没有提供解释标准和支持证据。ClinVar 工具的一个问题在于其数据库中大多数临床上重要的突变（83%）都是某个家庭中独一无二的，或者非常罕见的。

1. https://www.ncbi.nlm.nih.gov/clinvar/

复制代码

1. #ClinVar
   
2. wget ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_clinvar.vcf.gz

复制代码

4.3 dbVAR
       dbvar 数据库收录了来自多个研究项目的结构变异结构，在该数据库中，数据分为以下 3 个层级study,用 std 表示, 代表一个研究项目，该项目包含了一组结构变异分析结果，每个 study 有一个唯一的编号，如果来自 NCBI，则以 nstd 开头，如果来自 EBI,则以 estd 开头variant regions, 用 sv 表示，代表存在结构变异的基因组区域，同样有一个唯一的编号，如果来自 NCBI，则以 nsv 开头，如果来自 EBI,则以 esv 开头variant calls, 用 ssv 表示，代表一个具体的结构变异事件，包含了缺失，重复，CNV 等多种类型，同样有一个唯一的编号，如果来自 NCBI，则以 nssv 开头，如果来自 EBI,则以 essv 开头

1. https://www.ncbi.nlm.nih.gov/dbvar

复制代码

1. #dbVAR
   
2. wget ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_dbVar.gvf.gz

复制代码

4.4 其他数据库
       一些对于突变进行注释的软件，自带标准格式的数据库，流入 annovar，snpeff，vep 等，都带有各自的数据库。

1. annovar：http://annovar.openbioinformatics.org/en/latest/
   
2. snpeff：https://pcingola.github.io/SnpEff/
   
3. vep：https://github.com/Ensembl/ensembl-vep

复制代码