1. 測(cè)序數(shù)據(jù)可視化的核心痛點(diǎn)與整體思路做測(cè)序分析的人都有一個(gè)共同的痛點(diǎn)手里拿到的是幾十GB的BAM文件里面全是比對(duì)到參考基因組的reads但你想跟別人展示這個(gè)區(qū)域覆蓋度很高或者這個(gè)樣本在某個(gè)基因上有明顯的表達(dá)信號(hào)時(shí)總不能把BAM文件甩給對(duì)方讓他自己用IGV慢慢加載。BAM文件是二進(jìn)制格式記錄的是每一條read的比對(duì)位置、比對(duì)質(zhì)量、CIGAR字符串這些信息信息量巨大但極其不適合做全基因組層面的快速瀏覽和比較。這就是bigWig格式存在的意義。bigWig是一種索引化的二進(jìn)制格式專門用來存儲(chǔ)連續(xù)的基因組信號(hào)數(shù)據(jù)比如覆蓋深度、GC含量、保守性得分等。它最大的優(yōu)勢(shì)是支持隨機(jī)訪問——你打開一個(gè)bigWig文件想看chr1:1000000-2000000這個(gè)區(qū)域的信號(hào)它不需要從頭掃描整個(gè)文件而是通過索引直接跳到對(duì)應(yīng)位置讀取數(shù)據(jù)。這個(gè)特性讓它在基因組瀏覽器如UCSC Genome Browser、IGV中加載速度極快幾百M(fèi)B的bigWig文件可以秒開而同樣信息量的BAM文件可能需要幾十秒甚至更久。從BAM到bigWig的轉(zhuǎn)換本質(zhì)上是一個(gè)信息降維的過程把每條read的精確比對(duì)信息壓縮成每個(gè)堿基位置上的覆蓋深度或信號(hào)強(qiáng)度。這個(gè)過程中涉及幾個(gè)關(guān)鍵決策用什么工具做轉(zhuǎn)換、bin size設(shè)多大、要不要做歸一化、怎么處理PCR重復(fù)和低質(zhì)量reads。UCSC工具鏈提供了一套完整的解決方案核心工具包括bamCoverage來自deepTools雖然不是UCSC官方但常與UCSC工具鏈配合使用、bedGraphToBigWigUCSC官方工具以及bamToBed等輔助工具。整個(gè)流程可以概括為BAM文件 → 過濾和質(zhì)控 → 生成bedGraph中間文件 → 排序和索引 → 轉(zhuǎn)換為bigWig。每一步都有坑每一步的參數(shù)選擇都會(huì)影響最終可視化效果。我見過太多人直接拿原始BAM轉(zhuǎn)bigWig結(jié)果圖上一堆假信號(hào)或者bin size設(shè)得太大導(dǎo)致峰形完全失真。下面我會(huì)把整個(gè)流程拆開把每個(gè)環(huán)節(jié)的細(xì)節(jié)和踩過的坑都講清楚。2. 工具鏈選型與核心原理拆解2.1 為什么是UCSC工具鏈而不是其他方案市面上做BAM到bigWig轉(zhuǎn)換的工具不少常見的有deepTools的bamCoverage、UCSC的bedGraphToBigWig、MACS2的bdgcmp、還有bam2wig等。為什么我推薦UCSC工具鏈作為核心原因有三第一bedGraphToBigWig是UCSC官方維護(hù)的格式兼容性最好。bigWig格式本身就是UCSC定義的用官方工具轉(zhuǎn)換出來的文件在UCSC Genome Browser上加載不會(huì)出現(xiàn)任何格式兼容問題。我試過用某些第三方工具生成的bigWig在IGV里能打開但傳到UCSC Browser上就報(bào)錯(cuò)排查半天發(fā)現(xiàn)是索引格式有細(xì)微差異。第二UCSC工具鏈的輸入輸出非常明確。bedGraphToBigWig只接受排序好的bedGraph文件和chromosome size文件輸出就是標(biāo)準(zhǔn)的bigWig。這種單一職責(zé)的設(shè)計(jì)讓流程非?？煽爻隽藛栴}容易定位。第三性能足夠好。bedGraphToBigWig底層用C實(shí)現(xiàn)處理幾GB的bedGraph文件只需要幾分鐘內(nèi)存占用也低。相比之下某些Python實(shí)現(xiàn)的工具在處理大文件時(shí)內(nèi)存直接爆掉。當(dāng)然deepTools的bamCoverage也是很好的選擇它可以直接從BAM生成bigWig省去了中間步驟。但它的缺點(diǎn)是參數(shù)太多新手容易設(shè)錯(cuò)而且生成的bigWig在UCSC Browser上偶爾會(huì)有兼容性警告。所以我建議的流程是用bamCoverage或samtools生成bedGraph然后用UCSC的bedGraphToBigWig做最終轉(zhuǎn)換。這樣既利用了deepTools的靈活性又保證了格式的規(guī)范性。2.2 bedGraph中間格式的關(guān)鍵作用bedGraph是一種純文本格式每行四個(gè)字段染色體、起始位置、終止位置、信號(hào)值。比如chr1 1000 1001 5 chr1 1001 1002 8 chr1 1002 1003 12這種格式的好處是直觀、可讀、可編輯。你可以用head命令直接查看內(nèi)容用awk做過濾和計(jì)算用sort排序。在調(diào)試階段這種透明性非常重要。我經(jīng)常遇到的情況是轉(zhuǎn)換出來的bigWig信號(hào)不對(duì)這時(shí)候如果中間有bedGraph文件我可以直接檢查某個(gè)區(qū)域的原始信號(hào)值判斷是BAM過濾的問題還是轉(zhuǎn)換參數(shù)的問題。bedGraph的缺點(diǎn)是文件體積大。一個(gè)30x覆蓋度的全基因組bedGraph文件可能有幾十GB因?yàn)槊總€(gè)堿基位置都有一行記錄。所以bedGraph只是中間產(chǎn)物最終一定要轉(zhuǎn)成bigWig。但正是這個(gè)大文件給了你最后檢查數(shù)據(jù)質(zhì)量的機(jī)會(huì)。2.3 bin size的選擇邏輯bin size是bigWig轉(zhuǎn)換中最關(guān)鍵的參數(shù)之一。它決定了每個(gè)數(shù)據(jù)點(diǎn)代表多少個(gè)堿基。bin size1意味著每個(gè)堿基一個(gè)信號(hào)值分辨率最高但文件最大bin size10意味著每10個(gè)堿基合并成一個(gè)信號(hào)值文件小但分辨率降低。怎么選取決于你的應(yīng)用場(chǎng)景。如果你要看轉(zhuǎn)錄因子結(jié)合位點(diǎn)這種窄峰通常幾十到幾百bp寬bin size必須設(shè)小建議1-5。如果你要看組蛋白修飾的寬峰可能幾kb寬bin size可以設(shè)10-50。如果是看全基因組覆蓋度分布bin size設(shè)100甚至1000都沒問題。這里有個(gè)經(jīng)驗(yàn)公式bin size不要超過你預(yù)期最窄峰寬度的1/3。比如你預(yù)期最窄的峰是60bp那bin size最大設(shè)20。否則峰形會(huì)被平滑掉看起來像個(gè)小土包而不是尖峰。還有一個(gè)坑bin size設(shè)得太小會(huì)導(dǎo)致文件巨大且充滿噪聲。比如單堿基分辨率的bigWig在低覆蓋區(qū)域會(huì)出現(xiàn)大量0和1的交替看起來像條形碼。這時(shí)候適當(dāng)?shù)腷in size比如5-10可以平滑掉這種噪聲讓信號(hào)更連續(xù)。3. 從BAM到bedGraph的實(shí)操細(xì)節(jié)3.1 BAM文件的預(yù)處理與過濾拿到BAM文件后不要直接轉(zhuǎn)換。先做幾件事第一步檢查BAM文件是否排序和索引。用samtools quickcheck檢查文件完整性用samtools idxstats查看是否有索引。如果沒有索引先samtools index input.bam。未排序的BAM文件無法直接用于大多數(shù)轉(zhuǎn)換工具。第二步?jīng)Q定是否去除PCR重復(fù)。如果是全基因組測(cè)序或ATAC-seqPCR重復(fù)通常要去掉否則重復(fù)區(qū)域的信號(hào)會(huì)被高估。用samtools markdup或picard MarkDuplicates標(biāo)記重復(fù)然后用samtools view -F 0x400過濾掉。但如果是RNA-seqPCR重復(fù)可能代表高表達(dá)轉(zhuǎn)錄本不建議盲目去除。我一般會(huì)先看picard CollectDuplicateMetrics的結(jié)果如果重復(fù)率超過30%再考慮去除。第三步過濾低質(zhì)量reads。常用參數(shù)是-q 20比對(duì)質(zhì)量≥20和-F 0x904過濾掉未比對(duì)、次要比對(duì)和PCR重復(fù)。對(duì)于ATAC-seq還要考慮過濾線粒體reads用samtools view -F 0x904 input.bam | grep -v chrM。第四步?jīng)Q定是否提取特定區(qū)域。如果你只關(guān)心某些基因區(qū)域可以用bedtools intersect或samtools view -L regions.bed提取目標(biāo)區(qū)域的reads這樣后續(xù)轉(zhuǎn)換會(huì)快很多。注意過濾步驟的順序很重要。先過濾質(zhì)量再去除重復(fù)比反過來更高效因?yàn)檫^濾掉低質(zhì)量reads后重復(fù)標(biāo)記的計(jì)算量會(huì)小很多。3.2 用deepTools bamCoverage生成bedGraphbamCoverage是deepTools中最常用的工具基本命令如下bamCoverage -b input.bam -o output.bedGraph \ --binSize 10 \ --normalizeUsing RPKM \ --ignoreDuplicates \ --minMappingQuality 20 \ --extendReads 200 \ --outFileFormat bedgraph逐參數(shù)解釋--binSize 10每10bp一個(gè)bin。根據(jù)前面說的原則如果你看的是窄峰改成5或1。--normalizeUsing RPKM歸一化方法。RPKM適合比較不同樣本間的表達(dá)量但如果你只是看覆蓋度分布用CPMcounts per million更簡(jiǎn)單。還有BPMbins per million和RPGCreads per genomic content選擇取決于你的下游分析目的。--ignoreDuplicates忽略PCR重復(fù)。如果前面已經(jīng)用samtools去除了這里可以不加。--minMappingQuality 20最低比對(duì)質(zhì)量。對(duì)于unique比對(duì)20是個(gè)安全值對(duì)于允許multi-mapping的分析可以設(shè)0。--extendReads 200將reads延伸到200bp。這個(gè)參數(shù)對(duì)ATAC-seq和ChIP-seq很重要因?yàn)閷?shí)際信號(hào)來自片段而非read本身。對(duì)于RNA-seq通常不需要延伸或者延伸到平均片段長度。--outFileFormat bedgraph輸出bedGraph格式。這里有個(gè)容易忽略的點(diǎn)--extendReads的值怎么定對(duì)于ATAC-seq可以用bamPEFragmentSize計(jì)算實(shí)際片段長度分布取中位數(shù)。對(duì)于ChIP-seq通常取150-300。如果設(shè)得太大信號(hào)會(huì)過度平滑設(shè)得太小信號(hào)會(huì)碎片化。3.3 用samtools和bedtools手動(dòng)生成bedGraph如果你不想依賴deepTools也可以用samtools和bedtools手動(dòng)生成bedGraph。流程如下# 1. 將BAM轉(zhuǎn)換為BED格式 samtools view -b -q 20 -F 0x904 input.bam | \ bedtools bamtobed -i stdin reads.bed # 2. 生成覆蓋度bedGraph bedtools genomecov -i reads.bed -g chrom.sizes -bg coverage.bedGraph # 3. 如果需要?dú)w一化用awk計(jì)算 awk BEGIN{scale1000000/TotalReads} {print $1\t$2\t$3\t$4*scale} \ coverage.bedGraph normalized.bedGraph這種方法的優(yōu)點(diǎn)是每一步都透明可控缺點(diǎn)是速度比bamCoverage慢尤其是bedtools genomecov在處理大BAM文件時(shí)。我實(shí)測(cè)過一個(gè)30GB的BAM文件bamCoverage大約15分鐘完成而bedtools genomecov需要40分鐘以上。提示bedtools genomecov的-bg參數(shù)輸出的是bedGraph格式-bga會(huì)輸出所有位置包括0覆蓋的區(qū)域。如果你要做全基因組可視化用-bga如果只關(guān)心有信號(hào)的位置用-bg可以顯著減小文件。3.4 歸一化的選擇與計(jì)算歸一化是決定bigWig可比性的關(guān)鍵。假設(shè)你有兩個(gè)樣本一個(gè)測(cè)了20M reads一個(gè)測(cè)了40M reads如果不歸一化第二個(gè)樣本的所有信號(hào)都是第一個(gè)的兩倍看起來好像所有區(qū)域都富集了這顯然是錯(cuò)的。常用的歸一化方法方法全稱適用場(chǎng)景計(jì)算方式RPKMReads Per Kilobase per MillionRNA-seq表達(dá)量比較reads數(shù)/(基因長度×總reads數(shù)/1e6)CPMCounts Per Million通用覆蓋度比較reads數(shù)/(總reads數(shù)/1e6)BPMBins Per MillionbigWig專用bin內(nèi)reads數(shù)/(總bins數(shù)/1e6)RPGCReads Per Genomic Content全基因組覆蓋度reads數(shù)/(總reads數(shù)/基因組長度)對(duì)于bigWig可視化我通常推薦CPM或BPM。CPM簡(jiǎn)單直接BPM是deepTools的默認(rèn)推薦。RPGC適合比較不同基因組的樣本但計(jì)算稍復(fù)雜。如果你用bamCoverage直接用--normalizeUsing參數(shù)即可。如果手動(dòng)生成bedGraph需要自己計(jì)算scale factor# 獲取總reads數(shù) TotalReads$(samtools view -c -q 20 -F 0x904 input.bam) # 計(jì)算scale factor以CPM為例 ScaleFactor$(echo 1000000 / $TotalReads | bc -l) # 應(yīng)用歸一化 awk -v scale$ScaleFactor {print $1\t$2\t$3\t$4*scale} \ coverage.bedGraph normalized.bedGraph注意歸一化后的信號(hào)值可能是小數(shù)bedGraphToBigWig要求信號(hào)值為整數(shù)或浮點(diǎn)數(shù)。如果bedGraph中有科學(xué)計(jì)數(shù)法如1.2e-5需要先轉(zhuǎn)換為普通小數(shù)格式否則會(huì)報(bào)錯(cuò)。4. bedGraph到bigWig的轉(zhuǎn)換與優(yōu)化4.1 準(zhǔn)備chromosome size文件bedGraphToBigWig需要一個(gè)chromosome size文件格式為兩列染色體名和長度。這個(gè)文件必須與BAM文件的參考基因組一致。獲取方法# 從BAM文件頭獲取 samtools view -H input.bam | grep SQ | \ sed s/SQ\tSN://;s/\tLN:// chrom.sizes # 或者從參考基因組fasta獲取 samtools faidx reference.fa cut -f1,2 reference.fa.fai chrom.sizes注意chrom.sizes中的染色體名必須與bedGraph中的完全一致。如果BAM里是chr1而chrom.sizes里是1轉(zhuǎn)換會(huì)失敗或產(chǎn)生空文件。我遇到過好幾次這個(gè)問題排查半天才發(fā)現(xiàn)是命名不一致。4.2 排序bedGraph文件bedGraphToBigWig要求bedGraph按染色體和起始位置排序。如果bedGraph是bamCoverage生成的通常已經(jīng)排序好了。但如果是手動(dòng)生成的需要先排序# 按染色體和起始位置排序 sort -k1,1 -k2,2n input.bedGraph sorted.bedGraph這里有個(gè)坑sort命令的默認(rèn)排序是字典序?qū)τ谌旧w名如chr10和chr2字典序會(huì)排成chr10在chr2前面但bigWig要求按染色體在chrom.sizes中的順序排列。所以更安全的方法是用bedtools sortbedtools sort -i input.bedGraph -g chrom.sizes sorted.bedGraphbedtools sort會(huì)根據(jù)chrom.sizes中的順序排序確保與bigWig的要求一致。4.3 執(zhí)行轉(zhuǎn)換與驗(yàn)證轉(zhuǎn)換命令很簡(jiǎn)單bedGraphToBigWig sorted.bedGraph chrom.sizes output.bigWig但轉(zhuǎn)換完成后一定要驗(yàn)證。驗(yàn)證方法# 檢查bigWig基本信息 bigWigInfo output.bigWig # 提取某個(gè)區(qū)域的信號(hào)值 bigWigToBedGraph -chromchr1 -start1000000 -end1001000 \ output.bigWig /dev/stdoutbigWigInfo會(huì)輸出染色體數(shù)量、總數(shù)據(jù)點(diǎn)數(shù)、最小最大值等信息。如果輸出顯示chromCount: 0說明轉(zhuǎn)換失敗通常是chrom.sizes不匹配或bedGraph未排序。提示如果bedGraph文件很大超過10GBbedGraphToBigWig可能會(huì)因?yàn)閮?nèi)存不足而失敗。這時(shí)候可以先用split命令按染色體拆分bedGraph分別轉(zhuǎn)換后再用bigWigMerge合并。雖然麻煩但能解決問題。4.4 在基因組瀏覽器中加載與調(diào)整生成的bigWig可以直接拖入IGV或上傳到UCSC Genome Browser。在IGV中你可以調(diào)整顯示模式Bar chart適合看覆蓋度Heatmap適合看多個(gè)樣本的比較Line plot適合看連續(xù)信號(hào)。在UCSC Browser中通過add custom track上傳bigWig可以設(shè)置viewLimits控制顏色范圍。比如設(shè)置viewLimits0:50超過50的信號(hào)都顯示為最深顏色這樣可以避免個(gè)別極高信號(hào)點(diǎn)導(dǎo)致整體顏色過淺。注意UCSC Browser對(duì)上傳文件大小有限制通常不超過500MB。如果bigWig太大可以先用bigWigToBedGraph提取目標(biāo)區(qū)域再轉(zhuǎn)換回bigWig。5. 常見問題與排查技巧實(shí)錄5.1 轉(zhuǎn)換失敗與報(bào)錯(cuò)處理問題一bedGraphToBigWig報(bào)錯(cuò)Invalid coordinate原因bedGraph中有起始位置大于終止位置的行或者有負(fù)值。排查方法awk $2$3 || $20 || $30 input.bedGraph | head解決方法用awk過濾掉這些行或者檢查生成bedGraph的工具是否有bug。問題二bigWig文件為空或只有部分染色體原因chrom.sizes與bedGraph的染色體命名不一致或者bedGraph未排序。排查方法# 檢查bedGraph中的染色體名 cut -f1 input.bedGraph | sort -u # 檢查chrom.sizes中的染色體名 cut -f1 chrom.sizes | sort -u解決方法統(tǒng)一命名確保chrom.sizes包含bedGraph中所有染色體。問題三轉(zhuǎn)換后的bigWig在IGV中顯示異常原因可能是bin size太小導(dǎo)致噪聲或者歸一化參數(shù)設(shè)錯(cuò)。排查方法用bigWigToBedGraph提取一段區(qū)域檢查信號(hào)值是否合理。5.2 信號(hào)異常的診斷思路信號(hào)全為0檢查BAM文件是否有比對(duì)結(jié)果samtools flagstat檢查過濾參數(shù)是否過嚴(yán)比如-q 30可能過濾掉大部分reads檢查chrom.sizes是否與BAM參考基因組一致。信號(hào)出現(xiàn)周期性波動(dòng)這是bin size與read長度不匹配的典型表現(xiàn)。比如read長度150bpbin size設(shè)100就會(huì)出現(xiàn)每100bp一個(gè)峰。解決方法bin size設(shè)為read長度的約數(shù)或者用--extendReads平滑。信號(hào)在基因間區(qū)異常高可能是背景噪聲或比對(duì)錯(cuò)誤。檢查是否過濾了multi-mapping reads是否去除了重復(fù)序列區(qū)域。可以用blacklist文件過濾掉已知的假信號(hào)區(qū)域。不同樣本間信號(hào)不可比檢查歸一化方法是否一致測(cè)序深度是否相近。如果深度差異超過5倍建議用RPGC或先下采樣到相同深度。5.3 性能優(yōu)化與批量處理如果你有幾十個(gè)BAM文件要轉(zhuǎn)換手動(dòng)一個(gè)個(gè)跑效率太低。寫個(gè)循環(huán)腳本for bam in *.bam; do sample$(basename $bam .bam) bamCoverage -b $bam -o ${sample}.bedGraph \ --binSize 10 --normalizeUsing CPM \ --ignoreDuplicates --minMappingQuality 20 \ --extendReads 200 --outFileFormat bedgraph \ --numberOfProcessors 8 bedtools sort -i ${sample}.bedGraph -g chrom.sizes ${sample}.sorted.bedGraph bedGraphToBigWig ${sample}.sorted.bedGraph chrom.sizes ${sample}.bigWig rm ${sample}.bedGraph ${sample}.sorted.bedGraph done--numberOfProcessors 8可以顯著加速bamCoverage但bedGraphToBigWig是單線程的無法并行。如果bedGraph文件很大可以考慮按染色體拆分后并行轉(zhuǎn)換。提示批量處理時(shí)建議先用一個(gè)樣本測(cè)試整個(gè)流程確認(rèn)參數(shù)無誤后再跑全部。我吃過虧跑了20個(gè)樣本才發(fā)現(xiàn)歸一化參數(shù)設(shè)錯(cuò)全部重來。5.4 常見問題速查表問題現(xiàn)象可能原因排查方法解決方案轉(zhuǎn)換報(bào)錯(cuò)Invalid coordinatebedGraph有負(fù)值或起止顛倒awk檢查過濾異常行bigWig為空染色體命名不一致對(duì)比cut -f1結(jié)果統(tǒng)一命名IGV中信號(hào)全為0過濾參數(shù)過嚴(yán)samtools flagstat放寬過濾信號(hào)周期性波動(dòng)bin size與read長度不匹配檢查read長度調(diào)整bin size樣本間不可比歸一化不一致檢查歸一化參數(shù)統(tǒng)一歸一化轉(zhuǎn)換速度慢bedGraph未排序檢查排序用bedtools sort內(nèi)存不足bedGraph過大檢查文件大小按染色體拆分6. 從bigWig到可視化呈現(xiàn)的進(jìn)階技巧6.1 多樣本bigWig的合并與比較如果你有多個(gè)樣本的bigWig想在同一個(gè)圖中比較有兩種方法方法一用bigWigMerge合并。將多個(gè)bigWig合并成一個(gè)信號(hào)值相加或取平均。適合展示總信號(hào)強(qiáng)度。bigWigMerge sample1.bigWig sample2.bigWig merged.bedGraph bedGraphToBigWig merged.bedGraph chrom.sizes merged.bigWig方法二在IGV中疊加顯示。將多個(gè)bigWig加載到同一個(gè)IGV session設(shè)置不同的顏色和透明度可以直觀比較樣本間的差異。這種方法不需要合并文件更靈活。6.2 差異信號(hào)的bigWig生成如果你想展示兩個(gè)樣本間的差異信號(hào)比如處理組vs對(duì)照組可以先生成差異bedGraph再轉(zhuǎn)bigWig# 用bigWigCompare或手動(dòng)計(jì)算差異 bigWigCompare sample1.bigWig sample2.bigWig diff.bedGraph # 或者用awk計(jì)算log2比值 paste sample1.bedGraph sample2.bedGraph | \ awk {if($40 $80) print $1\t$2\t$3\tlog($4/$8)/log(2)} \ log2ratio.bedGraph差異bigWig在IGV中可以用Heatmap模式顯示紅色代表上調(diào)藍(lán)色代表下調(diào)非常直觀。6.3 在UCSC Genome Browser中創(chuàng)建自定義trackUCSC Browser支持通過URL加載bigWig格式如下track typebigWig nameMy Track bigDataUrlhttps://example.com/output.bigWig你可以把多個(gè)track寫在一個(gè)trackDb.txt文件中上傳到自定義track hub。這樣別人訪問你的hub URL就能看到所有track不需要手動(dòng)上傳文件。track hub的配置稍微復(fù)雜但一次配置好后非常方便。注意track hub需要HTTPS鏈接且服務(wù)器要支持byte-range請(qǐng)求。如果用自己的服務(wù)器確保配置了Accept-Ranges: bytes響應(yīng)頭。6.4 可視化效果的調(diào)優(yōu)經(jīng)驗(yàn)同樣的bigWig不同的顯示參數(shù)效果可能天差地別。幾個(gè)調(diào)優(yōu)經(jīng)驗(yàn)顏色范圍不要用默認(rèn)的自動(dòng)范圍。如果有個(gè)別極高信號(hào)點(diǎn)自動(dòng)范圍會(huì)把大部分區(qū)域壓成淺色。手動(dòng)設(shè)置viewLimits比如0:100讓主要信號(hào)區(qū)域有足夠的對(duì)比度。平滑窗口IGV支持Windowing Function可以設(shè)為mean或max。對(duì)于覆蓋度數(shù)據(jù)mean更平滑max更突出峰。我通常用mean看整體趨勢(shì)用max找精確峰位置。多track對(duì)齊在IGV中把多個(gè)track的Data Range設(shè)為相同值這樣顏色深淺可以直接比較。如果每個(gè)track自動(dòng)縮放顏色就沒有可比性了。導(dǎo)出高分辨率圖片IGV的Save Image功能可以導(dǎo)出PNG或SVG。SVG是矢量圖放大不模糊適合放到論文或報(bào)告里。導(dǎo)出時(shí)注意設(shè)置合適的分辨率通常300dpi足夠。6.5 實(shí)際項(xiàng)目中的流程整合在一個(gè)典型的測(cè)序分析項(xiàng)目中BAM到bigWig的轉(zhuǎn)換通常不是孤立的步驟而是整個(gè)流程的一部分。我通常把它整合到Snakemake或Nextflow流程中確?？芍貜?fù)性。一個(gè)簡(jiǎn)化的Snakemake規(guī)則rule bam_to_bigwig: input: bam aligned/{sample}.bam, chrom_sizes reference/chrom.sizes output: bigwig bigwig/{sample}.bigWig params: bin_size 10, normalize CPM shell: bamCoverage -b {input.bam} -o temp/{wildcards.sample}.bedGraph \ --binSize {params.bin_size} --normalizeUsing {params.normalize} \ --ignoreDuplicates --minMappingQuality 20 --extendReads 200 \ --outFileFormat bedgraph --numberOfProcessors 8 bedtools sort -i temp/{wildcards.sample}.bedGraph -g {input.chrom_sizes} \ temp/{wildcards.sample}.sorted.bedGraph bedGraphToBigWig temp/{wildcards.sample}.sorted.bedGraph \ {input.chrom_sizes} {output.bigwig} rm temp/{wildcards.sample}.bedGraph temp/{wildcards.sample}.sorted.bedGraph 這樣每次有新樣本只需要把BAM放到指定目錄運(yùn)行Snakemake就能自動(dòng)生成bigWig。流程化的好處是參數(shù)統(tǒng)一、結(jié)果可重復(fù)、出錯(cuò)容易追溯。我個(gè)人在實(shí)際操作中的體會(huì)是BAM到bigWig的轉(zhuǎn)換看似簡(jiǎn)單但細(xì)節(jié)決定成敗。bin size、歸一化方法、過濾參數(shù)這三個(gè)東西每次做新項(xiàng)目都要根據(jù)數(shù)據(jù)特點(diǎn)重新考慮不能無腦套用之前的參數(shù)。尤其是做多個(gè)樣本比較時(shí)歸一化方法必須統(tǒng)一否則圖做得再漂亮也是錯(cuò)的。還有一個(gè)容易被忽略的點(diǎn)是chrom.sizes文件我至少遇到過五次因?yàn)槿旧w命名不一致導(dǎo)致轉(zhuǎn)換失敗的情況現(xiàn)在每次都會(huì)先檢查一遍再跑流程。