BBQ(生信基础问题18,19,20)-----mapping专

今天我们来关注一下SAM文件的第五列MAPQ和第六列CIGAR(第一次看到这个词反应半天，雪茄是神马鬼哈哈)

1：什么是MAPQ？

根据SAM文件的官方定义：
  MAPQ: Mapping Quality. It equals -10 log10 Pr{mapping position is wrong}, rounded to the nearest integer. A value 255 indicates that the mapping quality is not available.
用公式来表示的话

MAPQ=-10 * log10{mapping出错的概率}

如何理解这个公式呢？

我们肯定是希望mapping的出错率越低越好的。举个例子，如果是在mapping的过程中，有100个碱基里面有1个碱基的错配的话，它的出错率是1%。这样算出来的MAPQ是20，所以，当MAPQ的值越大，它的出错率越低。而且一般最好的值是60。说明它的错误率是10^-6.

我们来看一个实际的例子。截图里MAPQ出现了3个不同的值：60,22,0
那么0就是说没有匹配上。如果是255的话，说明是占位符。

SAM文件截图

我们来看看BBQ的问题：
1.1 如果mapping的时候输入的是FASTA文件，那么MAPQ还有意义吗？为什么？

首先，我们回顾一下.Fa文件的格式： fasta文件格式，图片来自于BBQ17问
文件的第四列的信息是每个碱基的质量分数。而MAPQ是指mapping到reference序列上的匹配情况，二者是不同的概念。
1.2. 不同的比对软件比如bwa与bowtie2，计算出来的MAPQ意义相同吗？为什么？

bwa，bowtie，bowtie2都是基于BWT算法。但是它们有轻微的区别。BWA的准确率高，是SNP分析的首选比对软件。而Bowtie2的运算速度快，被广泛地应用于ChIP-seq, RNA-seq的分析当中。Bowtie2是支持局部比对，全局不对。而且还支持GAP.所以对于MAPQ的值可能会出现不同。
1.3. 请写出samtools view 命令 获得MAPQ 大于等于20的sam文件，假设原始的sam文件名为raw.sam，过滤后的sam文件名为filter_MAPQ20.sam

samtools view raw.sam -q 20  > filter_MAPQ20.sam

------------------------------------------------------分割线-------------------------------------------------

2：我们来看看第六列“雪茄”列---CIGAR

2.1 CIGAR是啥？
CIGAR = Concise Idiosyncratic Gapped Alignment Report（简洁独特的间隙对齐报告）

举个栗子：请大家看看截图的红箭头的部分，是24S126M和140M1I9M

SAM文件截图

M:Match 表示匹配上了
S:soft clipping (clipped sequences present in SEQ) 说明存在可变剪接。
I:Insertion 表示插入
所以对于截图里的第一个箭头就是有24个可变剪接的位置，126个匹配上了
对于第二个箭头来说有前140个匹配上了，有一个插入，后面9个可以匹配上。

---------------------------------------------------------分割线------------------------------------------------
我们来看BBQ的问题：
2.1. M,I,D,N分别是什么意思？如果1条序列的CIGAR=150M， 那么是不是可以说这150bp的区域中没有mismatch（错配）的现象？
M是match（匹配上或者没匹配上），所以150M不等于完全匹配上
I是Insertion(插入)
D是Deletion（删除）
N是可变剪接的位置（内含子等）

官方文档的解释

2.2 如果1条序列来自于成熟的mRNA，在mapping到基因组的时候会有什么问题？如果这条序列中间正好跨过了200bp的intron，前后各有75bp mapping到了exon上，那么这条序列的CIGAR值应该怎么写？
存在外显子和内含子的位置问题，因为成熟的mRNA都是只有外显子没有内含子的。但是基因组是外显子和内含子都有的。
75M200N75M
2.3. 根据下图提示，请理解clip的含义，无论是softclip还是hardclip。

引自 http://bioinformatics.cvr.ac.uk/blog/tag/cigar-string/

S：soft clipping (clipped sequences present in SEQ)
H：hard clipping (clipped sequences NOT present in SEQ)
clipped均表示一条read的序列被分开，之所以被分开，是因为read的一部分序列能匹配到第三列的RNAME序列上，而被分开的那部分不能匹配到RNAME序列上。

Soft&Hardde 例子

----------------------------------------------分割线---------------------------------------------------------

3:让人摸不着头脑的TLEN,SEQ,QUAL

我们来关注一下第7-11列的信息

表1 SAM格式的标准11列信息介绍

第7列，一般情况下是指Pair read的另一半的比对的参考基因组；，一般“＊”代表完全没有比对上，“＝”代表完全比对

第8列， 如果是双端测序，一般情况下是指Pair read的另一半的比对的参考基因组的坐标；

第9列，可以简单理解为这1对read比对到基因组上以后，上游第1个碱基到下游最后1个碱基的距离。如果用负号表示是下游的序列；如果是正数表示为上游的序列；如果是0表示只是单端比对上；

第10列，进行比对read的序列信息；

第11列，进行比对read的质量信息；

我们来看看BBQ的问题:

1. 图中第20行，第9列记录了TLEN值，列出算式计算TLEN值。

   
   sam图
   

   首先，它的TLEN是负值肯定是比对到了下游的序列上面了。那么这个计算就是：
   -（11123-10946+145）=-322

   
   计算过程

2. 如果使用FASTQ文件作为input，第11列的质量值是否还有意义？为什么？
   没有，因为.Fq自带碱基的质量信息

3. 有没有可能通过SAM文件，提取里面的序列信息并转换成FASTQ格式的文件？如果可能，请你写出程序思路。

samtools view -b -h -s  raw.sam > raw.bam
samtools bam2fq raw.bam > raw.fq

samtools view 帮助文档 samtools bam2fq帮助文档

Reference:
1:生物信息学100个基础问题 —— 第18题 高通量比对的质量
2:生物信息学100个基础问题 —— 第19题 SAM/BAM中的CIGAR值
3:生物信息学100个基础问题 —— 第20题 SAM/BAM中的其它重要信息列
3:生信：2：sam格式文件解读