专栏首页生物信息云转录组分析 | 使用SAMtools将SAM文件转换为BAM文件、排序、建立索引

转录组分析 | 使用SAMtools将SAM文件转换为BAM文件、排序、建立索引

我们前面经过质控,比对得到sam文件:

转录组分析 | fastqc进行质控与结果解读

转录组分析 | 使用trim-galore去除低质量的reads和adaptor

转录组分析 | 使用Hisat2进行序列比对

接下来为什么要进行格式转换?为了让计算机好处理。SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式。bam是sam的二进制格式,减少了sam文件的储存量。这在我之前的文章【生信中常见的数据文件格式】中也有介绍。

接下来,我们要做的事情就是使用SAMtools将SAM文件转换为BAM文件、排序、建立索引。

一.SAMtools介绍

SAMtools是一个用于操作sam和bam文件的工具合集。能够实现二进制查看、格式转换、排序及合并等功能,结合sam格式中的flag、tag等信息,还可以完成比对结果的统计汇总。同时利用linux中的grep、awk等操作命令,还可以大大扩展samtools的使用范围与功能。包含有许多命令。我这里主要介绍几个,重点是samtools view。

samtools更多信息:http://www.htslib.org/

1.samtools view

samtools view主要用来转换SAM/BAM/CRAM文件。将sam文件与bam文件互换;然后对bam文件进行各种操作,比如数据的排序(sort)和提取(这些操作 是对bam文件进行的,因而当输入为sam文件的时候,不能进行该操作);最后将排序或提取得到的数据输出为bam或sam(默认的)格式。如果没有指定选项或区域,则将指定的输入对齐文件(SAM、BAM或CRAM格式)中的所有对齐打印到SAM格式的标准输出(没有标头)。可以在输入文件名后指定一个或多个空格分隔的区域规范,以将输出限制为仅覆盖指定区域的那些对齐。使用区域规范需要一个协调排序和索引的输入文件(BAM或CRAM格式)。samtools view的参数很多,-b、-C、-1、-u、-h、-H和-C选项将更改缺省的无header SAM的输出格式,而-o和-u选项将设置输出文件名。-t和-T选项提供了额外的参考数据。当SAM输入不包含@SQ headers时,这两个选项中的一个是必需的,当编写CRAM输出时,-T选项是必需的。-L、-M、-r、-R、-d、-D、-s、-q、-L、-M、-f、-F和-G选项过滤将包含在输出中的对齐,只筛选那些匹配特定条件的对齐。-x和-B选项修改包含在每次对齐中的数据。如果数据文件夹不包含任何索引文件,可以使用-X选项允许用户指定定制的索引文件位置。

samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

参数介绍:

# samtools view 使用说明:
 -b output BAM
  默认下输出是 SAM 格式文件,该参数设置输出 BAM 格式
 -h print header for the SAM output
  默认下输出的 sam 格式文件不带 header,该参数设定输出sam文件时带 header 信息
 -H print header only (no alignments)
  仅仅输出文件的头
 -S input is SAM
  默认下输入是 BAM 文件,若是输入是 SAM 文件,则最好加该参数,否则有时候会报错。
 -u uncompressed BAM output (force -b)
  该参数的使用需要有-b参数,能节约时间,但是需要更多磁盘空间。
 -c Instead of printing the alignments, only count them and print the 
  total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ ,   are taken into account.
  过滤功统计功能
 -c print only the count of matching records
 -L FILE  output alignments overlapping the input BED FILE [null]
 -t FILE  list of reference names and lengths (force -S) [null]
  使用一个list文件来作为header的输入
 -T FILE  reference sequence file (force -S) [null]
  使用序列fasta文件作为header的输入
 -o FILE  output file name [stdout]
 -F INT   filtering flag, 0 for unset [0] 
  Skip alignments with bits present in INT [0]
  数字4代表该序列没有比对到参考序列上
  数字8代表该序列的mate序列没有比对到参考序列上
  过滤功能。如F12过滤只有双端map的
 -q INT   minimum mapping quality [0]
    比对的最低质量值,一般认为20就为unique比对了,可以结合上述-bF参数使用使用提取特定的比对结果

其他更多参数介绍,阅读:http://www.htslib.org/doc/samtools-view.html

2.samtools sort

samtools sort用来对SAM/BAM/CRAM文件进行排序,按最左坐标排序,或使用-n时按读取名称排序。默认情况下,排序后的输出被写到标准输出,或者在使用-o时写到指定的文件(out.bam)。此命令还将创建临时文件tmpprefixv .%d。当整个对齐数据无法装入内存时(通过-m选项控制),根据需要使用bam。

samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] [-T tmpprefix] [-@ threads] [in.sam|in.bam|in.cram]

重要参数:

-l INT 设置输出文件压缩等级。0-9,0是不压缩,9是压缩等级最高。不设置此参数时,使用默认压缩等级;
-m INT 设置每个线程运行时的内存大小,可以使用K,M和G表示内存大小。
-n 设置按照read名称进行排序;
-o FILE 设置最终排序后的输出文件名;
-O FORMAT 设置最终输出的文件格式,可以是bam,sam或者cram,默认为bam;
-T PREFIX 设置临时文件的前缀;
-@ INT 设置排序和压缩是的线程数量,默认是单线程。

3.samtools index

必须对bam文件进行默认情况下的排序后,才能进行index。否则会报错。建立索引后将产生后缀为.bai的文件,用于快速的随机处理。很多情况下需要有bai文件的存在,特别是显示序列比对情况下。比如samtool的tview命令就需要;gbrowse2显示reads的比对图形的时候也需要。IGV显示比对情况也需要。

samtools index [-bc] [-m INT] aln.sam|aln.bam|aln.cram [out.index]

参数:

-b 创建一个BAI索引。当不使用格式选项时,这是当前的默认设置。
-c 创建CSI索引。默认情况下,索引的最小间隔大小为2^14,与BAI格式使用的固定值相同。
-m INT 创建CSI索引,最小间隔大小为2^INT。
-@, --threads INT
    除了主线程[0]之外,要使用的输入/输出压缩线程数。

samtools sort命令按默认染色体位置排序,顺利建立Index,如果前面排序有出入,可能不能正确建立索引。

参考:http://www.htslib.org/doc/samtools-index.html

4.samtools flagstat

samtools flagstat用于给出BAM文件的比对结果。

samtools flagstat in.sam|in.bam|in.cram

参数:

-@ INT 设置读取文件时要使用的额外线程数。
-O FORMAT 设置输出格式。FORMAT可以设置为'default', 'json'或'tsv'来选择默认的,json或标签分隔值输出格式。如果不使用此选项,将选择默认格式。

参考:http://www.htslib.org/doc/samtools-flagstat.html


关于SAMtools工具的用法还很多,我就介绍上面几个,其他更多用法,阅读:http://www.htslib.org/doc/samtools.html

参考:https://blog.csdn.net/sinat_38163598/article/details/72910115

二.SAMtools实战

1.格式转换(samtools view)

我们先查看一下上一篇文章【转录组分析 | 使用Hisat2进行序列比对】中产生的sam文件。

ll -h ./cleandata/hisat2_mm10data/

先建立一个输出数据的文件夹。

mkdir ./cleandata/samtools_bam

同样的,我们可以单个转换,也可以写一个脚本批量转换,首先,我先介绍单个转换。

samtools view -S ./cleandata/hisat2_mm10data/CK4.sam -b > ./cleandata/samtools_bam/CK4.bam

这个过程的时间也是挺长的。

批量转换,创建脚本文件

 vim samtools_view_batch.sh

输入下面内容:

for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
  samtools view -S ./cleandata/hisat2_mm10data/${i}.sam -b -o ./cleandata/samtools_bam/${i}.bam
done

保存,然后执行脚本。

bash samtools_view_batch.sh

运行结束后我么你查看一下输入的数据:

ll -h ./cleandata/samtools_bam/

我们可以看到,转换后的数据只有32G,原来可是155G。相差5倍左右。

2 排序(samtools sort)

单个排序:

samtools sort -l 4 -o ./cleandata/samtools_bam/CK-4_sort.bam ./cleandata/samtools_bam/CK4.bam 

批量排序:创建脚本文件

vim samtools_sort.sh

输入下面内容:

#!/bin/bash
# This is for samtools_sort
for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
samtools sort -l 4 -o ./cleandata/samtools_bam/${i}_sort.bam ./cleandata/samtools_bam/${i}.bam
done

保存并执行脚本。

 bash samtools_sort.sh

执行结束后,我们查看一下文件:

ll -h cleandata/samtools_bam/*sort.bam

sort后文件几乎小了2倍。

3. 建立索引(samtools index)

samtools sort命令时,按默认染色体位置排序,顺利建立Index,如果前面排序有出入,可能不能正确建立索引。

这里我就一次建立索引了。

for i in CK-4 CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9;do samtools index ./cleandata/samtools_bam/${i}_sort.bam ./cleandata/samtools_bam/${i}_sort.bam.bai;done

4.查看reads比对情况(samtools flagstat)

samtools flagstat cleandata/samtools_bam/CK-4_sort.bam
(RNAseq1) [root@localhost RNAseq]# samtools flagstat cleandata/samtools_bam/CK-4_sort.bam
51977802 + 0 in total (QC-passed reads + QC-failed reads)
#总共的reads数
4063860 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
50969572 + 0 mapped (98.06% : N/A)
#总体上reads的匹配率
47913942 + 0 paired in sequencing
#有多少reads是属于paired reads
23956971 + 0 read1 
#reads1中的reads数
23956971 + 0 read2
#reads2中的reads数
45980036 + 0 properly paired (95.96% : N/A)
#完美匹配的reads数:比对到同一条参考序列,并且两条reads之间的距离符合设置的阈值
46308722 + 0 with itself and mate mapped
#paired reads中两条都比对到参考序列上的reads数
596990 + 0 singletons (1.25% : N/A)
#单独一条匹配到参考序列上的reads数,和上一个相加,则是总的匹配上的reads数。
232008 + 0 with mate mapped to a different chr
#paired reads中两条分别比对到两条不同的参考序列的reads数
210169 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一个,只是其中比对质量>=5的reads的数量

参考:

【1】https://blog.csdn.net/sinat_38163598/article/details/72910115

【2】http://www.htslib.org/doc/samtools.html

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-09-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Linux系统下Anaconda的安装和使用教程

    去官网下载:https://www.anaconda.com/products/individual

    DoubleHelix
  • PanGPCR | 预测多个潜在的GPCR靶标及其在组织中的表达位置,副作用以及GPCR药物的可能用途

    靶向G蛋白偶联受体(GPCR)(已知的最大治疗靶标)的药物发现具有挑战性。为了促进GPCR药物的快速发现和开发,Yufeng J Tseng等人构建了PanGP...

    DoubleHelix
  • 转录组分析 | 使用RSeQC软件对生成的BAM文件进行质控

    RSeQC是发表于2012年的一个RNA-Seq质控工具,属于python包。它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据,比如一些基本...

    DoubleHelix
  • Android 混合开发之JsBridge

    电商或者内容类APP中,H5通常都会占据一席之地,Native跟H5通信会必不可少,比如某些场景H5通知native去分享,native通知H5局部刷新等,An...

    看书的小蜗牛
  • samtools`markdup`操作的正确顺序

    具体实例在5 一步法找基因变异流程 samtoolsmarkdup操作的正确顺序

    Y大宽
  • 10行代码搞定【滚动回归】

    对于任意一天t,在[t - n, t]的区间内进行回归。如果数据一共有N天,那么就会得到N - n个数据点

    用户7652506
  • ESA2GJK1DH1K升级篇: 引入网页实现MQTT控制- 网页实现MQTT通信入门

      所以咱们会在数据篇或者安全篇做一套网页的管理软件,用来管理咱所有的MQTT设备

    杨奉武
  • ansible简易入门之playbook

    jeremyxu
  • Webpack实战-构建 Electron 应用

    Electron 可以让你使用开发 Web 的技术去开发跨平台的桌面端应用,由 Github 主导和开源,大家熟悉的 Atom 和 VSCode 编辑器就是使用...

    IMWeb前端团队
  • 总结了一些指针易出错的常见问题(五)

    指针与链表及其操作 //结构体定义 typedef struct _person{ char* firstname; char* lastna...

    互联网金融打杂

扫码关注云+社区

领取腾讯云代金券