[11] 7 数据可视化 & 主要涉及 bam 文件转变

2019-07-07 本文已影响39人热衷组培的二货潜

7 数据可视化

在整个分析过程中生成的数据的可视化对于获得数据质量的印象、理解结果以及将数据与其他信息集成是至关重要的。此外，运行 Focus 到感兴趣的点上。本章节，我们将生成能够在基因组浏览器中以reads 或 Peak 级别可视化数据的文件。

7.1 Reads 密度

Reads 的可视化可以在基因组比对之后进行（见章节 4 ）。包含 reads 比对信息的 BAM 文件可以直接导入到基因组浏览器中（ reads 比对到正向或反向链的情况明）（图 7.1 ）。然而，在具有高reads 密度的区域，例如 Peak ，其中数百个 reads 比对上面，这样的可视化不能很好地 scale 。因此，这种方法仅在单个 reads 级别调查信息时有用（比如：精确位置、错配）。
可以生成基因组覆盖文件来总结与基因组中的每个位置重叠的 reads ，而不是将每个单个 reads 可视化（图 7.1 ）。由于测序的 reads 仅代表由感兴趣的 TF 结合的共同纯化的DNA片段的末端（例如50bp ），每个 reads 在 3‘ 方向延伸到 200bp。200bp 表示在 ChIP-seq 流程中实验选择的片段大小的粗略估计（见章节 1.1 ）。这使信号平滑，并确保 Peak 的 Summit 和出现在结合位点之上。为了能够比较不同的样本，使用 Scale factor 将每个位置的 reads 数均一化为一百万条比对上的 reads （ 1,000,000 / 比对上的 reads 数 ）。reads 的基因组覆盖率以bedGraph格式（ chromosome / start / end / reads_density ）存储，并且可以压缩为 bigWig 格式。

图 7.1 在基因组浏览器中的数据可视化。

从样本 NRF1_ChIP_WT_1 集中于 mm10 基因组上的位置chr2:153652399的 score 为3的 Peak 的可视化。为由UCSC基因组浏览器生成显示比对到正向（深灰色）或反向链（浅灰色）的reads 密度、峰值区域和原始 reads 。
image

# Create a directory for the density tracks
mkdir tracks

sample=NRF1_CHIP_WT_1

# Chromosome sizes (to calculate read densities at all positions)
# To be downloaded from UCSC (genome.ucsc.edu/), Downloads,

Genome Data
# File is sorted according to chromosome and only conventional chromosomes are kept

wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/mm10.chrom.sizes 

sort -k1,1 mm10.chrom.sizes | grep -v chrUn |
grep -v random > mm10.chrom.sizes.tmp

mv mm10.chrom.sizes.tmp mm10.chrom.sizes

# Scale to normalise the number of reads at each position to 1,000,000 mapped reads in the library

# Count total number of mapped reads if the index BAI file is available (˜1 min)

samtools idxstats reads/${sample}.bam | 
awk '{if($1!="*"){total=total+$3}}END{print 1000000/total}'

# Otherwise count total number of mapped reads from the BAM file (˜ 2 min)
bamToBed -i reads/${sample}.bam | wc -l | awk '{print 1000000/total}'

scale=0.0322351

# One-line command
scale=$(samtools idxstats reads/${sample}.bam |\
awk '{if($1!="*"){total=total+$3}}END{print 1000000/total}')


# Extend reads checking they do not extend over chromosome limits

extend=200

# Extend reads from the forward strand (since files are sorted by starts, there is no need to sort again)

bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" \
-v E=$extend \
-v file=mm10.chrom.sizes
'BEGIN{while(getline<file){S[$1]=$2}}($6=="+"){if($2+E>S[$1]){print $1,$2,S[$1]}\
else{print $1,$2,$2+E}}’ > tracks/${sample}_forward.bed

# Extend reads from the reverse strand (since ends are not sorted if reads do not have the same length e.g. after trimming, we need to sort the output again by chromosome / start)

bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" -v E=$extend\
'($6=="-"){if($3-E<1){print $1,"1",$3}\
else{print $1,$3-E,$3}}' |\
sort -k1,1 -k2,2n > tracks/${sample}_reverse.bed

# Merge sorted file (sort -m much faster than just sort since input files are already sorted)

sort -m -k1,1 -k2,2n tracks/${sample}_forward.bed
tracks/${sample}_reverse.bed > tracks/${sample}.bed


# Generate bedGraph file from BED (read coverage at each non-zero position)

genomeCoverageBed -bg -i tracks/${sample}.bed -g
mm10.chrom.sizes -scale $scale > tracks/${sample}.bg

# Remove intermediate files
rm tracks/${sample}_forward.bed tracks/${sample}_reverse.bed tracks/${sample}.bed

# One-line command (˜6 min)
extend=200
sort -m -k1,1 -k2,2n <\
(bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" 
-v E=$extend 
-v file= mm10.chrom.sizes
'BEGIN{while(getline<file){S[$1]=$2}}($6=="+"){if($2+E>S[$1]){print $1,$2,S[$1]}else{print $1,$2,$2+E}}') <\
(bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" -v E=$extend
'($6=="-"){if($3-E<1){print $1,"1",$3}else{print
$1,$3-E,$3}}' |\
sort -k1,1 -k2,2n) |\
genomeCoverageBed -bg -i stdin -g mm10.chrom.sizes -scale $scale > tracks/${sample}.bg


# Generate bigWig file from bedGraph (˜2 min)
bedGraphToBigWig tracks/${sample}.bg mm10.chrom.sizes tracks/${sample}.bw

# Remove intermediate file
rm tracks/${sample}.bg

7.2 Peak 区间

可以在 Peak calling 步骤之后执行 Peak 区域的可视化（ 章节 6 ）。可以直接使用与峰值区域对应的简单 BED 文件（ (chromosome / start / end ），也可以将其压缩为 bigBed 格式。
```
# Generate a simple BED file (chromosome / start / end)

cut -f1-3 peaks/${sample}_peaks_peakzilla.bed >
peaks/${sample}_peaks_peakzilla_short.bed

# Generate bigBed file from BED
bedToBigBed peaks/${sample}_peaks_peakzilla_short.bed
mm10.chrom.sizes tracks/${sample}_peaks_peakzilla.bb
rm peaks/${sample}_peaks_peakzilla_short.bed
```
注意：在基因组浏览器中探索最佳和最差排序的 Peak 允许我们评估选定的阈值，并检查这些峰是否在其他样本中共享。然而，重要的是要记住，这种可视化并不考虑数据的所有潜在特征。因此，不值得调整初始 Peak calling 参数以适应任何一个给定区域，因为它可能不适合其他区域。

7.3 基因组浏览器

加州大学圣克鲁斯分校（ UCSC ）基因组浏览器（ ref1 ）可以在线使用，也可以在本地安装。Broad 研究所开发的综合基因组查看器（ IGV ）（ ref2 ）是一种流行的替代方法，可以在本地运行。在没有访问带有 Web 服务器的专用服务器的情况下工作时，这一点特别有用。

一旦在 UCSC 中选择了参考基因组，就可以获得各种预加载的数据tracks（例如：基因注释信息）。

# BAM files (with .bai file in same directory)
# track type=bam name="NRF1_CHIP_WT_1_reads"
bigDataUrl=http://your_server/NRF1_CHIP_WT_1.bam

# bigWig files
# track type=bigWig name="NRF1_CHIP_WT_1_density"
bigDataUrl=http://your_server/NRF1_CHIP_WT_1.bw

# bigBed files
# track type=bigBed name="NRF1_CHIP_WT_1_peaks" bigDataUrl=
http://your_server/NRF1_CHIP_WT_1_peaks_peakzilla.bb

注意：与本书相关的数据，包括要加载到 UCSC 基因组浏览器中的 URL，可以在线获得。

关于基因组浏览器使用不做赘述。

7.1 代码解析

1 创造 track 文件夹
```
mkdir tracks
```

2 从 UCSC 获取染色体大小文件

wget
http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/mm10.chrom.sizes

sort -k1,1 mm10.chrom.sizes | grep -v chrUn |\
grep -v random > mm10.chrom.sizes.tmp

mv mm10.chrom.sizes.tmp mm10.chrom.sizes

3 将每个文库的 reads 数放缩到 100万 reads 级别

# 这一步之前需为 bam 文件建立index，生成 .bai 文件。
# 计算比对上的总 reads 数目
samtools idxstats reads/${sample}.bam |\
awk '{if($1!="*"){total=total+$3}}END{print 1000000/total}'

# 或者，统计从 BAM 文件中比对的 reads 总数
bamToBed -i reads/${sample}.bam | wc -l |\
awk '{print 1000000/total}'

# 这里可以得到
scale=0.0322351

将步骤 3 合并为一步

scale=$(samtools idxstats reads/${sample}.bam | awk '{if($1!="*"){total=total+$3}}END{print 1000000/total}')

4 Extend reads 检查它们不会超出染色体长度范围

extend=200
#从正向向链扩展 reads ( 因为文件是按 Start 列排序的，所以不需要再次排序)
bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" 
-v E=$extend 
-v file=mm10.chrom.sizes
'BEGIN{
  while(getline<file){
  S[$1]=$2
      }
  }
  ($6=="+"){if($2+E>S[$1]){
      print $1,$2,S[$1]
      }
  else{
      print $1,$2,$2+E
      }
  }' > tracks/${sample}_forward.bed
  
  
# 从反向链延伸读取( 因为如果 reads 不具有相同的长度，例如修剪后，就不会对末端进行排序，我们需要再次按 染色体 / Start 对输出进行排序)
bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" 
-v E=$extend 
'($6=="-"){if($3-E<1){
  print $1,"1",$3
  }
  else{
  print $1,$3-E,$3
      }
  }' |\
sort -k1,1 -k2,2n > tracks/${sample}_reverse.bed

5 合并正向和反向的 bed 文件（ sort-m 比仅仅排序快得多，因为输入文件已经排序）

sort -m -k1,1 -k2,2n\
tracks/${sample}_forward.bed \
tracks/${sample}_reverse.bed \
> tracks/${sample}.bed

6 将 BED 文件转变为 bedGraph 文件（每个非零位点的 reads 覆盖度）

genomeCoverageBed -bg -i tracks/${sample}.bed\
-g mm10.chrom.sizes \
-scale $scale > tracks/${sample}.bg

# 删除中间文件
rm\
tracks/${sample}_forward.bed\
tracks/${sample}_reverse.bed\
tracks/${sample}.bed

4 5 6 一行命令

extend=200

sort -m -k1,1 -k2,2n <\
(bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" 
-v E=$extend 
-v file= mm10.chrom.sizes
'BEGIN{while(getline<file){S[$1]=$2}}
  ($6=="+"){
  if($2+E>S[$1]){
      print $1,$2,S[$1]
      }
     else{
      print $1,$2,$2+E
      }
  }'
) <\
(bamToBed -i reads/${sample}.bam |\
awk -v OFS="\t" 
-v E=$extend
'($6=="-"){
  if($3-E<1){
  print $1,"1",$3
  }
  else{
  print
      $1,$3-E,$3
      }
    }' |\
sort -k1,1 -k2,2n
) |\
genomeCoverageBed -bg -i stdin\
-g mm10.chrom.sizes -scale $scale 
> tracks/${sample}.bg

将 bedGraph文件转变为 bigWig 文件

bedGraphToBigWig tracks/${sample}.bg mm10.chrom.sizes
tracks/${sample}.bw

# 删除中间文件
rm tracks/${sample}.bg

ref1：The human genome browser at UCSC.

ref2：Integrative genomics viewer