[ Deeptools ] DBS(differential binding sites) 를 Deeptools의 ComputeMatrix와 plotHeatmap을 이용해서 시각화하기

DBA (Differential Binding sites) visualization

Chipseq 이나 Cut&Run, Cut&Tag 기술을 이용해서, target antibody가 붙는 sequence 정보를 얻을 수 있고 어느 지역에 더 많이 붙는지 DBA(differential binding sites analysis) 분석을 할 수 있다. 해당 글에서는 DBA(differential binding analysis) 로 확인된 지역을 Deeptools 의 computeMatrix 와 plotHeapmap을 이용하여 시각화 할 것이다.

 

 

Step 1. bigWigMerge (ucsc-bigwigmerge v377)

multiple bigWigs → A single Bedgraph file

bigWigMerge -max $bigwigs ${prefix}_merged.bg
# -max:merged value is maximum from input files rather than sum

LC_COLLATE=C sort -k1,1 -k2,2n ${prefix}_merged.bg > ${prefix}_merged.sort.bedgraph
# LC_COLLATE=C sort: 문자열 정렬 순서 결정을 보여준다. 해당 명령어는 column 1(-k1,1)을 첫번째,(-k2, 2n)을 두번째 기준으로 정렬한다. (2n 에서 n 은 numeric sort.)

 

Step 2. bedGraphToBigWig (ucsc-bedgraphtobigwig v377)

Step 1 에서 생성된 bedgraph 파일을 다시 bigwigs 파일로 변형 시켜 준다.

bedGraphToBigWig $bedgraph $<chrom_sizes> ${prefix}.bigWig

* BigWig 는 indexed binary format으로, Genome Browser server 나 IGV 같은 시각 Toll 에서 특정 부분을 시각화 할 수 있도록 만들어진 파일이다.

# bigWig file 
# bamCoverage 를 이용하여 BAM -> BigWig (Binary File)로 만들어진다. 
# 이때, mapping score/coverage = number of reads that overlap each genomic position 으로 계산되고, 해당 score은 CPM, RPM, FPKM 등과 같은 방법으로 normalized 할 수 있다.

# bedgraph file 예시 
scaffold_1      100     150     0.226466
scaffold_1      150     200     0.198157
scaffold_1      200     300     0.141541
scaffold_1      300     350     0.0566164
scaffold_1      350     600     0.0283082
scaffold_1      600     650     0.0566164
scaffold_1      650     700     0.0283082

 

 

Step 3. compute the count Matrix

Score per genome regions를 계산하여 computeMatrix (genomic regions(rows) x scores (column)) 를 생성한다. computeMatrix 는 2개의 모드로 사용이 가능하며 이후, matrix를 이용해서 plotHeatmap, plotProfiles 를 이용해서 시각화 가능하다. 

• reference-point 하나의 포인트를 기반으로 (유전자 시작이라던가 끝이라던가) signal distribution을 보여준다.
• scale-regions 모든 regions을 같은 사이즈로 scaling을 하여$서, 지역 간의 signal difference를 보여준다.

(1) reference-point 방법

# < computeMatrix 실행 > nextflow이용해서 변수설정이 이렇게 생겼지만 다르게해도 당연히 된다.
params {
	beforelen=3000
	afterlen=3000
}

params {
  bed=diffbound_regions.bed # differentally binding regions (both up/down regulated) 
}

computeMatrix \\
	-reference-point --referencePoint center \\
	-a ${params.beforelen} -b ${params.afterlen} \\
	--missingDataAsZero --skipZeros \\
	--averageTypeBins mean \\
	--sortRegions descend \\
	--regionsFileName $bed \\
    --scoreFileName $bigwig \\
    --outFileName ${prefix}.computeMatrix.mat.gz \\
    --outFileNameMatrix ${prefix}.computeMatrix.vals.mat.tab \\
    --outFileSortedRegions ${prefix}.sortedRegions.bed \\
    --numberOfProcessors $task.cpus

# reference-point [TSS|TES|center]: position within a BED region
# -a(--afterRegionsStartLength): distance upstream of the start site (TSS) of the regions
# -b(--beforeRegionsStartLength): distance downstream of the end site of the given regions
# --regionsFileName: regions to plot
# --scoreFileName: bigWig file(s) containing the scores to be plotted
# --outFileName: 나중에 plotHeatmap이나 plotProfile을 만들 때 필요한 gzipped matrix 저장위치

# < Ouput > 
# .computeMatrix.vals.mat.tab
#diffbound_regions.bed:1637	genes:1748
#downstream:3000	upstream:3000	body:0	bin size:10	unscaled 5 prime:0	unscaled 3 prime:0
diffbound_regions.bed:1637	genes:1748	A1_H3K4_IRI__CPM	A1_H3K4_IRI__CPM	A1_H3K4_IRI__CPM	A1_H3K4_IRI__CPM	A1_H3K4_IRI__CPM	A1_H3K4_IRI__CPM	A1_H3K4_I
0.1241	0.1489	0.1551	0.1551	0.1551	0.1551	0.1799	0.1861	0.1861	0.1861	0.1861	0.2357	0.2481	0.2481	0.2481	0.2481	0.2977	0.3102	0.3102	0.3102	0.3102	0.3598	0.3722	0.3722	0.3722	0.3722	0.4466	

# .sortedRegions.bed
#chrom	start	end	name	score	strand	thickStart	thickEnd	itemRGB	blockCount	blockSizes	blockStart	deepTools_group
scaffold_6	119064418	119069759	._r848	0.0	+	119064418	119069759	0	1	5341	119064417	spiny_diffbound_regions.bed
scaffold_15	82612901	82618242	._r1075	0.0	+	82612901	82618242	0	1	5341	82612893	spiny_diffbound_regions.bed
scaffold_8	107305850	107311191	._r568	0.0	+	107305850	107311191	0	1	5341	107305849	spiny_diffbound_regions.bed
scaffold_10	86210219	86215560	._r729	0.0	+	86210219	86215560	0	1	5341	86210211	spiny_diffbound_regions.bed
scaffold_13	10016502	10021843	._r95	0.0	+	10016502	10021843	0	1	5341	10016501	spiny_diffbound_regions.bed
scaffold_14	63047831	63053172	._r663	0.0	+	63047831	63053172	0	1	5341	63047825	spiny_diffbound_regions.bed
scaffold_14	14213455	14218796	._r676	0.0	+	14213455	14218796	0	1	5341	14213454	spiny_diffbound_regions.bed
## blockCount = 2 이면 원래 bed file의 region 이 2 개로 나눠졌다는 뜻이다.

(2) Scale-regions 방법

params {
  bed=diffbound_transcript_regions.bed # differentally binding regions that is overlapping with known genome regions 
}

computeMatrix \\ 
	scale-regions -a ${params.beforelen} -b ${params.afterlen} \\
    -m ${params.gene_body_len} \\
    --skipZeros --missingDataAsZero \\
	--averageTypeBins mean --sortRegions descend \\
    --regionsFileName $bed \\
    --scoreFileName $bigwig \\
    --outFileName ${prefix}.computeMatrix.mat.gz \\
    --outFileNameMatrix ${prefix}.computeMatrix.vals.mat.tab \\
    --outFileSortedRegions ${prefix}.sortedRegions.bed \\
    --numberOfProcessors $task.cpus

 

 

Step. 4 plotHeatmap 을 이용하여 시각화