有沒有想過這個(gè)問題：做研究的時(shí)候該用hg19還是hg38基因組？
人類基因組版本現(xiàn)狀
對(duì)于同一個(gè)物種，數(shù)據(jù)庫(kù)中存在不同的基因組版本，以人類（Homo Sapiens）為例，UCSC基因組瀏覽器中有多個(gè)版本：Dec. 2013 (GRCh38/hg38)、Feb. 2009 (GRCh37/hg19)、Mar. 2006 (NCBI36/hg38)等，括號(hào)前面的是組裝（Assembly）日期。是不是有點(diǎn)驚訝！都2022年了，默認(rèn)的基因組還是9年前的，更令人驚訝的是，好多網(wǎng)站現(xiàn)在默認(rèn)使用的還是hg19，也就是13年前的基因組版本，更離譜的是，GEO數(shù)據(jù)庫(kù)中存在大量hg19，甚至hg18的數(shù)據(jù)集。那么在GEO數(shù)據(jù)挖掘的在日常工作中，經(jīng)常會(huì)遇見以下兩種場(chǎng)景：1）hg19 -> hg38，例如文獻(xiàn)中使用的是hg19，你自己的測(cè)序數(shù)據(jù)是hg38
2）hg38 -> hg19，例如你師兄的師兄留給你的數(shù)據(jù)是hg19，而文獻(xiàn)是hg38
解決方案：UCSC提供的一個(gè)工具liftover
Liftover簡(jiǎn)介
http://genome.ucsc.edu/cgi-bin/hgLiftOver
將不同版本的染色體上的位置一一對(duì)應(yīng)，官方定義為：This tool converts genome coordinates and genome annotation files between assemblies.
不難想象，該工具至少需要3個(gè)參數(shù)：hg19的坐標(biāo)文件，hg38的坐標(biāo)文件，以及兩者之間關(guān)系的數(shù)據(jù)庫(kù)文件（chain文件）
bw文件簡(jiǎn)介
bw文件是bigwig的縮寫，存儲(chǔ)了坐標(biāo)及對(duì)應(yīng)信號(hào)信息。然而，bw是一種二進(jìn)制文件，不能用liftover直接處理。因此，要將hg19基因組的bw文件轉(zhuǎn)成hg38，需要找個(gè)代理，這個(gè)代理就是bedgraph文件，包含4列，例如
chr1 100 200 5.2
表示：1號(hào)染色體100到200區(qū)域中的信號(hào)是5.2
bedgraph格式可以被liftover直接處理。
圖1. 轉(zhuǎn)換示意圖
轉(zhuǎn)換代碼
前人栽樹后人乘涼，python 包CrossMap可以直接處理bw文件的轉(zhuǎn)化問題。
因此一行代碼的轉(zhuǎn)化過程如下：
1，安裝CrossMap
pip install CrossMap2，下載hg19-hg39轉(zhuǎn)化對(duì)應(yīng)的數(shù)據(jù)庫(kù)文件
http://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz
3，一行代碼轉(zhuǎn)化
CrossMap.py bigwig hg19ToHg38.over.chain input.bw output.bw然后就可以導(dǎo)入到IGV進(jìn)行查看和比較了。
當(dāng)然，也可以逐步進(jìn)行
bigWigToBedGraph input.bw input.bedGraph
liftOver input.bedGraph hg19ToHg38.over.chain input_hg38.bedgraph
fetchChromSizes hg38 > hg38.chrom.sizes
sort -k1,1 -k2,2n input_hg38.bedgraph > input_hg38.sorted.bedgraph
bedGraphToBigWig input_hg38.sorted.bedgraph hg38.chrom.sizes output.bw
其中bigWigToBedGraph、fetchChromSizes、bedGraphToBigWig都可以在UCSC下載
http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/
最后，讓我們來檢驗(yàn)下hg19和hg38的轉(zhuǎn)換效果吧：
圖2. Hg19可視化
圖3. Hg38可視化
一模一樣，肉眼看不出差別，說明結(jié)果正確。
注意：有些位點(diǎn)沒有對(duì)應(yīng)關(guān)系的話，會(huì)轉(zhuǎn)化失敗，因此最佳解決方案還是使用hg38基因組從原始數(shù)據(jù)重新處理。
回答開頭的問題：
現(xiàn)在包括UCSC、TCGA、Ensembl等大型數(shù)據(jù)庫(kù)均以hg38作為默認(rèn)基因組，因此，用hg38就對(duì)了，還在用hg19的研究者，請(qǐng)趕緊升級(jí)！
數(shù)據(jù)分析的時(shí)候，一定要看清楚，網(wǎng)上的數(shù)據(jù)到底是hg38還是hg19！因?yàn)槌舜笮蛿?shù)據(jù)庫(kù)外，其他的小型數(shù)據(jù)庫(kù)、網(wǎng)站經(jīng)常是發(fā)了文章就不再更新，甚至是發(fā)了文章6個(gè)月以后就找不到那種，一定要“三思而后行”，避免浪費(fèi)時(shí)間！
微生信助力發(fā)文章，谷歌引用590+，知網(wǎng)引用450+
