1. 如何選擇參考基因組？——李恒2017年年底博客
1.1 三種選擇
- 如果比對(duì)到GRCh37/hg19搀崭，使用：
ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz
- 如果比對(duì)到GRCh37/hg19，并且認(rèn)為包含
decoy序列能夠更準(zhǔn)確地進(jìn)行變異檢測(cè)咸这，使用：
- 如果比對(duì)到GRCh38/hg38苍姜，使用：
1.2 其他版本的GRCh37/GRCh38有什么問(wèn)題？
- 包含ALT contig序列
ALT contigs are large variations with very long flanking sequences nearly identical to the primary human assembly. Most read mappers will give mapping quality zero to reads mapped in the flanking sequences. This will reduce the sensitivity of variant calling and many other analyses. You can resolve this issue with an ALT-aware mapper, but no mainstream variant callers or other tools can take the advantage of ALT-aware mapping.
個(gè)人理解：上面說(shuō)的是一個(gè)多比對(duì)問(wèn)題曼月。新版的bwa可以處理ALT-aware的比對(duì)
ALT contig是什么谊却？
To better capture variation in the human genome across the world;
it(hg38最初版) contains more copies of some loci than hg19(最初版).
Patches（補(bǔ)度嶂纭）是什么哑芹？
patch releases do not change any previously existing sequences; they simply add new sequences for fix patches or alternate haplotypes that correspond to specific regions of the main chromosome sequences. For most users, the patches are unlikely to make a difference and may complicate the analysis as they introduce more duplication. (patches和ALT contigs的概念有些重合)
Alternate loci scaffolds:
- a scaffold that provides an alternate representation of a locus found in the primary assembly. These sequences do not represent a complete chromosome sequence although there is no hard limit on the size of the alternate locus; currently these are less than 1 Mb. These could either be NOVEL patch sequences, added through patch releases, or present in the initial assembly release.
- format: chr{chromosome number or name}_{sequence_accession}v{sequence_version}_alt
- e.g. chr6_GL000250v2_alt
所以ALT contig序列是為了反映人群多態(tài)性的一段替補(bǔ)序列，和原染色體位置對(duì)應(yīng)的序列之間有一定的差異捕透。放在ref中的隱患是人為增加了重復(fù)序列聪姿。
- 用很長(zhǎng)的N間隔這些ALT contig序列
增加了不必要的ref的size - 存在多處定位的參考序列
比如X染色體上的一些區(qū)域在Y染色體中也有(pseudoautosomal region, PAR)。在標(biāo)準(zhǔn)的流程中乙嘀，這些區(qū)域是沒(méi)法進(jìn)行變異檢測(cè)的（怎么理解末购？上面第一點(diǎn)也提到了）。正確的做法是將Y染色體中的這類區(qū)域進(jìn)行mask虎谢。 - 沒(méi)有使用rCRS的線粒體序列
rCRS在群體遺傳學(xué)中廣泛用到盟榴，但是官方的GRCh37的線粒體序列比rCRS長(zhǎng)2bp，構(gòu)建線粒體的系統(tǒng)發(fā)生樹應(yīng)該使用rCRS婴噩。GRCh38使用的是rCRS擎场。 - 將模棱兩可的堿基編碼為N
http://biocorp.ca/IUB.php - 使用accession numbers而不是chrXXX作為染色體名稱
- 不包含未精確定位的contig序列
hg19/chromFa.tar.gz from UCSC: 1, 3, 4 and 5.
hg38/hg38.fa.gz from UCSC: 1, 3 and 5.
GCA_000001405.15_GRCh38_genomic.fna.gz from NCBI: 1, 3, 5 and 6.
Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz from EnsEMBL: 3.
Homo_sapiens.GRCh38.dna.toplevel.fa.gz from EnsEMBL: 1, 2 and 3.
共性問(wèn)題是3
翻譯參考[1]
1.3 針對(duì)1.1前兩條的說(shuō)明
human_g1k_v37.fasta包含了人的24條染色體水平的序列，線粒體序列几莽，以及沒(méi)有定位的contig序列迅办。
hs37d5.fa包含了以上，以及NC_007605章蚣、hs37d5：
NC_007605是人皰疹病毒的序列站欺，Human herpesvirus 4 type 1，也稱EBV；
hs37d5就是decoy序列矾策，由很多短序列間隔一定數(shù)量的N組成磷账。
這些短序列是什么？有什么作用蝴韭？
這些短序列可以簡(jiǎn)單認(rèn)為是來(lái)源于人但是ref里面不包含（組裝技術(shù)的局限）够颠。怎么找到的？對(duì)比已有的克隆序列和ref序列確定的榄鉴。參考基因組并不完整履磨，總有一些區(qū)域不包含在已有的參考基因組中，這種情況下庆尘，來(lái)源于這些區(qū)域的reads要么比對(duì)不上剃诅，要么錯(cuò)誤地比對(duì)到其他地方（造成假陽(yáng)性）。如果能把那些已經(jīng)確定的但是ref沒(méi)有的序列加上去驶忌，就能有所改善矛辕。
EBV序列和decoy序列的作用類似，都是盡可能讓reads比對(duì)到它真實(shí)的地方付魔。只不過(guò)EBV序列不屬于human genome聊品，但細(xì)胞里面可能含有，提取DNA測(cè)序可能就測(cè)到了几苍。
參考[2]
1.4 個(gè)人搜集
Ensembl
可以下到最新版
ftp://ftp.ensembl.org/pub/release-98/fasta/homo_sapiens/dna/
ftp://ftp.ensembl.org/pub/release-98/gtf/homo_sapiens/
GATK
https://software.broadinstitute.org/gatk/download/bundle
包括SNP, InDel這類為變異檢測(cè)提供參考的文件翻屈。
現(xiàn)在好像進(jìn)不了
NCBI
最新版
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13
2. 認(rèn)識(shí)GRCh38/hg38
- 包含alternate contigs，能反映群體的單體型
比如HLA alternate妻坝，這段序列多態(tài)性非常高伸眶，一個(gè)ref不一定具有代表性。 - The GRCh38 analysis set hard-masks regions and provides decoy contigs for optimal read mapping（見(jiàn)下文NCBI部分）
- The challenge alternate contigs presents is a familiar one
- Latest versions of BWA-MEM handle GRCh38 alternate contig mappings
這一點(diǎn)是說(shuō)最新版的bwa可以支持比對(duì)到這類ALT區(qū)域了刽宪，針對(duì)李恒說(shuō)的1.2第一點(diǎn)厘贼。 - Alt-handling requires the SAM format ALT index file
- New Tutorial#8017 shows how to map to GRCh38 with alt-handling and then some
- Simulate read mapping for your favorite alternate haplotype
- Our production workflow for single sample variant calling on GRCh38 is public and uses shiny new features
- Finally, there is no better time than now to start learning WDL
參考[3]
- Inclusion of model centromere sequences
3. Masked, soft-masked and unmasked genomes
Ensembl上的參考基因組有以上三種形式，分別表示
- 將重復(fù)區(qū)域和低復(fù)雜度區(qū)域替換為N
信息缺失 - 將重復(fù)區(qū)域和低復(fù)雜度區(qū)域替換為小寫
和第三種沒(méi)什么區(qū)別圣拄，因?yàn)檐浖笮懸恢聦?duì)待嘴秸。
But it should be noted that most of the alignment tools do not take into account soft-masked regions, for example BWA, tophat, bowtie2 tools always use all bases in alignment weather they are in lowercase nucleotide or not. - 不做處理
推薦用這種
參考[4]
4. UCSC
以這個(gè)hg38文件夾為例，在這個(gè)文件夾中有幾個(gè)小文件夾：
4.1 initial/
GRCh38最初的版本庇谆，包含the original alternate sequences (261) and no fix sequences岳掐，和NCBI上的GCA_000001405.15有區(qū)別嗎？
hg38.chromFa.tar.gz - The assembly sequence in one file per chromosome.
Repeats from RepeatMasker and Tandem Repeats Finder (with period of 12 or less) are shown in lower case; non-repeating sequence is shown in upper case.
hg38.fa.gz - "Soft-masked" assembly sequence in one file.
Repeats from RepeatMasker and Tandem Repeats Finder (with period of 12 or less) are shown in lower case; non-repeating sequence is shown in upper case.
4.2 latest/
GRCh38最新的版本族铆，注意這里的最新是UCSC的最新岩四，可能NCBI更新更快。比如此時(shí)UCSC最新是GRCh38.p12哥攘，但NCBI已經(jīng)到GRCh38.p13了剖煌。
還有幾個(gè)重要文件：
4.3 mrna.fa.gz
Human mRNA from GenBank. This sequence data is updated regularily via automatic GenBank updates.
4.4 refMrna.fa.gz
RefSeq mRNA from the same species as the genome. This sequence data is updated regularily via automatic GenBank updates.
4.3和4.4可能是一樣的
4.5 標(biāo)準(zhǔn)的基因組序列文件
For most users, this will be the file "latest/hg38.fa.gz" in this directory.
http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/latest/
4.6 分析集序列文件
if you need a genome file for alignment or variant calling, please read the section "Analysis set" below.
http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/analysisSet/
下面兩個(gè)較上面多了alt-scaffolds材鹦；2,4是soft-masked的；第2個(gè)文件同the NCBI file GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
4.7 與NCBI ENSEMBL的比較
主要染色體的序列在幾個(gè)數(shù)據(jù)庫(kù)中是一樣的耕姊，但染色體名稱不同桶唐。
"chr1" at UCSC, "NC_000001.11" at NCBI, and "1" at the ENSEMBL;
soft-masked的堿基不完全一樣，因?yàn)槿哌\(yùn)行Repeatmasker設(shè)置參數(shù)不一樣
5. NCBI
最新版是GRCh38.p13茉兰，GCF_000001405.39
同樣的尤泽，包含主要染色體，未定位的scaffold序列规脸，patch和alt contig序列坯约。但是patch和alt contig序列是直接追加的，沒(méi)有在兩端加N莫鸭，因此解壓出來(lái)3.1G左右闹丐。
gff3文件信息很豐富，條目包括：CDS, enhancer, exon, gene, lnc_RNA, miRNA, rRNA, scRNA, snRNA, snoRNA, tRNA, TATA_box等等被因。
下面那個(gè)紅框是專門為分析流程準(zhǔn)備的文件卿拴，包含bwa, hisat2, bowtie等的索引，也就是分析集文件梨与，如下：
6. ensembl
官網(wǎng)：http://asia.ensembl.org/index.html
ftp://ftp.ensembl.org/pub/release-98/gff3/homo_sapiens/
這幾個(gè)有什么區(qū)別堕花？
- 98.chr.gtf.gz 只包含主要染色體（有MT）的注釋
- 98.gtf.gz 包含主要染色體以及一些沒(méi)有定位的染色體片段上的注釋信息
- 98.chr_patch_hapl_scaff.gtf.gz 除了上面的，還有patches和ALT contig的注釋
條目包含：CDS, exon, 5'-UTR, 3'-UTR, stop_codon等等常見(jiàn)的條目粥鞋。
ftp://ftp.ensembl.org/pub/release-98/fasta/homo_sapiens/dna/
6.1 toplevel & primary_assembly
toplevel
This includes chromsomes, regions not assembled into chromosomes and N padded haplotype/patch regions.
這個(gè)文件解壓出來(lái)有60G缘挽，遠(yuǎn)大于3G，增加了很多無(wú)用的N陷虎。比如一個(gè)MHC基因座（位于chr6）的多態(tài)性序列到踏，兩端會(huì)增加很多N杠袱，直到序列約等于chr6的長(zhǎng)度尚猿，再添加到ref中。
primary_assembly
Primary assembly contains all toplevel sequence regions excluding haplotypes and patches. This file is best used for performing sequence similarity searches where patch and haplotype sequences would confuse analysis. If the primary assembly file is not present, that indicates that there are no haplotype/patch regions, and the 'toplevel' file is equivalent.
7. 思考
盡管在1.中楣富，李恒老師列出了一些ref的不足凿掂，但實(shí)際分析中選不同的參考基因組影響不太大（他的博客也提到了這一點(diǎn)）。另外纹蝴，NCBI和ENSEMBL更新快于UCSC庄萎。
盡管patches和ALT contigs會(huì)不斷更新，GRCh的次版本（比如GRCh38.p13）也會(huì)隨之更新塘安，但主要染色體的gff/gtf應(yīng)該是不變的糠涛，因此不用擔(dān)心換了ref的次版本而找不到注釋文件，因?yàn)榇蠖鄶?shù)時(shí)候我們只看主要染色體上的注釋兼犯。
hg38提及decoy比較少忍捡，可能組裝過(guò)程中已經(jīng)整合進(jìn)去了一部分集漾，另外在NCBI分析集中提到了decoy（見(jiàn)上文）
7.1 需要探究的問(wèn)題？
在top level中砸脊，patch和ALT contig是如何存放的具篇？
每一個(gè)patch或者alt contig序列都是前后增加了很多N，然后追加到primary中凌埂。
NCBI和ENSEMBL是否都含有ncRNA的注釋驱显？（指第3列，ENSEMBL第9列也有miRNA, lncRNA等信息）
NCBI有瞳抓，ENSEMBL沒(méi)有埃疫。
或者這樣問(wèn)：它們的gtf有什么區(qū)別？NCBI的條目更多孩哑。
另外：NCBI上面的染色體命名不好認(rèn)熔恢，基因、轉(zhuǎn)錄本的ID沒(méi)有ENSEMBL全臭笆。
gencode的注釋和上文的gtf有關(guān)嗎憾朴？
- GENCODE最新版(Release 32)的注釋就是ENSEMBL最新版的默認(rèn)注釋集(GRCh38.p13)。另外坪蚁，GENCODE還提供了lncRNA, tRNA的gtf/gff文件蒙兰，ENSEMBL上面沒(méi)有這些單獨(dú)的文件，它是將type信息全都放在第9列茵乱。
- 而USCS將GENCODE作為一個(gè)子集track：
- NCBI的注釋是它自己的茂洒，叫“NCBI RefSeq”。在UCSC中也能下載這個(gè)子集track的gff/gtf瓶竭。
8. reference
[1] https://lh3.github.io/2017/11/13/which-human-reference-genome-to-use
[2] https://www.biostars.org/p/73100/
[3] https://software.broadinstitute.org/gatk/blog?id=8180
[4] https://genestack.com/blog/2016/07/12/choosing-a-reference-genome/
