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