全基因组 - 人类基因组变异分析(PacBio) -- minimap2 + Sniffles2
原创全基因组 - 人类基因组变异分析(PacBio) -- minimap2 + Sniffles2
原创
一、Minimap2 使用
官方网站：https://github.com/lh3/minimap2
1. 软件安装
首先从github官网上下载minimap2的二进制文件压缩包，minimap2-2.26_x64-linux.tar.bz2，然后上传到服务器上。
# minimap2，v2.26压缩包解压缩
$ tar -xjvf minimap2-2.26_x64-linux.tar.bz2
# -x 解压
# -j 有bz2属性的
# -v 显示所有过程
# -f 使用档案名字，切记，这个参数是最后一个参数，后面只能接档案名。
#将minimap2加入绝对路径
$ 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 ~/.bashrc2. 软件使用
- 数据格式转化
公共数据只有 .bam 格式，所以要先将.bam 文件转化为.fastq文件才能输入minimap2；如果直接获得.fastq文件则可以省略此步转化。
数据格式转化所用的程序为BAM2fastx（PacBio官方工具），PacBio将一系列工具，包括对.bam文件进行索引的pbindex，都放在pbtk（pb tool kit）中，所以运行以下命令全部安装：
BAM2fastx tools ：https://github.com/PacificBiosciences/bam2fastx
# conda安装pbtk
$ conda install -c bioconda pbtk德系犹太人家系：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
对.bam文件进行索引（生成.bam.pbi文件），然后进行.bam 文件到.fastq文件的转换：
# generates out.fastq.gz
$ bam2fastq -o out in_1.bam in_2.bam in_3.xml in_4.bam
#因为运行bam2fastq 需要对bam文件先进行索引
$ 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程序运行完得到以下文件
$ 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- 运行minimap2
#数据是PacBio-HiFi-CCS数据
$ minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or later)
#实际运行
$ 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
#因为samtools=1.9没有sort没有--write-index选项对于minimap2的参数：
-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线程数.
-x STR preset (always applied before other options; see minimap2.1 for details) .
因为minimap2输出的是**.sam**文件格式，所以使用**samtools**将**.sam**转换成**.bam**，并且使用**samtools sort** 对 .bam 进行排序。下面是samtools的参数：
-@ 线程数.
-b output BAM 默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式.
-S input is SAM 默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错.
二、Sniffles2 使用
官方网站：https://github.com/fritzsedlazeck/Sniffles
1. 软件安装
# pip 安装
$ pip install sniffles
# conda 安装
$ conda install sniffles=2.22. 软件使用
sniffles2使用分为四种场景：
1. 单样本鉴定结构变异
- 单样本
.bam文件
#vcf
$ sniffles --input sorted_indexed_alignments.bam --vcf output.vcf
#vcf.gz
$ sniffles --input sorted_indexed_alignments.bam --vcf output.vcf.gz- 单样本
.cram文件
$ sniffles --input sample.cram --vcf output.vcf.gz- 单样本生成单个
.snf文件，后期用于多样本鉴定结构变异
$ sniffles --input sample1.bam --snf sample1.snf- 同时生成
.vcf和.snf文件，.snf后期用于多样本鉴定结构变异
$ sniffles --input sample1.bam --vcf sample1.vcf.gz --snf sample1.snf- 指定串联重复区域以及参考基因组序列。指定串联重复区域能提高这一区域的检测性能。
$ sniffles --input sample1.bam --vcf sample1.vcf.gz --tandem-repeats tandem_repeats.bed --reference genome.fa --mosaic2. Mosaic SV Calling (Non-germline or somatic SVs) 马赛克模式结构变异检测，适用于非胚系或者体细胞结构变异
#只需加入--mosaic
$ sniffles --input mapped_input.bam --vcf output.vcf --mosaic3. 多样本联合变异
#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.vcf4. 对指定变异进行检测
$ sniffles --input sample.bam --genotype-vcf input_known_svs.vcf --vcf output_genotypes.vcf实际运行
#单样本分别检测变异
$ 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整个joint calling的时间很短，不到20秒。
3. 最终结果
最后得到整合的vcf文件 multisample.vcf.gz。
参考文献：
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。