对 SAM 或 BAM 格式的文件进行快速的统计分析:
samtools flagstat xx.bam
示例:
415088184 + 0 in total (QC-passed reads + QC-failed reads)
55152890 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
408418313 + 0 mapped (98.39% : N/A)
359935294 + 0 paired in sequencing
179967647 + 0 read1
179967647 + 0 read2
340743416 + 0 properly paired (94.67% : N/A)
352376046 + 0 with itself and mate mapped
889377 + 0 singletons (0.25% : N/A)
9397130 + 0 with mate mapped to a different chr
4761549 + 0 with mate mapped to a different chr (mapQ>=5)
示例说明:
总共有 415,088,184 条 reads,其中全部通过质量控制(QC-passed reads)。
没有任何 QC 失败的 reads。
共有 55,152,890 条次要对(secondary alignments)。
没有补充对(supplementary alignments)。
没有重复 reads。
共有 408,418,313 条 reads 被成功对齐(mapped),占总 reads 数的 98.39%。
共有 359,935,294 条 reads 是成对测序的。
其中有 179,967,647 条是 read1,179,967,647 条是 read2。
共有 340,743,416 条 reads 被正确配对(properly paired),占成对 reads 的 94.67%。
有 352,376,046 条 reads 自身和其 mate 均被成功对齐。
共有 889,377 条单一(singleton) reads,占总 reads 的 0.25%。
有 9,397,130 条 reads 的 mate 被成功对齐到了不同的染色体。
有 4,761,549 条 reads 的 mate 被成功对齐到了不同的染色体且 mapQ 值大于等于 5。
查看 bam 文件的头信息
samtools view -H xx.bam
包括文件头 (@HD)、染色体信息 (@SQ)、读组信息 (@RG)、以及处理 BAM 文件的程序信息 (@PG)。
@HD: 这个记录指定了 BAM 文件的版本 (VN)、排序方式 (SO) 等信息。
@SQ: 每个@SQ记录描述了一个染色体的信息,包括染色体名字 (SN) 和长度 (LN)。
@RG: 每个@RG记录描述了一个读组的信息,包括唯一标识符 (ID)、样本名称 (SM)、测序平台 (PL)、文库名称 (LB)、测序单元 (PU) 等。
@PG: 每个@PG记录描述了一个处理程序的信息,包括程序名字 (PN)、版本 (VN)、命令行参数 (CL) 等。
每一个条目之间必须制表符分隔!
报错:
A USER ERROR has occurred: Argument --emit-ref-confidence has a bad value: Can only be used in single sample mode currently. Use the --sample-name argument to run on a single sample out of a multi-sample BAM file.
原因 bam 文件中头信息有问题:
@RG ID:TRJ_NamRoo_1 SM:TRJ_NamRoo_1 PL:ILLUMINA LB:TRJ_NamRoo_1 PU:L
@RG ID:TRJ_NamRoo_2 SM:TRJ_NamRoo_2 PL:ILLUMINA LB:TRJ_NamRoo_2 PU:L
@RG ID:TRJ_NamRoo_3 SM:TRJ_NamRoo_3 PL:ILLUMINA LB:TRJ_NamRoo_3 PU:L
@RG ID:TRJ_NamRoo_4 SM:TRJ_NamRoo_4 PL:ILLUMINA LB:TRJ_NamRoo_4 PU:L
SM:一列的样本信息必须是唯一指对,这里将每一行均改为 TRJ_NamRoo 即可。
可以通过 samtools reheader new.header input.bam > output.bam 解决。
samtools 进行 sam 转 bam,比 gatk 运行速度快:
samtools view -b input.sam -o input.bam
多线程运行:
samtools view -@ 4 -b input.sam -o input.bam
samtools 截取特定染色体:
samtools view -b input.bam chr1 -o chr1.bam
samtools 截取参考基因组中的特定序列:
samtools faidx ref.fa chr1:100-500 > test.fa