Samtools对比对结果的处理

原创已于 2025-07-22 11:29:11 修改 · 1.5k 阅读

26 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#linux #服务器 #生物信息学

于 2025-07-22 10:49:59 首次发布

汇佰生物-生信介绍专栏收录该内容

9 篇文章

订阅专栏

该文章已生成可运行项目，

常用序列比对软件的区别和选择(获取sam文件)-CSDN博客https://blog.csdn.net/Avalon96/article/details/149324810?spm=1001.2014.3001.5501上一篇我们介绍了如何通过比对获取sam文件，接下来我们会介绍如何对sam文件进行后续处理。

1、Sam文件格式详解

头部行以'@'开始，后面跟着两个字母的代码标识头部的类型。常见的头部类型包括：

@HD：文件级别的元数据，如版本、排序顺序。
@SQ：参考序列字典，包括参考序列名称（SN）和长度（LN）。
@RG：读段（read）所属的测序文库信息。
@PG：使用的程序信息。
@CO：注释信息。

比对部分每行有11个必选字段，字段之间用制表符分隔。这些字段按顺序分别是：

QNAME：查询序列（即测序读段）的名称。
FLAG：一个整数，表示比对的属性（如是否配对、是否比对上等），通过位掩码表示。
RNAME：参考序列的名称（即比对到的染色体或contig名称），如果未比对则为'*'。
POS：比对上的最左侧位置（从1开始计数），如果未比对则为0。
MAPQ：比对质量（Mapping Quality），表示该比对正确的概率的Phred分数。
CIGAR：一个字符串，表示比对的细节（如匹配、插入、缺失等）。
RNEXT：配对读段中另一端的参考序列名称，如果没有则为'*'，如果是同一个参考序列则为'='。
PNEXT：配对读段中另一端的比对位置。
TLEN：模板长度（Template length），即插入片段长度。
SEQ：查询序列的碱基序列，如果未存储则为'*'（但通常都会存储）。
QUAL：查询序列的质量值，使用ASCII编码的Phred质量分数。

2、Samtools的安装

samtools/samtools: Tools (written in C using htslib) for manipulating next-generation sequencing datahttps://github.com/samtools/samtools推荐在conda中安装Samtools，安装方法如下

conda install samtools -c bioconda

samtools help #查看安装是否成功

3、Samtools基本命令

view sam与bam格式转换

此命令可以将sam格式转为bam格式，bam格式的优点在于二进制文件占用磁盘空间小，参与后续分析运行速度快。

同时在这步还可以实现一些筛选的操作，比如可以提取比对/未比对到参考基因组的序列，到bam文件中。

# sam格式转bam格式
samtools view -bS sample.sam  -o sample.bam

# bam格式转sam格式
samtools view -h sample.bam -o sample.sam

# 提取比对到参考基因组上的序列
samtools view -bF 4 sample.bam  -o sample.map.bam

# 提取没有比对到参考基因组上的序列
samtools view -bf 4 sample.bam  -o sample.unmap.bam

# 提取两条reads都比对到参考基因组上的序列
samtools view -bF 12 sample.bam  -o sample.paired.bam

# 提取比对到Chr1染色体的序列
samtools view sample.bam chr1 > scaffold1.sam

sort 对bam文件进行排序

上一步生成的bam文件一般需要排序后才可使用

samtools sort sample.sam -o sample_sort.bam

index bam文件建索引

此命令对排序后的bam文件建立索引，获得bai文件。bai经常和bam文件绑定使用，所以也是排序后需要运行的一个必须步骤。

# 此命令将自动生成索引文件 sample_sort.bam.bai
samtools index sample_sort.bam

faidx Fasta文件建索引

此命令将根据读取的Fasta文件建立索引，获得fai文件。常用于对参考基因组建索引。

samtools faidx genome.fasta

flagstat 查看比对统计结果

以现有bam文件为例，展示统计结果如下

samtools flagstat sample_sort.bam

7199749 + 0 in total (QC-passed reads + QC-failed reads)

7199454 + 0 primary ##总reads数

295 + 0 secondary

0 + 0 supplementary

0 + 0 duplicates

0 + 0 primary duplicates

6536962 + 0 mapped (90.79% : N/A)

6536667 + 0 primary mapped (90.79% : N/A) ##比对率

7046435 + 0 paired in sequencing

3523084 + 0 read1

3523351 + 0 read2

6404378 + 0 properly paired (90.89% : N/A)

6409438 + 0 with itself and mate mapped

4615 + 0 singletons (0.07% : N/A)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

depth 比对深度统计

统计参考基因组上每一个位点的测序深度，获得一个位点-深度列表，可以用于后续可视化绘图。

samtools depth sample.bam > depth.txt

4、使用示例——新冠数据比对

现有新冠测序下机数据，sample_R1.fq.gz、sample_R2.fq.gz。

新冠参考基因组，MN908947.3.fasta。

人参考基因组，hg38.fasta。

以下是分析流程以及步骤说明：

序列质量

1、查看原始数据质量，结果输出到raw_fastqc_result目录中。

fastqc sample_R1.fq.gz sample_R2.fq.gz -o raw_fastqc_result

2、序列质控，以8线程运行fastp程序，去除质量不好的序列。

fastp -w 8 -i sample_R1.fq.gz -I sample_R2.fq.gz -o fastp_output.clean.fastq  -O fastp_output.clean.2.fastq

3、查看质控后序列质量，结果输出到clean_fastqc_result目录中。

fastqc fastp_output.clean.fastq fastp_output.clean.2.fastq -o clean_fastqc_result

去宿主

以下去宿主步骤为通用步骤，兼容任何可以得到sam文件的比对软件。如果只想用bowtie2去宿主，可以使用bowtie2自带的去宿主命令：

bowtie2 -x human -1 reads_R1.fastq -2 reads_R2.fastq -p 8 -S hg.sam --un-conc reads_host_removed

4、人参考序列建索引

bowtie2-build hg38.fasta hg38

5、将序列比对到人参考基因组上。

bowtie2 -p 8 -x hg38 -1 fastp_output.clean.fastq -2 fastp_output.clean.2.fastq -S hg38_output.sam

6、使用samtools将sam转bam格式并排序

samtools view -Sb hg38_output.sam -o hg38_output.bam
samtools sort hg38_output.bam -o hg38_output_sort.bam

7、使用samtools提取未比对的序列，并转为fastq格式

samtools view -b -f 4 hg38_output_sort.bam > unmapped.bam
samtools bam2fq -1 unmapped_1.fastq -2 unmapped_2.fastq unmapped.bam

如果后续分析中，出现双端序列不匹配的问题可以使用bbmap中的序列修复处理。

# bbmap安装 conda install bbmap -c bioconda
repair.sh in=unmapped_1.fastq in2=unmapped_2.fastq out=unmapped_1_repaired.fastq out2=unmapped_2_repaired.fastq

参考基因组比对

8、新冠参考基因组建索引

bowtie2-build MN908947.3.fasta sars-cov-2

9、将去宿主后的序列比对到新冠参考基因组上

bowtie2 -p 8 -x sars-cov-2 -1 unmapped_1_repaired.fastq  -2 unmapped_2_repaired.fastq  -S output.sam

10、使用samtools将sam转bam格式并排序

samtools view -Sb output.sam -o output.bam
samtools sort output.bam -o output_sort.bam

11、使用samtools查看比对情况

# 比对率计算
samtools flagstat output_sort.bam

# 深度计算
samtools depth output_sort.bam -o depth.txt

病原微生物分析系统——新冠模块

由深圳臻合智造生物科技有限公司开发的病原微生物分析系统，在新冠分析模块对于比对率和深度有直观的数据统计展示。在后续的变异和进化分析中也进行了可视化统计。

深度可视化

比对覆盖率可视化

变异统计可视化

进化树可视化

病原微生物分析系统

本文章已经生成可运行项目