donglin.li

主页/Home 履历/CV 帖文/Posts 深色/Dark

🧬 · Hi-C笔记：从原始数据到接触矩阵

Posted on 2024/8/26

    本文为Hi-C数据处理的初学者笔记，包含AI生成的代码解释
    使用的操作系统为macOS(Apple Silicon)，但是应该适用于其他系统
    原教程：https://github.com/hms-dbmi/hic-data-analysis-bootcamp

    基础指令

    mkdir #创建目录（即文件夹）
cd #进入文件夹
cd .. #返回上一级
ls #列出文件夹内容
rm #删除文件
pwd #显示当前目录
head # 显示文件的开头部分
| #管道符号，将上一个指令的输出作为下一个指令的输入
\  #换行

docker exec #在Docker容器中执行指令

    前期准备

    下载Docker Desktop

    安装以下镜像

    docker pull duplexa/4dn-hic
docker pull higlass/higlass-docker

    4DN Hi-C

    HiGlass

    • 可视化工具

    将原始数据下载到以下路径：
    /Users/lidonglin/Documents/HiC_Exercise

    cd /Users/lidonglin/Documents/HiC_Exercise
mkdir data
cd data/
# input fastq files
wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/input_R1.fastq.gz
wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/input_R2.fastq.gz
gunzip input_R1.fastq.gz
gunzip input_R2.fastq.gz
# bwa genome index
wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/hg38.bwaIndex.tgz
tar -xzf hg38.bwaIndex.tgz
rm hg38.bwaIndex.tgz
# chromsizes
wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/hg38.mainonly.chrom.size
# prebaked output files
wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/prebaked.tgz
tar -xzf prebaked.tgz
rm prebaked.tgz
# move back a directory
cd ..

    如果没有wget可以使用homebrew安装，或者使用curl -o代替

    现在该目录中应该有这些文件

    • FASTQ 格式的 Hi-C 输入数据

      • input_R1.fastq
      • input_R2.fastq

      在文库构建之后，通常会使用Illumina测序平台对文库进行双端测序（paired-end sequencing），这意味着每个DNA片段的两端都会被测序，从而生成两个FASTQ文件。这两个文件中的每一条序列是配对的，代表了来自同一个DNA片段的两端

    • bwaIndex (hg38，人) - 适用于 BWA 的参考基因组格式

      • hg38.fasta.amb
      • hg38.fasta.ann
      • hg38.fasta.bwt
      • hg38.fasta.pac
      • hg38.fasta.sa

    • chrom size (hg38) - 列出染色体及其大小的文本文件

      • hg38.mainonly.chrom.size

    除此之外，其他的文件是教程作者准备的中间输出文件，可以暂时归档。

    如果要自己准备这些文件，以果蝇(dm6)的基因组为例

    在UCSC Genome Browser找到这个链接：https://hgdownload.soe.ucsc.edu/goldenPath/dm6/bigZips/

    可以下载dm6.fa.gz和dm6.chrom.sizes

    然后，按照下文的方法将目录挂载到容器，进入目录后执行以下命令

    gunzip dm6.fa.gz #解压
bwa index dm6.fa #创建索引

    稍等片刻即可

    将准备好的data目录挂载到容器

    docker run -it -v /Users/lidonglin/Documents/HiC\_Exercise/data:/d/:rw duplexa/4dn-hic bash

    显示目录内容，并进入d目录

    ls
ls /d
cd /d/

    检查文件

    head input_R1.fastq
head hg38.mainonly.chrom.size

    Alignment 序列比对

    工具简介：

    • bwa：用于将FASTQ格式的序列读段（reads）比对到参考基因组上。这是一个常用的高效比对工具，特别适用于高通量测序数据。
    • samtools：一个用于处理和操作SAM/BAM格式文件的工具包，支持格式转换、排序、索引等多种操作。

    用法：

    bwa mem -SP5M -t8 <index> <fastq1> <fastq2> | samtools view -bhS - > <outbam>

    指令：

    bwa mem -SP5M -t8 hg38.fasta input\_R1.fastq input\_R2.fastq | samtools view -bhS - > output.bam

    解释

    • bwa mem：BWA的一个子命令，用于执行序列比对，通常用于处理短读段数据。
    • -SP5M：这是BWA的Hi-C特定选项，用于处理Hi-C数据。它可能包括特殊的参数设置，以适应Hi-C数据的特点（如相对较长的片段和复杂的配对关系）。
    • -t8：指定BWA使用8个CPU核心来进行并行处理，以加速比对过程。
    • |：管道符号，将BWA的输出直接传递给下一个命令，即samtools。
    • samtools view -bhS -：
      • view 是samtools的一个子命令，用于将SAM格式转换为BAM格式。
      • -bhS 选项：
        • -b 表示输出BAM格式文件（比SAM格式更为紧凑，适合大规模数据处理）。
        • -h 保留SAM头部信息。
        • -S 表示输入是SAM格式（BWA mem的默认输出）。
        • 表示从标准输入（即来自BWA mem的输出）读取数据。

    Bam文件是存储读数的比对结果的二进制文件，可以通过以下指令查看文件，按Q退出查看模式

    samtools view output.bam | less -S

    Filtering 过滤

    Parse & Sort

    Parse: 分析，生成.pairsam文件
    Sort：排序
    可能存在的配对类型：https://pairtools.readthedocs.io/en/latest/formats.html#pair-types

    用法：

    samtools view -h <input\_bam> | pairsamtools parse -c <chromsize> -o <parsed\_pairsam>
pairsamtools sort --nproc 8 -o <sorted\_pairsam> <parsed\_pairsam>

    指令：

    #这一步会生成 parsed.pairsam.gz（中间文件）：将 BAM 文件转换成 PairsAM 格式并进行压缩
samtools view -h output.bam | pairsamtools parse -c hg38.mainonly.chrom.size -o parsed.pairsam.gz

#这一步会生成 sorted.pairsam.gz （输出文件）：将 PairsAM 文件进行排序，并将其保存为压缩文件。排序后的 PairsAM 文件可以方便后续的 Hi-C 数据分析操作
pairsamtools sort --nproc 8 -o sorted.pairsam.gz parsed.pairsam.gz

#这一步会解压并打开输出的文件,按Q可以退出
gunzip -c parsed.pairsam.gz |less -S

    Dedup & Select & Split

    Dedup：去重
    Select：过滤
    Split：生成pairs或bam文件
    用法：

    pairsamtools dedup --mark-dups -o <deduped_pairsam> <sorted_pairsam>
pairsamtools select <condition> -o <filtered_pairsam> <deduped_pairsam> 
pairsamtools split --output-pairs <output_pairs> <filtered_pairsam>

    指令：

    pairsamtools dedup --mark-dups -o deduped.pairsam.gz sorted.pairsam.gz
pairsamtools select '(pair_type == "UU") or (pair_type == "UR") or (pair_type == "RU")' -o filtered.pairsam.gz deduped.pairsam.gz 
pairsamtools split --output-pairs output.pairs.gz filtered.pairsam.gz

#显示输出文件
gunzip -c  output.pairs.gz |less -S

    Pairs文件命令示例:

    • bgzip output.pairs: 如果你的文件未压缩，可以使用 bgzip 对其进行 gzip 压缩。
    • pairix -f output.pairs.gz: 使用 pairix 工具对 pairs 文件进行索引，以便快速查询。（输出px2文件）
    • pairix output.pairs.gz ‘chr2:1-60000000|chr4:5000000-6000000’: 查询 pairs 文件中位于 chr2:1-60000000 和 chr4:5000000-6000000 区域之间的接触事件。

    Binning 分区块

    注：后来发现所使用的docker镜像中的cooler版本太老了，为了方便后续使用cooltools分析数据，换用随cooltools安装的新版cooler，详见对应文章。安装好后可以继续以下步骤，不需要在docker容器内进行。

    cooler文档：https://cooler.readthedocs.io/en/latest/

    这一步决定最终接触矩阵的分辨率

    https://cooler.readthedocs.io/en/latest/cli.html#cooler-cload-pairix

    用法：

    cooler cload pairix <chrsize\_file>:<bin\_size> <pairs\_file> <out\_prefix>.cool

    指令：

    cooler cload pairix hg38.mainonly.chrom.size:500000 output.pairs.gz output.cool

    输出为.cool文件

    Normalization 正则化

    https://cooler.readthedocs.io/en/latest/cli.html#cooler-balance

    用法：

    cooler balance <cool\_file>

    指令：

    cooler balance output.cool

    Aggregation 汇总

    https://cooler.readthedocs.io/en/latest/cli.html#cooler-zoomify

    用法：

    cooler zoomify <cool\_file>

    指令：

    cooler zoomify output.cool

    输出mcool文件
    到这一步就可以在Docker中关闭4dn-hic的容器了

    HiGlass 可视化

    文档：https://docs.higlass.io/tutorial.html

    HiGlass-Docker：https://github.com/higlass/higlass-docker

    首次运行higlass-docker：

    docker run --detach \
           --publish 8888:80 \
           --volume ~/hg-data:/data \
           --volume ~/hg-tmp:/tmp \
           --name higlass-container \
           higlass/higlass-docker

    现在可以在http://localhost:8888/app打开HiGlass

    将文件注入HiGlass容器

    #将mcool文件拷贝到容器的tmp目录
docker cp /Users/lidonglin/Documents/HiC_Exercise/data/output.mcool higlass-container:/tmp/output.mcool

#检查文件是否在tmp目录
docker exec higlass-container ls /tmp

#将该文件注入
docker exec higlass-container python higlass-server/manage.py ingest_tileset --filename /tmp/output.mcool --filetype cooler --datatype matrix

    现在在HiGlass页面中点击“+” - Add Track - center - 选中output.mcool即可

    HiGlass可视化结果

    Juicer / Juicebox

    使用HiCLift将.cool转换为.hic: HiCLift GitHub Repository

    在(base)环境中执行了以下指令：

    HiCLift --input <COOL_FILE> --input-format cooler --out-pre <OUTPUT_FILENAME> --output-format hic --out-chromsizes <CHROMSIZE_FILE> --in-assembly <ASSEMBLY_NAME_hg38> --out-assembly <ASSEMBLY_NAME_hg38> --memory 2G

    过程中安装了一堆东西

    pip install hiclift
pip install hicstraw
pip install cooler
pip install pairtools
conda install mamba (好像没用)
pip install path/to/downloaded/kerneltree.tar.gz  (installed from PyPI website by downloading .tar.gz and modifying setup.py license=["GPL2"] to license="GPL2")
conda install cxx-compiler
Java (install from official website)
brew install htslib
conda install pairix

    相关工具

    Powered by Hexo
    Theme based on Minima by Adi Sakti Jrs
    Published on GitHub