據(jù)可視化與peak質(zhì)控全指南)
1. 項目概述為什么IGV是組學數(shù)據(jù)可視化的第一塊敲門磚先交代一下背景。我最早接觸IGVIntegrative Genomics Viewer綜合基因組學查看器是在做ChIP-seq數(shù)據(jù)分析的時候。當時上游的call peak流程跑完了拿到了一個幾百甚至幾千個顯著富集區(qū)間的列表但光看數(shù)字根本判斷不了數(shù)據(jù)質(zhì)量好不好——比對率是達標的、peak數(shù)量也不少、motif分析也做了可這些區(qū)間到底是不是真的富集覆蓋度長什么樣有沒有明顯的假陽性這些問題只有把reads比對情況和peak區(qū)間疊加到參考基因組上親眼看一看才能確認。IGV就是干這個事的。它是一款輕量級、跨平臺Windows/macOS/Linux都能跑的桌面端基因組瀏覽器由Broad Institute開發(fā)。它最核心的價值就一句話把幾百GB的比對文件、幾十MB的peak區(qū)間文件、以及基因注釋信息整合到同一個視圖里讓你在幾秒鐘內(nèi)完成對某個基因組區(qū)域的“肉眼質(zhì)檢”。這篇內(nèi)容主要面向三類人一是剛跑完ChIP-seq、DAP-seq或ATAC-seq流程但不知道怎么把結(jié)果“看”出來的新手二是已經(jīng)被各種IGV報錯卡住、網(wǎng)上搜不到有效解決方案的實操黨三是需要在文章審稿階段給關(guān)鍵peak區(qū)間做截圖展示的同學。我盡量把從安裝到出圖的完整路徑講清楚重點放在最容易踩坑的細節(jié)上。2. 核心概念與數(shù)據(jù)格式2.1 三種測序技術(shù)的數(shù)據(jù)特點差異ChIP-seq、DAP-seq和ATAC-seq雖然最后都落到“看峰”這一步但三種技術(shù)的數(shù)據(jù)特征有明顯差異理解這些差異才能選對可視化參數(shù)。ChIP-seq染色質(zhì)免疫沉淀測序針對的是蛋白-DNA互作抗體富集特定轉(zhuǎn)錄因子或組蛋白修飾所以數(shù)據(jù)特點是peak比較尖銳對于轉(zhuǎn)錄因子或者呈現(xiàn)較寬的富集域?qū)τ贖3K27ac這類組蛋白修飾。看這類數(shù)據(jù)時IGV中通常需要同時查看input或IgG對照因為很多區(qū)域存在非特異富集。DAP-seqDNA親和純化測序是體外表達轉(zhuǎn)錄因子蛋白然后與基因組DNA孵育沒有染色質(zhì)背景信噪比通常比ChIP-seq更高。但它有一個特殊問題因為沒有體內(nèi)交聯(lián)步驟peak邊界常常不夠清晰而且容易在基因組高拷貝區(qū)域出現(xiàn)假陽性所以檢查時尤其要關(guān)注重復區(qū)域的reads覆蓋是否均勻。ATAC-seq轉(zhuǎn)座酶可及染色質(zhì)測序關(guān)注的是開放染色質(zhì)區(qū)域數(shù)據(jù)特點一是peak寬度變化大從幾百bp到幾kb二是需要特別留意線粒體reads比例是否過高三是TSS區(qū)域應(yīng)該有明顯的reads堆積信號??梢暬瘯r建議把不同的生物學重復樣本同時載入用于快速核對peak的可重復性。2.2 看懂IGV中的文件格式IGV里最常打交道的有三種文件格式很多新人一上來就混淆。BAM是比對文件也就是reads比對回參考基因組的結(jié)果。染色質(zhì)測序的數(shù)據(jù)可視化基礎(chǔ)就是它。BAM需要配套BAM.BAI索引文件才能被IGV正確讀取這個索引就是BAM文件的“目錄”沒有它IGV就像去圖書館卻沒有檢索系統(tǒng)只能傻等甚至直接報錯。BED是最常用的peak區(qū)間格式由call peak軟件如MACS2輸出。標準的BED最少包含三列染色體號、起始位置、終止位置。IGV會自動把BED顯示成黑色條塊。如果文件帶有score列IGV還可以用顏色深淺來表示信號強弱。BWBigWig是覆蓋度文件本質(zhì)上是把每個基因組位置的reads堆疊深度轉(zhuǎn)換成連續(xù)信號。它的優(yōu)勢是文件體積小、加載快適合同時看多個樣本的整體信號分布。IGV中BW信號默認顯示為“山峰”狀圖形這也是很多人最直觀感受到的“峰”的視覺印象。另外還有一個常見問題很多人分不清BED和narrowPeak的區(qū)別。narrowPeak是MACS2專門輸出的peak格式比BED多了10列包含peak中心位置、富集倍數(shù)、qvalue等信息。IGV可以直接識別narrowPeak但如果你想要的是最簡展示效果建議用cut命令把前3列或前6列提取出來轉(zhuǎn)成標準BED。2.3 參考基因組的選擇邏輯參考基因組選錯是IGV使用中最隱蔽的坑。IGV讀取BAM文件時會嚴格檢查BAM頭部記錄的染色體名稱是否與當前選中的基因組版本一致。如果你用的BAM是hg19GRCh37比對的但IGV里選擇的是hg38GRCh38那么所有reads都會對不上位置顯示“No sequence available”或者直接空白。類似地植物數(shù)據(jù)中常見的是TAIR10擬南芥、IRGSP水稻日本晴等版本。這里有個判斷技巧從BAM文件的頭部信息或比對軟件的輸出了報告中確認具體版本號不要靠猜。注意有些標準化處理流程會直接把染色體名稱中的“chr”前綴去掉比如人類數(shù)據(jù)從chr1變成1這會導致IGV無法正確匹配染色體。一定要保持染色體命名規(guī)則前后一致。3. 實操過程從安裝到5分鐘出圖3.1 下載安裝與環(huán)境檢查IGV的安裝非常簡單但有幾個細節(jié)值得留意。去IGV官網(wǎng)下載頁面選擇對應(yīng)系統(tǒng)版本W(wǎng)indows用戶直接運行exe安裝包macOS用戶下載dmg文件后拖入Applications目錄即可Linux用戶需要下載壓縮包后解壓通過igv.sh腳本啟動。另外強調(diào)一點IGV要求本機裝有Java運行環(huán)境且不同版本的IGV對Java版本要求不同——老版本2.11及更早通常需要Java 8或11新版本2.12建議Java 17以上。啟動時如果報UnsupportedClassVersionError多半是Java版本不匹配導致的直接更新Java或改用IGV自帶JRE的版本可以解決。Linux服務(wù)器上如果沒有圖形界面也可以使用IGV的批處理模式——啟動時加-b參數(shù)指定一個批處理腳本能自動完成載入文件、跳轉(zhuǎn)到指定位置、截圖等操作。這個方法對批量出圖非常有用我會在后面詳細展開。3.2 數(shù)據(jù)載入的標準流程IGV啟動后首先在左上角的基因組下拉框中選擇正確的參考基因組版本。如果你的物種不是默認列表里的需要通過Genomes - Load Genome from File...導入fasta文件或.fasta.fai索引文件。注意導入自定義基因組后IGV需要建立索引耗時取決于基因組大小通常幾十秒到幾分鐘。接下來載入比對文件File - Load from File...選擇BAM文件IGV會自動找到同目錄下的.bai索引。如果BAM和BAI不在同一目錄會報錯解決辦法稍后詳細說明。然后是peak文件同樣通過Load from File載入narrowPeak或BED。為了讓peak顯示更直觀可以在peak軌道的track名上右鍵選擇Expanded顯示模式讓每個peak顯示成獨立的條塊而不是緊湊堆疊在一起。最后是BW覆蓋度文件載入后IGV默認顯示為“Signal”模式。不同樣本可以用不同顏色區(qū)分——右鍵軌道選擇Set Track Color...方便做樣本間對比。3.3 快速導航到目標區(qū)域有幾種方式可以快速跳轉(zhuǎn)到感興趣的區(qū)域。最簡單的辦法是在IGV頂部搜索框直接輸入基因名比如TP53或染色體坐標比如chr3:178,916,000-178,920,000。搜索框支持模糊匹配輸入基因名后會在下拉框中提示匹配到的注釋轉(zhuǎn)錄本。第二種方式是從peak文件直接跳轉(zhuǎn)。在IGV的peak軌道上點擊任意一個peak條塊視圖會自動居中到該peak位置。如果你用的是narrowPeak格式且包含score信息鼠標懸停時還能看到富集倍數(shù)和qvalue。第三種方式是梭巡模式。菜單欄View - Preferences - Advanced中開啟Allow one-click navigation后在top的坐標軸上單擊可以直接跳轉(zhuǎn)到對應(yīng)位置。實操中我習慣的做法是先用grep或awk從narrowPeak文件中提取某個感興趣的peak坐標然后直接在搜索框中輸入坐標跳轉(zhuǎn)。這樣可以避免因為基因名有多個轉(zhuǎn)錄本注釋而跳轉(zhuǎn)到不想要的位置。3.4 常用展示模式調(diào)整IGV的可視化效果可以從幾個維度調(diào)整。首先是read顯示模式。默認情況下IGV對BAM軌道使用Squished模式也就是壓縮顯示。對ChIP-seq這類數(shù)據(jù)我通常調(diào)整為Expanded模式這樣可以看到每條read的實際比對方向和read結(jié)構(gòu)方便判斷鏈特異性或結(jié)構(gòu)變異信號。但要注意Expanded模式在reads很多的區(qū)域會非常占內(nèi)存加載可能變慢。第二條是展示reads的結(jié)構(gòu)信息。在BAM軌道右鍵選擇Color By - Read Strand可以用顏色區(qū)分正負鏈選擇Color By - Read Group可以區(qū)分多個樣本混合比對時的來源。對ATAC-seq數(shù)據(jù)尤其建議查看insert size分布幫助判斷核小體定位信號是否正常。第三條是peak信號范圍調(diào)整。BW軌道的縱軸范圍默認是自動縮放的但有時候自動縮放會導致信號看起來噪聲很大。右鍵選擇Set Data Range...手動設(shè)置一個合適的信號值范圍比如0到50然后保存為默認設(shè)置這能顯著提升出圖美觀度。另外IGV支持筆記Region Navigator功能。在坐標軸左側(cè)的梯形圖標上按住拖動可以標記自定義區(qū)域。標記后可以在Regions列表中快速跳轉(zhuǎn)方便對一批候選peak進行逐個審閱。3.5 出圖保存與文章級截圖技巧對發(fā)文章來說IGV截圖幾乎成了組學數(shù)據(jù)展示的標準動作。IGV可以直接File - Save Image...保存為PNG或SVG格式。這里有幾個讓截圖“看起來專業(yè)”的經(jīng)驗PNG適合直接提交到雜志系統(tǒng)分辨率建議設(shè)300dpi以上。SVG適合后續(xù)用AI或Inkscape二次編輯靈活性更高但需要矢量工具支持。在Save Image對話框中可以設(shè)置圖片寬度和縮放比例通常2倍縮放足夠清晰。為了美化截圖有幾個常規(guī)操作在View - Preferences - Tracks里關(guān)閉不必要的外顯子/內(nèi)含子結(jié)構(gòu)顯示只保留基因模型和關(guān)鍵軌道在BAM軌道右鍵選擇Collapsed模式讓大量reads聚合成堆疊條帶視覺上更干凈peak軌道用紅色或深色填充色BW信號用藍色較直觀將坐標軸字體調(diào)大保證縮小后仍能看清位置信息。實操心得截圖前一定要先View - Preferences - General中勾選Show base modification和Show all tracks這些基礎(chǔ)選項并確認Gene track處于可見狀態(tài)。很多初學者的截圖問題不是IGV不會用而是注釋軌道隱藏了導致圖片里看不到基因結(jié)構(gòu)。4. 常見問題與排查技巧實錄4.1 文件加載類報錯報錯1Could not find index file這個錯誤幾乎每個人都會遇到。IGV載入BAM文件時會在同一目錄下尋找同名.bai文件。如果你的.bai文件不存在或與BAM文件名不完全一致比如BAM叫sample.bam索引叫sample.bam.bai是正常的但如果叫sample.bai則IGV認不出來就會提示出錯。解決辦法分兩步先用samtools index sample.bam生成索引如果索引已存在但文件名不匹配用samtools index -b sample.bam重命名為sample.bam.bai再試。另一個可能被忽略的問題是BAM文件本身需要是排序過的——samtools sort排序后再建立索引。未排序的BAM即使有索引IGV加載時也會出現(xiàn)位置錯亂。報錯2No data is available for the selected region這個報錯信息很誤導人。它不一定是說這個區(qū)域沒數(shù)據(jù)而是說IGV找不到可以展示的track。最常見的場景是載入了BAM但沒載入任何注釋軌道或者參考基因組版本錯誤導致reads實際比對位置不在當前視圖。排查時先確認以下幾點BAM軌道是否處于可見狀態(tài)左側(cè)復選框勾上了嗎當前坐標是否在染色體的實際范圍內(nèi)BAM文件頭部是否包含SQ行以及對應(yīng)染色體名稱。最后提醒自己檢查參考基因組版本是否與比對版本一致。報錯3Track loaded but nothing displays有時候BAM載入了IGV不報錯但界面上一片空白。這種情況首先檢查軌道顯示范圍軌道是否存在但在縮放過小的視圖中把視圖縮放到單基因區(qū)域再看。如果仍無顯示用IGV底部的“Region”插件——右鍵點擊BAM軌道選擇View in GeneLists輸入一個已知高覆蓋區(qū)域如ACTB驗證數(shù)據(jù)可顯示性。4.2 展示異常類排查問題1所有reads顯示成灰色或全鏈同一方向大部分是因為染色鏈分析設(shè)置。在BAM軌道右鍵選擇Color By - Read Strand然后檢查BAM文件本身是否混入了未處理的adapter序列或者比對方向異常。ChIP-seq數(shù)據(jù)一般不存在鏈特異性所以read也不要刻意按鏈著色。問題2peak顯示在基因間區(qū)且覆蓋度極低這類peak通常是MACS2開默認參數(shù)跑出來的低質(zhì)量peak。先檢查BW信號是否在比對文件載入后正常顯示如果BW信號微弱但peak文件顯示大量區(qū)間說明call peak時可能使用了過低的qvalue閾值??梢杂胊wk $95篩選富集倍數(shù)更高的peak再看。問題3IGV載入大文件時內(nèi)存爆滿人類全基因組數(shù)據(jù)按30×覆蓋度BAM文件很容易超過100GB。IGV在載入時如果內(nèi)存設(shè)置過低會直接卡死。解決辦法是在啟動腳本或快捷方式的參數(shù)中添加-Xmx4g甚至-Xmx8g來提高最大堆內(nèi)存。對于只做peak檢查的場景還可以用samtools view -b截取特定染色體區(qū)間生成一個臨時小BAM快速載入查看不必每次都加載全基因組數(shù)據(jù)。4.3 常見問題速查表現(xiàn)象可能原因解決方案找不到.bai索引未生成索引或名稱不匹配samtools index生成標準名索引比對區(qū)域空白參考基因組版本不一致核對BAM頭部SQ信息與IGV選擇版本BAM能載入但無reads視圖區(qū)間過大或軌道被隱藏縮小坐標范圍檢查軌道左端復選框BW信號為全0或全滿數(shù)據(jù)范圍參數(shù)異常手動設(shè)置Data Range或重新生成BW文件啟動時報Java錯誤JDK版本不匹配更新到Java 17或使用自帶JRE的IGVpeak與gene重疊很少peak質(zhì)量差或參考注釋不全檢查注釋來源改用最新GTF注釋截圖文字過小保存參數(shù)設(shè)置不當提高分辨率至300dpi并開啟2×縮放自定義基因組無法載入fasta缺少.fai索引用samtools faidx生成索引4.4 批量出圖讓IGV自動化處理多個區(qū)間如果文章需要展示幾十個peak的IGV截圖手動一個個跳轉(zhuǎn)再保存顯然是低效的。IGV的批處理模式非常適合這個場景。寫一個簡單的批處理腳本# 建立IGV批處理命令文件 cat igv_batch.txt EOF snapshotDirectory /path/to/output_dir genome hg38 load /path/to/sample.bam load /path/to/sample_peaks.narrowPeak load /path/to/sample.bw goto chr3:178,916,000-178,920,000 snapshot TP53_peak.png goto chr7:55,025,000-55,030,000 snapshot EGFR_peak.png exit EOF # 運行IGV批處理 igv -b igv_batch.txtsnapshotDirectory指定截圖輸出目錄goto命令直接跳轉(zhuǎn)到目標區(qū)域snapshot保存圖片。這個腳本可以有多個goto-snapshot組合一次運行生成幾十張圖。對于大多數(shù)情況IGV默認使用當前加載的track設(shè)置所以先手動把BAM、BW、peak的樣式調(diào)整好再在批處理腳本中只load文件和goto坐標就能保持統(tǒng)一風格。實操心得批量運行前務(wù)必要在腳本中提前snapshotDirectory否則IGV會默認把截圖保存到當前工作目錄最后找起來非常麻煩。5. 三種數(shù)據(jù)類型的經(jīng)典IGV可視化實操案例5.1 ChIP-seq看轉(zhuǎn)錄因子結(jié)合位點和input對照場景你想展示某個轉(zhuǎn)錄因子在一個靶基因啟動子區(qū)域的結(jié)合信號。除了特異性抗體樣本還必須展示input對照。實際操作中我通常先載入三個trackIP樣本的BAM/BW、input樣本的BW、以及MACS2輸出的narrowPeak。查看時注意input的BW信號如果在一個區(qū)域也普遍偏高說明這個區(qū)域本身是開放染色質(zhì)或高拷貝區(qū)IP樣本在該區(qū)域的peak可能是背景而不是真實結(jié)合。對比手段除了肉眼觀察還可以在narrowPeak上懸停查看fold enrichment值。在一個區(qū)域如果IP信號高但peak文件標注的差異倍數(shù)只有1.5左右那很可能是不太可靠的peak需要在文章中隱掉或者標注低置信度。補充技巧當IP樣本的BW信號顯示出一條平滑的“山脊”同時narrowPeak的區(qū)間覆蓋該“山脊”的頂部區(qū)域這是典型的高質(zhì)量結(jié)合峰形。如果peak出現(xiàn)在信號平坦區(qū)域甚至信號為零的區(qū)域那就是peak calling的假陽性一定要在最終結(jié)果中去掉。5.2 DAP-seq注意體外實驗的景觀型富集場景DAP-seq因為沒有體內(nèi)復雜背景peak通常更干凈但也更容易在轉(zhuǎn)座子等重復區(qū)域產(chǎn)生虛假結(jié)合。如果你看一個DAP-seq結(jié)果的IGV截圖時發(fā)現(xiàn)大量reads堆積在轉(zhuǎn)座子注釋區(qū)域首先別懷疑實驗做壞更可能是在比對階段沒有對多拷貝區(qū)域做特殊處理。處理辦法比對時加--multimapping過濾或使用MAPQ10的閾值過濾多比對reads在IGV中也可以通過Filter Alignments功能過濾低質(zhì)量reads后再查看。DAP-seq還有另一個典型特征個別轉(zhuǎn)錄因子的結(jié)合基序非常短導致peak很短而且信號集中看起來像尖銳的單峰甚至是SNP信號。這時候放大到base level看reads信息如果發(fā)現(xiàn)只有少數(shù)reads的比對質(zhì)量很低那基本可以判定不是真正的結(jié)合信號。5.3 ATAC-seq驗證開放區(qū)域、檢查核小體周期性信號場景ATAC-seq的可視化有兩個重點一是看TSS附近的開放狀態(tài)二是評估核小體排布。建議載入兩個trackreads的BAM和MACS2輸出的broadPeak。IGV中縮小到TSS±2kb范圍可以看到明顯的reads堆積和中間缺口——這是開放染色質(zhì)的特征。更進一步將坐標縮短到±400bp如果看到所有reads都集中在中心說明該位點主要是開放區(qū)沒有穩(wěn)定核小體結(jié)合。ATAC-seq還有一個需要留意的信號如果BAM軌道按insert size著色Color By - Insert Size你應(yīng)該能看到小于100bp的短片段主要富集在真正的開放區(qū)域而大于150bp的長片段在核小體占位區(qū)域有信號。這個檢查能有效判斷數(shù)據(jù)中開放區(qū)域和核小體區(qū)域的相對信號比值是否合理。6. 進階操作讓IGV出圖更高效的小技巧6.1 自定義會話管理與配置保存當你同時分析多個樣本時手動重復載入文件和設(shè)置樣式會浪費大量時間。IGV支持保存會話文件.xml格式在當前配置滿意時選擇File - Save Session...下次直接Load Session...即可恢復所有軌道和視圖設(shè)置。會話文件是一個純文本XML里面會記錄已加載的文件路徑、軌道樣式、顏色、縮放位置等信息。對于同一項目下的所有樣本可以存一份標準的會話模板載入新數(shù)據(jù)時只替換文件路徑即可非常方便。6.2 利用IGV的Region Navigator批量審閱peakIGV右上角有一個Region Navigator面板支持導入自定義區(qū)域列表比如所有高置信度peak的中心±1kb坐標。用這個面板可以快速在一個個區(qū)域之間切換同時保持固定的顯示寬度非常適合做人工peak審閱。做法是準備一個簡單的BED格式文件三列分別為chr、start、end然后通過Region Navigator的“Import Regions”功能導入。之后用鍵盤上下鍵就可以在區(qū)域間切換結(jié)合箭頭縮放的設(shè)置可以快速瀏覽幾百個peak而不用動手輸入坐標。6.3 帶狀圖和sashimi圖對RNA-seq或剪接分析有需求的同學可能見過IGV的sashimi plot功能。這個是展示junction reads的好工具不過對ChIP-seq/ATAC-seq而言用到的場景較少。類似地IGV的Alternate View帶狀圖模式適合查看結(jié)構(gòu)變異如果你用的是長reads測序數(shù)據(jù)或需要查看大片段缺失可以切換到這個模式看reads如何在兩個距離很遠的區(qū)域間連接。這在ATAC-seq或ChIP-seq的質(zhì)控中偶爾會遇到比如線粒體reads比例異常或者外源DNA污染時。6.4 不同系統(tǒng)下的效率對比與推薦以我個人的使用習慣日常探索性分析在本地圖形界面IGV中完成批量出圖會用到批處理腳本而真正高并發(fā)的質(zhì)檢流程會考慮在服務(wù)器上使用IGV的Web服務(wù)接口進行無頭渲染。不過對于大多數(shù)做ChIP-seq/DAP-seq/ATAC-seq的項目來說本地IGV完全夠用尤其現(xiàn)在SSD普及之后即使50GB級別的大BAM文件跳轉(zhuǎn)和縮放都相當流暢。7. 寫在最后IGV的使用門檻其實很低絕大多數(shù)人真正卡住的原因不是軟件有多復雜而是對數(shù)據(jù)格式和邏輯理解不透徹。掌握了BAM/BW/BED這幾類文件的關(guān)系理解了參考基因組版本一致性的重要性剩下就是在實踐中不斷積累一些“看數(shù)據(jù)”的手感。我個人在帶新人時最推薦的路徑是先用一個已發(fā)表文章中的公共數(shù)據(jù)GEO上非常多完整跑一遍從載入到出圖的流程同時對照MACS2輸出的峰值列表逐一審視等有了“正常數(shù)據(jù)長什么樣”的直覺后再切換到自己的數(shù)據(jù)做質(zhì)檢馬上就能發(fā)現(xiàn)潛在問題。這種數(shù)據(jù)敏感度是IGV練出來的核心能力比任何參數(shù)技巧都重要。如果看完這篇內(nèi)容你也能獨立完成一次“從BAM到發(fā)表級截圖”的完整流程那這篇文章的價值就達到了。當然組學數(shù)據(jù)分析中還有更多細節(jié)等待你逐一探索愿你少踩一些我當年踩過的坑。