BBQ(生信基础问题17)----mapping专题:SAM和BAM了解一下?

我们利用bowtie建立了index的话,接下来是不是可以直接进行比较了呢?那么mapping(是指把数据回贴到参考基因组的这个过程)是如何实现的呢?

----------------------------------------------我是给出代码的分割线-----------------------------------------

cd /asnas/sunyl_group/liull/snptest/ETYY_1F/bwa

bwa mem  -t 5 -R "@RG\tID:L2\tSM:L2\tLB:WGS\tPL:Illumina" /asnas/sunyl_group/liull/Database/hg38/chroms/hg38.fa ./ETYY_1F_D17065907_H2N2HALXX_L2_1_paired.fq  ./ETYY_1F_D17065907_H2N2HALXX_L2_2_paired.fq > ETYY_1F_D17065907_H2N2HALXX_L2_paired.sam
samtools sort ./ETYY_1F_D17065907_H2N2HALXX_L1_paired.sam > ./ETYY_1F_D17065907_H2N2HALXX_L1_paired_sorted.bam

1:先cd进入到我的fq的文件夹中 2:利用的是bwa 这个软件,mem模式(这个模式是适用于)70bp-1Mbp的序列比对的,适用于二代测序的结果 3:-t 是规定了线程数 4:-R是给出了头的信息

-------------------------------------分割线------------------------------------------------------ 我们来仔细看一下代码。我们发现生成文件是2种类型的:bam/sam 问题来了,什么是Bam文件,什么是Sam文件呢? 今天我们就要谈一谈SAM/BAM文件格式。 SAM(The Sequence Alignment / Map format)格式,即序列比对文件的格式,详细介绍文档:http://samtools.github.io/hts-specs/SAMv1.pdf

首先先说二者之间的关系,BAM文件是SAM文件的压缩格式,压缩以后可以节省空间,排好序的BAM文件还可以提供随机访问功能,性能优良。但是BAM文件和SAM文件储存的内容是完全一样的。我们以后还要单独再说BAM文件的操作方法,今天我们把重点放在文件中的内容上。

SAM文件的全称是:Sequence Alignment Map,它设计之初就是为了存储mapping结果的。一个标准的SAM文件由2部分组成,第1部分是以“@”开头的头部,在文件的最前面;第2部分就是紧跟在头部后面的比对结果文件。 我随便打开了一个我的Sam文件:可以看到sam文件是以@开头的,记录了比对的信息结果。

sam文件head

然后我们再接着往下翻,看看红框的部分,我们可以执行的比对的代码。问题来了,这个多列,每一列的信息是什么呢?

sam文件的“正文”

BBQ问题:1:SAM文件的头部内容中常见的标志符号有@HD,@SQ,@PG,请问这三者后面跟随的信息分别是什么意思?

@HD VN:1.0 SO:unsorted (排序类型) 头部区第一行:VN是格式版本;SO表示比对排序的类型,有unknown(default),unsorted,queryname和coordinate几种。 @SQ SN:1 LN:249250621 (序列ID及长度) 参考序列名,这些参考序列决定了比对结果sort的顺序,SN是参考序列名;LN是参考序列长度;每个参考序列为一行。 例如:@SQ SN:NC_000067.6 LN:195471971 @RG ID:sample01 (样品基本信息) Read Group。1个sample的测序结果为1个Read Group;该sample可以有多个library的测序结果,可以利用bwa mem -R 加上去这些信息。

@RG信息

ID:样品的ID号 SM:样品名 LB:文库名 PU:测序以 PL:测序平台 这些信息可以在形成sam文件时加入,ID是必须要有的后面是否添加看分析要求 @PG 写清楚了测序的平台,名字以及命令等。

@PG信息

这里的ID是bwa,PN是bwa,VN是0.7.10-r789版本。CL可以认为是运行程序@RG是上面RG表示的内容,后面是程序内容,这里的@GR内容是可以自己在运行程序是加入的 ---------------------------------------------分割线---------------------------------------------------------

标准SAM文件中的11列内容代表的含义

我们主要关注前面4列: 第一列:比对片段的(template)的编号 第二列:FLAG:标识信息,这个信息量很大

官方文档给出的FLAG解释

举例说明,比如说实例中的99=64+32+2+1, 也就是这个记录所代表的read是来自于双端测序R1,且匹配到了正链。而147=128+16+2+1则是表示这个记录来自于双端测序的R2,完全匹配到负链。 有一个网站可以直接根据数字计算匹配信息:https://www.plob.org/article/1697.html 第三列: 参考序列的编号,匹配到的染色体信息,比对到了第几号染色体 第四列:比对上的位置,注意是从1开始计数。

BBQ问题: 如果1条序列的FLAG=83 (图3标号38的行)请解释其比对含义。

FLAG结果

解释

--------------------------------------------分割线-------------------------------------- 附:

给出BWA的部分帮助文档信息:

1: Index  

bwa   index   [ - p   p r e f i x ]   [ - a   a l g o T y p e ]   < i n . d b . f a s t a >
index数据库序列以FASTA格式。

选项

-p STR   输出数据库的前缀[与db 文件名相同]

-a STR   算法用于构建BWT  index。可以使用的选项:

      is   IS线性时间算法用于构建suffix array。需要5.37N内存,N是数据库的大小。IS算法相对较快,但是无法处理数据库大于2GB的数据。因为IS算法比较简单,作为默认值。目前IS算法的脚本由Yuta Mori从新植入。

      Bwtsw  BWT-SW中使用的算法。这个算法主要是针对于人类基因组。

2: bwa   mem   [ - a C H M p P ]   [ - t   n T h r e a d s ]   [ - k   m i n S e e d L e n ]   [ - w   b a n d W i d t h ]   [ - d   z D r o p o f f ]   [ - rs e e d S p l i t R a t i o ]   [ - c   m a x O c c ]   [ - A   m a t c h S c o r e ]   [ - B   m m P e n a l t y ]   [ - O   g a p O p e n P e n ]   [ - Eg a p E x t P e n ]   [ - L   c l i p P e n ]   [ - U   u n p a i r P e n ]   [ - R   R G l i n e ]   [ - v   v e r b o s e L e v e l ]   d b . p r e f i xr e a d s . f q   [ m a t e s . f q ]

           BWA-MEM 算法比对70bp-1Mbp的输入序列。简要的说,算法主要是通过最大精确匹配作为种子比对,然后基于Smith-Waterman算法进行仿射空位罚分。

           如果mate.fq文件是缺失的和选项-p并未设置。这个命令说明数据为单端测序。如果mates.fa存在,命令行假设reads.fq的第i行与mates.fq的第i行形成read对。如果-p被使用,命令行假设reads.fq的2i行和2i+1行形成read对。这类文件被称为interleaved。在这种例子中,mates.fq文件被忽略。在paired-end 模式中,mem命令行会推断从一批reads中推断reads的方向和插入大小的分布。
选项

 -t  INT  线程数目

 -k  INT  最小种子长度。少于INT的匹配将会被忽略。匹配的速度通常对于这个值不敏感,除非明显偏离20.  [19]

  -w  INT  空值宽度。必要的说,gaps长于INT将不会被发现。需要注意最大gap长度同时受到评分矩阵和hit长度所影响。并不只由这个选项确定。[100]

 -d  INT  off-diagonal-X-dropoff (z-dropoff)。如果最好和目前的延伸分数差距大于 |i-j|*A+INT,将停止延伸,其中i和j是被比对和参考基因组中的位置。A是匹配得分。Z-dropoff 类似于BLAST 中的X-dropoff,除了该算法中并没有空格罚分。Z-dropoff不仅避免了不必要的延伸,同时减少了在较差的延伸比对中的比对。 [100]

 -r  FLOAT 引发长度大于minSeedLen *FLOAT的重新搜索。这是启发式算法调节算法性能的关键参数。更大的值产生更少的seeds,导致更快的比对速度但是更低的准确性。[1.5]

   -c  INT 丢弃大于INT出现次数的MEM。这是一个不敏感参数。

 -p  在paired-end 模式中,运行SW搜索得到缺失的命中。

 -A   INT  匹配得分。 [1]

 -B   INT 错配得分。序列的错误率估计方法:
     {0.75*exp[-log(4)*B/A]}.    [4]

 -o   INT  空值罚分。 [-6]

 -E   INT  空值延伸罚分。一个长度为K 的罚分为 O+K*E

 -L   INT          裁剪罚分。
-U INT  对于未配对read对罚分。对于未配对的read对BWA—MEM以scoreRead1+scoreRead2-INT进行评分。评分scoreRead1+scoreRead2-insertPenalty。比较这两种评分从而确定是否应该强制配对。 [9]

-p   假设第一个输入文件为interleaved 配对FASTA\Q文件。

-R STR 完成read group header行 ’\t‘可以在字符串中使用,将会在SAM文件中转换成SAM文件。read group ID也会附在输出文件每一个reads中。  

-T INT 不要输出比对分数低于INT的比对。这个结果只影响最终结果。

-a 输出所有的比对以单端或未配对双端测序方式。

-c 将FASTA/Q的comment 追加到SAM输出中。选项可以将reads meta信息转移到SAM输出。注意FASTA/Q comment必须符合SAM特定要求。不正确的格式将导致不争取的SAM输出。

-H  使用大写H在SAM输出文件中,这个选项可以显著的减少输出文件的复杂度。当比对长或Bac序列时。

-M 标记短split hit为第二个。

-v INT 控制输出的冗长程度。这个选项并未在BWA完全被支持。理想的,值0 表示不输出到stderr。1表示只输出error。2表示warning和errror。3表示所有信息。4表示对于debug的更高信息。当选项是4时候,输出并不是SAM。 [3]

Reference: 1:推荐一个SAM文件中flag含义解释具 2:生物信息学100个基础问题 —— 第17题 高通量测序的回贴问题 II (SAM文件与BAM文件) 3:生信:2:sam格式文件解读 4:http://samtools.github.io/hts-specs/SAMv1.pdf

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

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券