一、Minimap2 使用

官方網(wǎng)站：https://github.com/lh3/minimap2

1. 軟件安裝

首先從github官網(wǎng)上下載minimap2的二進(jìn)制文件壓縮包，minimap2-2.26_x64-linux.tar.bz2，然后上傳到服務(wù)器上。

# minimap2，v2.26壓縮包解壓縮
$ tar -xjvf minimap2-2.26_x64-linux.tar.bz2

# -x 解壓
# -j 有bz2屬性的
# -v 顯示所有過(guò)程
# -f 使用檔案名字，切記，這個(gè)參數(shù)是最后一個(gè)參數(shù)，后面只能接檔案名。

#將minimap2加入絕對(duì)路徑
$ echo 'PATH=/path/to/minimap2-2.26_x64-linux:$PATH' >> ~/.bashrc && source ~/.bashrc

$ echo 'PATH=/mnt/data/home/mli/Desktop/Software/minimap2-2.26_x64-linux:$PATH' >> ~/.bashrc && source ~/.bashrc

2. 軟件使用

• 數(shù)據(jù)格式轉(zhuǎn)化

公共數(shù)據(jù)只有 .bam 格式，所以要先將.bam 文件轉(zhuǎn)化為.fastq文件才能輸入minimap2；如果直接獲得.fastq文件則可以省略此步轉(zhuǎn)化。

數(shù)據(jù)格式轉(zhuǎn)化所用的程序?yàn)?code>BAM2fastx（PacBio官方工具），PacBio將一系列工具，包括對(duì).bam文件進(jìn)行索引的pbindex，都放在pbtk（pb tool kit）中，所以運(yùn)行以下命令全部安裝：

BAM2fastx tools ：https://github.com/PacificBiosciences/bam2fastx

# conda安裝pbtk
$ conda install -c bioconda pbtk

Example Datasets
德系猶太人家系：HG002（子）、HG003（父）、HG004（母），屬于個(gè)人基因組計(jì)劃中的樣本。

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

對(duì).bam文件進(jìn)行索引（生成.bam.pbi文件），然后進(jìn)行.bam 文件到.fastq文件的轉(zhuǎn)換：

# generates out.fastq.gz
$ bam2fastq -o out in_1.bam in_2.bam in_3.xml in_4.bam

#因?yàn)檫\(yùn)行bam2fastq 需要對(duì)bam文件先進(jìn)行索引
$ pbindex m84010_220919_232145_s1.hifi_reads.bam
$ bam2fastq -o m84010_220919_232145_s1.hifi_reads m84010_220919_232145_s1.hifi_reads.bam &

$ pbindex m84010_220919_235306_s2.hifi_reads.bam &
$ bam2fastq -o m84010_220919_235306_s2.hifi_reads m84010_220919_235306_s2.hifi_reads.bam &

$ pbindex m84011_220902_175841_s1.hifi_reads.bam &
$ bam2fastq -o m84011_220902_175841_s1.hifi_reads m84011_220902_175841_s1.hifi_reads.bam &

將fastq.gz 文件重新命名:

# bam2fastq程序運(yùn)行完得到以下文件
$ ls
m84010_220919_232145_s1.hifi_reads.fastq.gz，m84010_220919_235306_s2.hifi_reads.fastq.gz，m84011_220902_175841_s1.hifi_reads.fastq.gz

#修改名字
mv m84010_220919_232145_s1.hifi_reads.fastq.gz HG004.fastq.gz
mv m84010_220919_235306_s2.hifi_reads.fastq.gz HG003.fastq.gz
mv m84011_220902_175841_s1.hifi_reads.fastq.gz HG002_1.fastq.gz

• 運(yùn)行minimap2

#數(shù)據(jù)是PacBio-HiFi-CCS數(shù)據(jù)
$ minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or later)

#實(shí)際運(yùn)行
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG002_1.fastq.gz | samtools view -@ 12 -bS | samtools sort -@ 12 -o HG002_1.minimap2.align.bam
#或者可以用以下命令一步完成sam到bam，排序和索引
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG002_1.fastq.gz | samtools sort -@ 12 -o HG002_1.minimap2.align.bam --write-index -

#最終版本, samtools=1.18
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG002_1.fastq.gz | samtools sort -@ 12 -o HG002_1.minimap2.align.bam --write-index -
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG003.fastq.gz | samtools sort -@ 12 -o HG003.minimap2.align.bam --write-index -
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG004.fastq.gz | samtools sort -@ 12 -o HG004.minimap2.align.bam --write-index -

#最終版本， samtools=1.9
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG002_1.fastq.gz | samtools view -@ 12 -bS | samtools sort -@ 12 -o HG002_1.minimap2.align.bam
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG003.fastq.gz | samtools view -@ 12 -bS | samtools sort -@ 12 -o HG003.minimap2.align.bam
$ minimap2 --MD -t 12 -ax map-hifi Human_ref/GRCh38.fa HG004.fastq.gz | samtools view -@ 12 -bS | samtools sort -@ 12 -o HG004.minimap2.align.bam
#因?yàn)閟amtools=1.9沒(méi)有sort沒(méi)有--write-index選項(xiàng)

對(duì)于minimap2的參數(shù)：
-a Generate CIGAR and output alignments in the SAM format. Minimap2 outputs in PAF by default.
--MD Output the MD tag (see the SAM spec).（后面sniffles軟件需要MD tag）.
-t Number of threads CPU線程數(shù).
-x STR preset (always applied before other options; see minimap2.1 for details) .
因?yàn)閙inimap2輸出的是.sam文件格式，所以使用samtools將.sam轉(zhuǎn)換成.bam，并且使用samtools sort 對(duì) .bam 進(jìn)行排序。下面是samtools的參數(shù)：
-@ 線程數(shù).
-b output BAM 默認(rèn)下輸出是 SAM 格式文件，該參數(shù)設(shè)置輸出 BAM 格式.
-S input is SAM 默認(rèn)下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數(shù)，否則有時(shí)候會(huì)報(bào)錯(cuò).

二、Sniffles2 使用

官方網(wǎng)站：https://github.com/fritzsedlazeck/Sniffles

sniffles2流程圖

1. 軟件安裝

# pip 安裝
$ pip install sniffles

# conda 安裝
$ conda install sniffles=2.2

2. 軟件使用

sniffles2使用分為四種場(chǎng)景：

1. 單樣本鑒定結(jié)構(gòu)變異

• 單樣本.bam文件

$ sniffles --input sorted_indexed_alignments.bam --vcf output.vcf

$ sniffles --input sorted_indexed_alignments.bam --vcf output.vcf.gz

• 單樣本.cram文件

$ sniffles --input sample.cram --vcf output.vcf.gz

• 單樣本生成單個(gè).snf文件，后期用于多樣本鑒定結(jié)構(gòu)變異

$ sniffles --input sample1.bam --snf sample1.snf

• 同時(shí)生成.vcf和.snf文件，.snf后期用于多樣本鑒定結(jié)構(gòu)變異

$ sniffles --input sample1.bam --vcf sample1.vcf.gz --snf sample1.snf

• 指定串聯(lián)重復(fù)區(qū)域以及參考基因組序列。指定串聯(lián)重復(fù)區(qū)域能提高這一區(qū)域的檢測(cè)性能。

$ sniffles --input sample1.bam --vcf sample1.vcf.gz --tandem-repeats tandem_repeats.bed --reference genome.fa --mosaic

2. Mosaic SV Calling (Non-germline or somatic SVs) 馬賽克模式結(jié)構(gòu)變異檢測(cè)，適用于非胚系或者體細(xì)胞結(jié)構(gòu)變異

#只需加入--mosaic
$ sniffles --input mapped_input.bam --vcf output.vcf --mosaic

3. 多樣本聯(lián)合變異

#Step 1. Create .snf for each sample: 
$ sniffles --input sample1.bam --snf sample1.snf

#Step 2. Combined calling:
$ sniffles --input sample1.snf sample2.snf ... sampleN.snf --vcf multisample.vcf

4. 對(duì)指定變異進(jìn)行檢測(cè)

$ sniffles --input sample.bam --genotype-vcf input_known_svs.vcf --vcf output_genotypes.vcf

實(shí)際運(yùn)行

#單樣本分別檢測(cè)變異
$ sniffles -t 12 --input HG002_1.minimap2.align.bam --vcf HG002_1.vcf.gz --snf HG002_1.snf --tandem-repeats Human_ref/human_GRCh38_no_alt_analysis_set.trf.bed
$ sniffles -t 12 --input HG003.minimap2.align.bam --vcf HG003.vcf.gz --snf HG003.snf --tandem-repeats Human_ref/human_GRCh38_no_alt_analysis_set.trf.bed
$ sniffles -t 12 --input HG004.minimap2.align.bam --vcf HG004.vcf.gz --snf HG004.snf --tandem-repeats Human_ref/human_GRCh38_no_alt_analysis_set.trf.bed

#joint calling
$ sniffles --input HG002_1.snf  HG003.snf HG004.snf --vcf multisample.vcf.gz

整個(gè)joint calling的時(shí)間很短，不到20秒。

Joint calling整個(gè)運(yùn)行過(guò)程

3. 最終結(jié)果

最后得到整合的vcf文件 multisample.vcf.gz。

整合vcf部分結(jié)果

參考文獻(xiàn)：

1. 生信分析|Minimap2+sniffles calling SVs