全基因组 - 人类基因组变异分析（PacBio) （3）-- pbmm2

原创

三代测序说

修改于 2023-11-23 08:10:36

6781

修改于 2023-11-23 08:10:36

文章被收录于专栏：三代测序-说三代测序-说

一. 长读段比对 (Long-Read Mapping)常用的比对软件

长读段比对算法与一代/二代测序数据的比对算法有很大的不同，因为长读段通常更长、包含更多错误和变异，并且需要更复杂的比对策略。

PacBio三代测序数据有很多比对软件可供选择，例如：

LAST：https://github.com/mcfrith/last-genome-alignments

BlasR：https://github.com/mchaisso/blasr

BLASR是第一个针对PacBio序列的比对工具，2012年发表在《BMC Bioinformatics》期刊上，由PacBio研究团队开发，最新版本停在2019.4.11，目前Google引用次数为1170（截止2023.10.25）。

BWA-MEM：https://github.com/lh3/bwa

BWA老牌比对软件了。BWA-MEM 是一种新的比对算法，用于将测序 reads 或者组装后 contigs 比对至大型参考基因组，例如人参考基因组。它该算法对测序错误有良好的稳定性，适用的reads 长度范围较广，从70bp至几Mb。

GraphMap：https://github.com/isovic/graphmap

MECAT：https://github.com/xiaochuanle/MECAT

中山大学研究团队开发的新工具MECAT，集超快比对、校正、组装于一体，可提高三代测序数据序列比对，校正和组装的运算速度，降低计算资源的消耗。

minimap2：https://github.com/lh3/minimap2

Minimaps2是李恒大神在2018年发表在Bioinformatics上的一款针对三代数据开发的比对工具。

NGMLR: https://github.com/philres/ngmlr

NextGenMap-LR(ngmlr)主要用于三代测序的长reads（PacBio 、Oxford Nanopore）与参考基因组的比对。

pbmm2: https://github.com/PacificBiosciences/pbmm2 等等。

现在较为常用的是pbmm2、minimap2、NGMLR，其中pbmm2是PacBio官方基于minimap2进行优化的版本。

二. pbmm2的使用教程

在得到sample.CCS.bam文件后，因为HiFi数据质量较高，一把不需要额外的质控步骤，就可以将HiFi数据和下载的参考基因组序列进行比对了。

1.参考基因组的获取

分析前，除测序数据外，我们还需准备对应物种的参考基因组fasta文件。对此可以根据自己研究的需要，在NCBI、Ensembl、UCSC等常见数据库中进行下载。

主流的基因注释版本有三种：RefSeq/Ensembl/UCSC

Refseq=NCBI；Ensembl=Gencode

Ensemble注释更全面，Refseq适合那些不那么复杂的注释。

对于人类和小鼠来说，我们还可以从Gencode数据库中进行下载。Gencode综合HAVANA和Ensembl数据库中的信息，通过实验手段加以验证，从而构建了一个高质量的注释信息数据库。Gencode的注释来源于两部分。分别是Ensembl-Havana团队生成的手动基因注释和Ensembl-genebuild的自动基因注释。

gencode中标识HAVANA来源的，这表示它是人工注释的。但是这些注释也有可能是由于Havana手动注释和Ensembl自动注释合并的结果;

而如果标识的是ENSEMBL，则表明这条注释是由的确是Ensembl自动注释得到的。

以下载人类参考基因组为例：

Gencode：https://www.gencodegenes.org/human/

点击图1 所示的 fasta 连接即可下载完整的人类参考基因组（GRCh38.p14）。

Ensembl: https://useast.ensembl.org/Homo_sapiens/Info/Index

点击图2 所示的 Download DNA sequence (FASTA) 连接进入下载页面，进一步选择Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz , 即可下载完整的人类参考基因组（GRCh38.p14）。

注释：

Ensembl提供的参考基因组有2种组装形式和3种重复序列处理方式，分别是primary,toplevel,unmasked(dna),soft-masked(dna_sm),masked(dna_rm)。一般选择dna.primary和dna_sm.primary。 Toplevel: 这样的数据包含了所有的序列区域（比如染色体、非染色体以及用大量N填充的单倍型haplotypes或基因组补丁patches区域），比如：Homo_sapiens.GRCh37.dna.toplevel.fa.gz Primary Assembly：在上面toplevel的基础上，排除了单倍型或基因组补丁区域。如果看到目录中不存在这种类型的数据（比如这里果蝇就没有，而人类的基因组数据就存在），那么就意味着基因组不包含单倍型或基因组补丁区域，其实也就是等同于Toplevel.

dna：原原本本的DNA序列。

dna_rm（hard_mask）：利用RepeatMasker工具将重复区域和低复杂度区域标记成一串N，会造成信息的丢失。

dna_sm（soft_mask）：所有重复区域和低复杂度区域替换为小写的碱基

RefSeq：https://www.ncbi.nlm.nih.gov/datasets/taxonomy/9606/

点击图3 中的Download即可下载人类基因参考组（GRCh38.p14）。

UCSC genome browser: https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/

点击图4 中hg38.fa.gz即可下载人类基因参考组（GRCh38.p14）。

2. pbmm2安装

#使用conda安装pbmm2
$ conda install -c bioconda pbmm2

#安装版本
v1.13.0

3. pbmm2使用

建立人类参考基因组索引

Index: index reference and store as .mmi file Usage: pbmm2 index options <ref.fa|xml> <out.mmi> --preset 比对模式默认为CCS, 如果为CCS - HiFi数据，则不用设置

#建立基因组文件夹
$ mkdir human_ref
#将参考基因组GRCh38.fa拷入human_ref 文件夹, 进入human_ref
$ cd  human_ref

#进入参考基因组所在文件夹，对参考序列进行索引
$ pbmm2 index GRCh38.fa GRCh38.mmi   #此过程在我们计算服务上比较快

#如果运行时间较长，使用 nohup 加 & 将程序不挂断运行并放入后台
$ nohup pbmm2 index GRCh38.fa GRCh38.mmi &

数据与人类参考基因组进行比对

Usage: pbmm2 align options <ref.fa|xml|mmi> <in.bam|xml|fa|fq> out.aligned.bam|xml 除了PacBio默认的 bam文件， fastq / fasta 作为输入文件也是可以的。

#比对参考序列，-j -J 如果在服务器空闲时可以不指定
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.log --sample sample_name \
human_ref/GRCh38.mmi sample.ccs.bam sample.align.bam 

# --log-level INFO，给出基本mapping统计结果和时间
# --sample 指定样品名称；-j 指定比对线程；-J 8 指定排序步骤线程，8为最大；--sort 输出排序后的bam，--preset选择参数集-默认

pbmm2 index 帮助文档:

pbmm2 index - Index reference and store as .mmi file

Usage:
  pbmm2 index [options] <ref.fa|xml> <out.mmi>

  ref.fa|xml                STR   Reference FASTA, ReferenceSet XML
  out.mmi                   STR   Output Reference Index

Parameter Set Option:
  --preset                  STR   Set alignment mode. See below for preset parameter details. Valid choices: (SUBREAD,
                                  CCS, ISOSEQ, UNROLLED). [CCS]

Parameter Override Options:
  -k                        INT   k-mer size (no larger than 28). [-1]
  -w                        INT   Minimizer window size. [-1]
  -u,--no-kmer-compression        Disable homopolymer-compressed k-mer (compression is active for SUBREAD & UNROLLED
                                  presets).

  -h,--help                       Show this help and exit.
  --version                       Show application version and exit.
  -j,--num-threads          INT   Number of threads to use, 0 means autodetection. [0]
  --log-level               STR   Set log level. Valid choices: (TRACE, DEBUG, INFO, WARN, FATAL). [WARN]
  --log-file                FILE  Log to a file, instead of stderr.

Alignment modes of --preset:
    SUBREAD     : -k 19 -w 10
    CCS or HiFi : -k 19 -w 10 -u
    ISOSEQ      : -k 15 -w 5  -u
    UNROLLED    : -k 15 -w 15

Copyright (C) 2004-2023     Pacific Biosciences of California, Inc.
This program comes with ABSOLUTELY NO WARRANTY; it is intended for
Research Use Only and not for use in diagnostic procedures.

pbmm2 align帮助文档:

pbmm2 align - Align PacBio reads to reference sequences

Usage:
  pbmm2 align [options] <ref.fa|xml|mmi> <in.bam|xml|fa|fq|gz|fofn> [out.aligned.bam|xml]

  ref.fa|xml|mmi             STR    Reference FASTA, ReferenceSet XML, or Reference Index
  in.bam|xml|fa|fq|gz|fofn   STR    Input BAM, DataSet XML, FASTA, or FASTQ
  out.aligned.bam|xml        STR    Output BAM or DataSet XML

Basic Options:
  --chunk-size               INT    Process N records per chunk. [100]

Sorting Options:
  --sort                            Generate sorted BAM file.
  -m,--sort-memory           STR    Memory per thread for sorting. [768M]
  -J,--sort-threads          INT    Number of threads used for sorting; 0 means 25% of -j, maximum 8. [0]

Parameter Set Options:
  --preset                   STR    Set alignment mode. See below for preset parameter details. Valid choices:
                                    (SUBREAD, CCS, HIFI, ISOSEQ, UNROLLED). [CCS]

General Parameter Override Options:
  -k                         INT    k-mer size (no larger than 28). [-1]
  -w                         INT    Minimizer window size. [-1]
  -u,--no-kmer-compression          Disable homopolymer-compressed k-mer (compression is active for SUBREAD & UNROLLED
                                    presets).
  -A                         INT    Matching score. [-1]
  -B                         INT    Mismatch penalty. [-1]
  -z                         INT    Z-drop score. [-1]
  -Z                         INT    Z-drop inversion score. [-1]
  -r                         INT    Bandwidth used in chaining and DP-based alignment. [-1]
  -g                         INT    Stop chain enlongation if there are no minimizers in N bp. [-1]

Gap Parameter Override Options (a k-long gap costs min{o+k*e,O+k*E}):
  -o,--gap-open-1            INT    Gap open penalty 1. [-1]
  -O,--gap-open-2            INT    Gap open penalty 2. [-1]
  -e,--gap-extend-1          INT    Gap extension penalty 1. [-1]
  -E,--gap-extend-2          INT    Gap extension penalty 2. [-1]

IsoSeq Parameter Override Options:
  -G                         INT    Max intron length (changes -r). [-1]
  -C                         INT    Cost for a non-canonical GT-AG splicing (effective in ISOSEQ preset). [-1]
  --no-splice-flank                 Do not prefer splice flanks GT-AG (effective in ISOSEQ preset).

Read Group Options:
  --sample                   STR    Sample name for all read groups. Defaults, in order of precedence: SM field in
                                    input read group, biosample name, well sample name, "UnnamedSample".
  --rg                       STR    Read group header line such as '@RG\tID:xyz\tSM:abc'. Only for FASTA/Q inputs.

Identity Filter Options (combined with AND):
  -y,--min-gap-comp-id-perc  FLOAT  Minimum gap-compressed sequence identity in percent. [70]

Output Options:
  -l,--min-length            INT    Minimum mapped read length in basepairs. [50]
  -N,--best-n                INT    Output at maximum N alignments for each read, 0 means no maximum. [0]
  --strip                           Remove all kinetic and extra QV tags. Output cannot be polished.
  --split-by-sample                 One output BAM per sample.
  --unmapped                        Include unmapped records in output.
  --bam-index                STR    Generate index for sorted BAM output. Valid choices: (NONE, BAI, CSI). [BAI]
  --short-sa-cigar                  Populate SA tag with short cigar representation.

Input Manipulation Options (mutually exclusive):
  --median-filter                   Pick one read per ZMW of median length.
  --zmw                             Process ZMW Reads, subreadset.xml input required (activates UNROLLED preset).
  --hqregion                        Process HQ region of each ZMW, subreadset.xml input required (activates UNROLLED
                                    preset).

Sequence Manipulation Options:
  --collapse-homopolymers           Collapse homopolymers in reads and reference.

  -h,--help                         Show this help and exit.
  --version                         Show application version and exit.
  -j,--num-threads           INT    Number of threads to use, 0 means autodetection. [0]
  --log-level                STR    Set log level. Valid choices: (TRACE, DEBUG, INFO, WARN, FATAL). [WARN]
  --log-file                 FILE   Log to a file, instead of stderr.

Alignment modes of --preset:
    SUBREAD     : -k 19 -w 19    -o 5 -O 56 -e 4 -E 1 -A 2 -B 5 -z 400 -Z 50  -r 2000   -g 5000
    CCS or HiFi : -k 19 -w 19 -u -o 6 -O 26 -e 2 -E 1 -A 1 -B 4 -z 400 -Z 50  -r 2000   -g 5000
    ISOSEQ      : -k 15 -w 5  -u -o 2 -O 32 -e 1 -E 0 -A 1 -B 2 -z 200 -Z 100 -r 200000 -g 2000 -C 5 -G 200000
    UNROLLED    : -k 15 -w 15    -o 2 -O 32 -e 1 -E 0 -A 1 -B 2 -z 200 -Z 100 -r 2000   -g 10000

Copyright (C) 2004-2023     Pacific Biosciences of California, Inc.
This program comes with ABSOLUTELY NO WARRANTY; it is intended for
Research Use Only and not for use in diagnostic procedures.

5. 公共数据演示：

(1) 从gencode数据库下载人类参考基因组, 进行pbmm2索引。

PacBio推荐人类参考基因组（详细参照李恒博客），所以采用推荐基因组进行后续分析。

李恒2017年年底博客 https://lh3.github.io/2017/11/13/which-human-reference-genome-to-use

推荐的人类参考基因组：GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

#新建文件夹
$ mkdir Human_ref

#下载参考基因组 GRCh38
$ wget -c -t 0 ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

#对参考基因组进行解压
$ gunzip  GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

#对参考基因组从新命名
$ mv GCA_000001405.15_GRCh38_no_alt_analysis_set.fna GRCh38.fa

#pbmm2 index
$ nohup pbmm2 index GRCh38.fa GRCh38.mmi

(2) 下载示例数据。

Example Datasets

如图6所示：下载示例人类基因组数据。

德系犹太人家系：HG002（子）、HG003（父）、HG004（母），属于个人基因组计划中的样本。

HG002_1：m84011_220902_175841_s1.hifi_reads.bam

HG003：m84010_220919_235306_s2.hifi_reads.bam

HG004：m84010_220919_232145_s1.hifi_reads.bam

(3) 示例数据与人类参考基因组进行比对。

#HG002_1
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.HG002_1.log --sample HG002_1 \
Human_ref/GRCh38.mmi m84011_220902_175841_s1.hifi_reads.bam HG002_1.align.bam 

#HG003
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.HG003.log --sample HG003 \
Human_ref/GRCh38.mmi m84010_220919_235306_s2.hifi_reads.bam HG003.align.bam 

#HG004
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.HG004.log --sample HG004 \
Human_ref/GRCh38.mmi m84010_220919_232145_s1.hifi_reads.bam HG004.align.bam