打嗝

一个轻量级库，可生成和处理与冷却器兼容的 2Dbedgraph 或instaGRAAL格式的 Hi-C 接触图。它本质上是yahcp管道、hicstuff库和3C 教程和DADE 管道中说明的额外功能的合并，所有这些都打包在一起以提供额外的便利。

目标是使 Hi-C 矩阵的生成和操作尽可能简单，并适用于任何生物体。

目录

• 安装
• 用法
  • 完整的管道
  • 单个组件
• 图书馆
• 连接模块
• 文件格式
• 贡献

安装

要安装稳定版本：

pip3 install -U hicstuff

或者，对于最新的开发版本：

    pip3 install -e git+https://github.com/koszullab/hicstuff.git@master#egg=hicstuff

外部依赖

该实用程序需要 Bowtie2 和/或 minimap2 以及 samtools pipeline。

您可以通过 conda 包管理器安装它们：

conda install -c bioconda minimap2 bowtie2 samtools

或者，在 ubuntu 上，您还可以通过 APT 安装它们以及其他依赖项：

apt-get install samtools bowtie2 minimap2 libbz2-dev liblzma-dev

码头工人安装

dockerhub 上提供了一个预构建的 docker 映像，可以使用以下命令运行：

docker run koszullab/hicstuff

用法

hicstuff 命令行界面由多个子命令组成。您始终可以通过运行以下命令获取所有可用命令的摘要：

hicstuff --help

Simple Hi-C pipeline for generating and manipulating contact matrices.

usage:
    hicstuff [-hv] <command> [<args>...]

options:
    -h, --help                  shows the help
    -v, --version               shows the version

The subcommands are:
    digest          Digest genome into a list of fragments.
    distancelaw     Analyse and plot distance law.
    filter          Filters Hi-C pairs to exclude spurious events.
    iteralign       Iteratively aligns reads to a reference genome.
    pipeline        Hi-C pipeline to generate contact matrix from fastq files.
    rebin           Bin the matrix and regenerate files accordingly.
    subsample       Bootstrap subsampling of contacts from a Hi-C map.
    view            Visualize a Hi-C matrix.  

完整的管道

可以使用该hicstuff pipeline命令一次运行管道的所有组件。这允许从单个命令中的读取生成接触矩阵。默认情况下，输出是与 GRAAL 兼容的 COO 稀疏矩阵格式，但它可以是 2D 床位图或酷文件，而不是使用该--matfmt选项。更详细的文档可以在 readthedocs 网站上找到：https ://hicstuff.readthedocs.io/en/latest/index.html

usage:
    pipeline [--quality-min=INT] [--size=INT] [--no-cleanup] [--start-stage=STAGE]
             [--threads=INT] [--aligner=bowtie2] [--matfmt=FMT] [--prefix=PREFIX]
             [--tmpdir=DIR] [--iterative] [--outdir=DIR] [--filter] [--enzyme=ENZ]
             [--plot] [--circular] [--distance-law] [--duplicates] [--read-len=INT]
             [--centromeres=FILE] [--remove-centromeres=INT] --genome=FILE <input1> [<input2>]

arguments:
    input1:             Forward fastq file, if start_stage is "fastq", sam
                        file for aligned forward reads if start_stage is
                        "bam", or a .pairs file if start_stage is "pairs".
    input2:             Reverse fastq file, if start_stage is "fastq", sam
                        file for aligned reverse reads if start_stage is
                        "bam", or nothing if start_stage is "pairs".

例如，要使用 8 个线程使用 minimap2 运行管道并在目录中生成 instagraal 格式的矩阵out：

hicstuff pipeline -t 8 -a minimap2 -e DpnII -o out/ -g genome.fa reads_for.fq reads_rev.fq

如果您已经对齐读取，hicstuff 管道也可以将 bam 文件作为输入。例如，要生成一个具有 5kb 固定 bin 大小的酷格式矩阵：

# Note the bam files have to be name-sorted, this can be done using samtools
samtools sort -n aligned_for.bam -o namesorted_for.bam
samtools sort -n aligned_rev.bam -o namesorted_rev.bam
hicstuff pipeline -S bam -e 5000 -M cool -o out/ -g genome.fa namesorted_for.bam namesorted_rev.bam

管道也可以使用hicstuff.pipeline子模块从 python 运行。例如，这将使用 bowtie2（默认）使用迭代对齐运行管道并保留所有中间文件。有关使用 API 的更多示例，请参阅API 演示

from hicstuff import pipeline as hpi

hpi.full_pipeline(
    'genome.fa', 
    'end1.fq', 
    'end2.fq', 
    no_cleanup=True
    iterative=True
    out_dir='out', 
    enzyme="DpnII")

流水线的一般步骤如下：

单个组件

对于更高级的用法，可以在命令行上独立使用不同的脚本来执行管道的各个部分。本自述文件包含快速说明和示例用法。要获得任何子命令的详细说明，可以使用hicstuff <subcommand> --help.

迭代对齐

将 fastq 文件中的读取截断为 20 个碱基对，并迭代扩展和重新对齐未映射的读取，以优化 3C 库中唯一对齐读取的比例。

usage:
    iteralign [--aligner=bowtie2] [--threads=1] [--min_len=20]
              [--tempdir DIR] --out_sam=FILE --genome=FILE <reads.fq>

基因组的消化

根据限制酶或固定块大小将 fasta 文件消化成片段。在名为“info_contigs.txt”和“fragments_list.txt”的目标目录中生成两个输出文件

usage:
    digest [--plot] [--figdir=FILE] [--circular] [--size=INT]
           [--outdir=DIR] --enzyme=ENZ <fasta>

例如，用 MaeII 和 HinfI 消化酵母基因组并显示片段长度的直方图：

hicstuff digest --plot --outdir output_dir --enzyme MaeII,HinfI Sc_ref.fa

过滤 3C 事件

根据默认情况下从库中自动估计的最小距离阈值，从库中过滤虚假 3C 事件，例如循环和未剪辑。还可以绘制 3C 库统计信息。该模块将具有 9 列的对文件作为输入（readID、chr1、pos1、chr2、pos2、strand1、strand2、frag1、frag2）并对其进行过滤。

usage:
    filter [--interactive | --thresholds INT-INT] [--plot]
           [--figdir FILE] <input.pairs> <output.pairs>

查看联系地图

将 Hi-C 矩阵文件可视化为接触频率的热图。允许通过对矩阵进行分箱和归一化来调整可视化，并将输出图像保存到磁盘。如果未指定输出，则以交互方式显示输出。如果提供了两个接触图，将显示第一个除以第二个的对数比。

usage:
    view [--binning=1] [--despeckle] [--frags FILE] [--trim INT] [--n-mad FLOAT]
         [--normalize] [--max=99] [--output=IMG] [--cmap=CMAP] [--dpi=INT]
         [--transform=FUN] [--circular] [--region=STR] <contact_map> [<contact_map2>]

arguments:
    contact_map             Sparse contact matrix in bg2, cool or graal format
    contact_map2            Sparse contact matrix in bg2, cool or graal format,
                            if given, the log ratio of contact_map/contact_map2
                            will be shown.

例如，要从以 10kb 重组的全基因组 Hi-C 矩阵查看 1 号染色体的 1Mb 区域：

    hicstuff view --normalize --binning 10kb --region chr1:10,000,000-11,000,000 --frags fragments_list.txt contact_map.tsv

图书馆

hicstuff 程序的所有组件都可以用作 python 模块。请参阅rethedocs上的文档。该库的预期联系人地图格式是一个简单的 CSV 文件，该库处理的对象是简单的numpy数组。hicstuff 的各种子模块包含各种实用程序。

import hicstuff.digest # Functions to work with restriction fragments
import hicstuff.iteralign # Functions related to iterative alignment
import hicstuff.hicstuff # Contains utilities to modify and operate on contact maps as numpy arrays
import hicstuff.filter # Functions for filtering 3C events by type (uncut, loop)
import hicstuff.view # Utilities to visualise contact maps
import hicstuff.io # Reading and writing hicstuff files
import hicstuff.pipeline # Generation and processing of files to generate matrices.

连接模块

此处描述的所有步骤在运行hicstuff pipeline. 但是如果你想手动连接不同的模块，可以使用一些 python 脚本来处理中间输入和输出文件。

对齐读数

您可以使用自己喜欢的读取映射软件、使用命令行实用程序hicstuff iteralign或使用align_reads子模块中的辅助函数独立生成 SAM 文件hicstuff.pipeline。例如，要使用 minimap2（而不是默认的 bowtie2）执行迭代对齐：

使用python函数：

from hicstuff import pipeline as hpi

hpi.align_reads("end1.fastq", "genome.fasta", "end1.bam", iterative=True, minimap2=True)

使用命令行工具：

hicstuff iteralign --minimap2 --iterative -f genome.fasta -o end1.sam end1.fastq

从对齐中提取联系人

的输出hicstuff iteralign是一个 SAM 文件。为了检索 Hi-C 对，您需要在两个 fastq 文件上分别运行 iteralign，并使用pipelinehicstuff 的子模块将生成的对齐文件处理为按名称排序的 BAM 文件，如下所示。

from hicstuff import pipeline as hpi
import pysam as ps
# Sort alignments by read names and get into BAM format
ps.sort("-n", "-O", "BAM", "-o", "end1.bam.sorted", "end1.sam")
ps.sort("-n", "-O", "BAM", "-o", "end2.bam.sorted", "end2.sam")
# Combine BAM files
hpi.bam2pairs("end1.sorted.bam", "end2.sorted.bam", "output.pairs", "info_contigs.txt", min_qual=30)

这将生成一个“pairs”文件，其中包含所有读取对，其中两个读取都与至少 30 的映射质量对齐。

将每次读取归因于限制片段

为了建立联系矩阵，我们需要将每个读取归因于基因组中的一个片段。这是通过针对基因组中的限制性位点列表对每个读取位置执行二进制搜索来完成的。

from hicstuff import digest as hcd
from Bio import SeqIO

# Build a list of restriction sites for each chromosome
restrict_table = {}
for record in SeqIO.parse("genome.fasta", "fasta"):
    # Get chromosome restriction table
    restrict_table[record.id] = hcd.get_restriction_table(
        record.seq, enzyme, circular=circular
    )

# Add fragment index to pairs (readID, chr1, pos1, chr2,
# pos2, strand1, strand2, frag1, frag2)
hcd.attribute_fragments("output.pairs", "output_indexed.pairs", restrict_table)

过滤对

然后可以在命令行中使用hicstuff filter命令或在 python 中使用hicstuff.filter子模块过滤生成的对文件。否则，可以直接从未过滤的对构建矩阵。

在命令行上过滤：

hicstuff filter output_indexed.pairs output_filtered.pairs

在python中过滤：

from hicstuff import filter as hcf

uncut_thr, loop_thr = hcf.get_thresholds("output_indexed.pairs")
hcf.filter_events("output_indexed.pairs", "output_filtered.pairs", uncut_thr, loop_thr)

请注意，命令和 python 函数都有各种选项来生成图形或调整过滤阈值。这些选项可以使用显示hicstuff filter -h

矩阵生成

然后可以使用 python 子模块生成 Hi-C 稀疏接触矩阵hicstuff.pipeline。该矩阵可以以 GRAAL 兼容的 COO 格式、bedgraph2 或酷格式生成。

from hicstuff import pipeline as hpi

n_frags = sum(1 for line in open(fragments_list, "r")) - 1
hpi.pairs2matrix("output_filtered.pairs", "abs_fragments_contacts_weighted.txt", 'fragments_list.txt', mat_fmt="GRAAL")

文件格式

• 对文件：此格式用于管道中的所有中间文件，也被hicstuff filter. 它是一种制表符分隔的格式，包含有关 Hi-C 对的信息。它具有由 4D Nucleome 数据协调和集成中心定义的官方规范。
• 2D 床位图：这是hicstuff pipeline稀疏矩阵的可选输出格式。它每行有两个片段，以及它们一起被发现的次数。它具有以下字段：chr1、start1、end1、chr2、start2、end2、发生
  • 这些文件可以由cooler加载，其中chrom.sizescooler load -f bg2 <chrom.sizes>:<binsize> in.bg2.gz out.cool是一个制表符分隔的文件，每行都有染色体名称和长度，binsize是矩阵中bin的大小。
• GRAAL 稀疏矩阵：这是一个简单的制表符分隔文件，包含 3 列：frag1、frag2、contacts。id 列对应于限制片段的绝对 id（0 索引）。第一行是包含矩阵中行数、列数和非零条目数的标题。例子：

564	564	6978
0	0	3
1	2	4
1	3	3

• Fragments_list.txt：此制表符分隔文件提供有关限制片段位置、大小和 GC 内容的信息。请注意，坐标是 0 点碱基对，与具有 1 点碱基对的对格式不同。例子：
  • id：染色体内基于 1 的限制性片段索引。
  • chrom：染色体标识符。顺序应与 info_contigs.txt 或对文件中的相同。
  • start_pos：从 0 开始的片段开始，以碱基对表示。
  • end_pos：从 0 开始的片段结尾，以碱基对表示。
  • size：片段的大小，以碱基对为单位。
  • gc_content：片段中 G 和 C 核苷酸的比例。

id	chrom	start_pos	end_pos	size	gc_content
1	seq1	0	21	21	0.5238095238095238
2	seq1	21	80	59	0.576271186440678
3	seq1	80	328	248	0.5201612903225806

• info_contigs.txt：此制表符分隔文件提供有关 contig 的信息，例如限制片段的数量和大小。例子：
  • contig：染色体鉴定。成对文件或 Fragments_list.txt 中的顺序应相同。
  • 长度：染色体长度，以碱基对计。
  • n_frags：染色体中限制性片段的数量。
  • cumul_length：先前染色体的累积长度，以碱基对为单位。

contig	length	n_frags	cumul_length
seq1	60000	409	0
seq2	20000	155	409

贡献

欢迎所有贡献，以错误报告、建议、文档或拉取请求的形式。在记录函数时，我们对文档字符串使用numpy 标准。

我们使用的代码格式化标准是black， --line-length=79 遵循 PEP8 建议。我们使用pytest和pytest-doctest插件pytest-pylint作为我们的测试框架。理想情况下，新函数应该有相关的单元测试，放在tests文件夹中。要测试代码，您可以运行：

pytest --doctest-modules --pylint --pylint-error-types=EF --pylint-rcfile=.pylintrc hicstuff tests

引文

请使用官方 DOI 引用 hicstuff，如下所示：

Cyril Matthey-Doret、Lyam Baudry、Amaury Bignaud、Axel Cournac、Remi-Montagne、Nadège Guiglielmoni、Théo Foutel Rodier 和 Vittore F. Scolari。2020. hicstuff：用于生成和处理 Hi-C 数据的简单库/管道。泽诺多。http://doi.org/10.5281/zenodo.4066363

中文条目：

@software{cyril_matthey_doret_2020_4066351,
  author       = {Cyril Matthey-Doret and
                  Lyam Baudry and
                  Amaury Bignaud and
                  Axel Cournac and
                  Remi-Montagne and
                  Nadège Guiglielmoni and
                  Théo Foutel-Rodier and
                  Vittore F. Scolari},
  title        = {hicstuff: Simple library/pipeline to generate and handle Hi-C data },
  month        = oct,
  year         = 2020,
  publisher    = {Zenodo},
  version      = {v2.3.1},
  doi          = {10.5281/zenodo.4066351},
  url          = {http://doi.org/10.5281/zenodo.4066363}
}