生信软件25 - 三代测序数据灵敏比对工具ngmlr

1. ngmlr简介

CoNvex Gap-cost Ments for Long Reads（ngmlr）是一种长reads比对工具，可以将PacBio或Oxford Nanopore灵敏地与（大）参考基因组（比如人类参考基因组）对齐，能快速和正确地比对reads，包括那些跨越（复杂）结构变异的reads。

Ngmlr使用结构变异（SV）感知的k-mer搜索来找到reads的近似映射位置，然后使用Smith-Waterman比对算法来计算最终比对。Ngmlr使用凸间隙成本模型（convex gap cost model），该模型对较长间隙的间隙扩展的罚分小于较短间隙的间隙扩展，以计算精确的比对。

差距模型允许ngmlr同时考虑测序误差和真实的基因组变异，并使其在更精确地鉴定源于结构变异的断裂点的位置方面特别有效。k-mer搜索有助于检测和拆分不能线性比对的reads，使ngmlr能够可靠地将读段与广泛的不同结构变异（包括嵌套SV）（例如侧翼缺失的倒位）进行比对。

NGMLR和Sniffles可以自动过滤错误事件并对低覆盖率数据，从而降低了阻碍长reads在临床和研究环境中应用的高成本。

参考文献：
Accurate detection of complex structural variations using single-molecule sequencing
https://www.nature.com/articles/s41592-018-0001-7

2. ngmlr安装

# 下载安装包
wget https://github.com/philres/ngmlr/releases/download/v0.2.7/ngmlr-0.2.7-linux-x86_64.tar.gz

# 解压
tar xvzf ngmlr-0.2.7-linux-x86_64.tar.gz

cd ngmlr-0.2.7/

3. PacBio数据运行

三代fastq数据质控参考：三代测序PacBio&ONT reads过滤和修剪-Chooper和NanoFilt

结构变异Calling参考：三代测序结构变异分析 - 单样本Germline SV calling和多样本SV Calling

# 数据质控
gunzip -c sample.fastq.gz | NanoFilt -q 10 -l 1000 > sample.clean.fastq

# 比对
ngmlr -t 4 -r reference.fasta \
-q sample.clean.fastq -o sample.sam
# -t: 线程数
# -r: 参考基因组fasta路径
# -q: clean fastq路径
# -o: 比对输出sam文件路径

# Sam2bam
samtools view -bS sample.sam > sample.bam

# 排序
samtools sort sample.bam -o sample.sorted.bam
samtools index sample.sorted.bam

# 查看bam
samtools view sample.sorted.bam|less -S

4. Oxford Nanopore运行

gunzip -c sample.fastq.gz | NanoFilt -q 10 -l 1000 > sample.clean.fastq

# 参数-x: 三代数据类型，默认PacBio
ngmlr -t 4 -r reference.fasta \
-q sample.clean.fastq -o test.sam -x ont


# Sam2bam
samtools view -bS sample.sam > sample.bam

# 排序
samtools sort sample.bam -o sample.sorted.bam
samtools index sample.sorted.bam

# 查看bam
samtools view sample.sorted.bam|less -S

主要参数说明

Usage: ngmlr [options] -r <reference> -q <reads> [-o <output>]

Input/Output:
    -r <file>,  --reference <file>
        (required)  Path to the reference genome (FASTA/Q, can be gzipped)
    -q <file>,  --query <file>
        Path to the read file (FASTA/Q) [/dev/stdin]
    -o <string>,  --output <string>
        Path to output file [stdout]
    --skip-write
        Don't write reference index to disk [false]
    --bam-fix
        Report reads with > 64k CIGAR operations as unmapped. Required to be compatible with the BAM format [false]
    --rg-id <string>
        Adds RG:Z:<string> to all alignments in SAM/BAM [none]
    --rg-sm <string>
        RG header: Sample [none]
    --rg-lb <string>
        RG header: Library [none]
    --rg-pl <string>
        RG header: Platform [none]
    --rg-ds <string>
        RG header: Description [none]
    --rg-dt <string>
        RG header: Date (format: YYYY-MM-DD) [none]
    --rg-pu <string>
        RG header: Platform unit [none]
    --rg-pi <string>
        RG header: Median insert size [none]
    --rg-pg <string>
        RG header: Programs [none]
    --rg-cn <string>
        RG header: sequencing center [none]
    --rg-fo <string>
        RG header: Flow order [none]
    --rg-ks <string>
        RG header: Key sequence [none]

General:
    -t <int>,  --threads <int>
        Number of threads [1]
    -x <pacbio, ont>,  --presets <pacbio, ont>
        Parameter presets for different sequencing technologies [pacbio]
    -i <0-1>,  --min-identity <0-1>
        Alignments with an identity lower than this threshold will be discarded [0.65]
    -R <int/float>,  --min-residues <int/float>
        Alignments containing less than <int> or (<float> * read length) residues will be discarded [0.25]
    --no-smallinv
        Don't detect small inversions [false]
    --no-lowqualitysplit
        Split alignments with poor quality [false]
    --verbose
        Debug output [false]
    --no-progress
        Don't print progress info while mapping [false]