samtools 一些使用命令

2024-05-30  本文已影响0人  风知秋

对 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

上一篇下一篇

猜你喜欢

热点阅读