一睁枕、sam涝缝、bam簡(jiǎn)介

1、sam

• Sequence Alignment/Map format譬重，直譯就是序列比對(duì)格式拒逮；

• 序列比對(duì)是測(cè)序結(jié)果分析中最基礎(chǔ)的一步操作，它主要基于參考序列信息回答了每條reads來(lái)自人的哪條染色體的哪個(gè)區(qū)域臀规；

  
  
  sequence mapping

• 目前的比對(duì)軟件很多滩援，例如bwa，star塔嬉，bowtie2玩徊，算法也各有差異，但究其比對(duì)結(jié)果都是以sam純文本格式展現(xiàn)谨究；

• 在sam格式中恩袱，對(duì)比對(duì)情況、比對(duì)質(zhì)量胶哲、PE雙端比對(duì)等各種比對(duì)可能結(jié)果都有全面的記錄畔塔，具體會(huì)在下面介紹。

2鸯屿、bam

• 簡(jiǎn)單來(lái)說(shuō)澈吨，bam格式就是sam的二進(jìn)制版本；
• 因?yàn)榧陌冢瑂am最大的缺點(diǎn)就是文件體積太大谅辣，而bam格式可完整保留sam信息同時(shí)，可明顯降低文件大猩裟铡桑阶；所以一般僅保留bam文件進(jìn)行后續(xù)分析柏副。
• samtools可實(shí)現(xiàn)sam與bam間的轉(zhuǎn)換，以及查看bam內(nèi)容蚣录，具體操作示例見(jiàn)筆記最后

3搓扯、序列比對(duì)小例子

• 如下圖所示，準(zhǔn)確來(lái)說(shuō)是5條reads的比對(duì)結(jié)果包归，其中具體比對(duì)情況為--

  
  
  image.png
  
  

  （1）r001/1與r001/2是PE測(cè)序的pair/mate reads，現(xiàn)在遇到的基本都是雙端測(cè)序铅歼；
  （2）r002公壤、r003、r004都是單獨(dú)的一條read(在這個(gè)圖片結(jié)果是這樣椎椰，也許對(duì)應(yīng)的pair reads比對(duì)到很遠(yuǎn)的位置厦幅，甚至其它染色體上)；
  （3）r003出現(xiàn)兩次是因?yàn)槠淇杀葘?duì)到參考序列的不同位置慨飘，但是注意 -r003表示是其reverse reads比對(duì)上的确憨；
  （4）reads序列中的大寫字母表示與參考序列對(duì)應(yīng)堿基相同，小寫字母則表示不一致瓤的；星號(hào)或點(diǎn)號(hào)表示考慮insert/delete情況休弃，reads才能比對(duì)到參考序列中；

• 而sam格式有系統(tǒng)的規(guī)則圈膏，可簡(jiǎn)潔塔猾、全面的描述上述比對(duì)結(jié)果，如下圖所示稽坤。

  
  
  image.png

在PE雙端測(cè)序中 pair reads與mate pairs概念基本相同丈甸，區(qū)別見(jiàn)https://www.biostars.org/p/77293/

2、sam格式剖析

image.png

2.1 header部分（optimal）

• 以@開頭為標(biāo)志尿褪，從不同角度為比對(duì)結(jié)果補(bǔ)充信息睦擂；

• 常見(jiàn)的有以下四類

  
  
  sam header
• （1）@HD 如果有，必須放在第一行
  VN表示sam格式的版本號(hào)杖玲，SO說(shuō)明比對(duì)結(jié)果是否排序了
• （2）@SQ 說(shuō)明參考基因組的情況顿仇，一般會(huì)有很多行
  SN表示染色體名字，對(duì)應(yīng)比對(duì)結(jié)果的RNAME列摆马；LN表示該染色體的長(zhǎng)度
• （3）@PG則說(shuō)明比對(duì)軟件情況
  ID表示程序ID號(hào)(一般與PN相同)夺欲，PN表示軟件名，VN表示軟件版本號(hào)今膊，CL表示產(chǎn)生該sam的命令語(yǔ)句
• （4）@RG主要說(shuō)明測(cè)序樣本信息
  ID些阅，在一條lane只有一個(gè)sample情況下，可用flowcell+lane序號(hào)表示斑唬；SM表示sample名市埋；LB表示原始測(cè)序文庫(kù)名黎泣；PL表示測(cè)序平臺(tái)

2.2 alignment section（required）

• alignment section是sam的主體部分，其中那個(gè)每行代表一次比對(duì)結(jié)果缤谎，以及對(duì)應(yīng)詳細(xì)的比對(duì)信息抒倚。

• 如下圖，每條alignment record可分為9大部分坷澡，11+n列(metadata不固定)

  
  
  alignment record

第一列：read name

在雙端測(cè)序托呕，一般至少有兩條相同記錄的read name（pair reads），但有時(shí)會(huì)有很多條(>2)記錄是因?yàn)?/p>

• 首先是Chimeric alignment情況频敛，就是類似最一開始介紹的r003的比對(duì)情況项郊，由于序列的特性，可能比對(duì)到不同的比較合適區(qū)域斟赚；
• 然后是Multiple mapping情況着降，這主要是由于reads序列與基因組的重復(fù)序列特征導(dǎo)致的（在人基因組中約占20%）。
  
  read name

第二列：flag

• 簡(jiǎn)單來(lái)說(shuō)就是設(shè)置該條reads的比對(duì)屬性拗军；
• 表示某條reads是否有12條特征中的1至多條(多選題)任洞，結(jié)果用10進(jìn)制數(shù)值加和表示
  
  image.png
  
  簡(jiǎn)單理解所有選項(xiàng)含義---
  1表示是雙端測(cè)序；
  2表示pair reads均比對(duì)到了參考序列(good!)
  4表示這條read沒(méi)比對(duì)到參考序列
  8表示這條read的mate read沒(méi)比對(duì)到參考序列
  16表示這條read是reverse(反向)比對(duì)到參考序列的
  32表示這條read的mate read是reverse(反向)比對(duì)到參考序列的'
  64表示這條read是pair read的第一條read(左邊的)
  128表示這條read是pair read的第二條read(右邊的)
  256表示是這條read的Multiple mapping reads比對(duì)情況

其余三種情況比較少見(jiàn)发侵，不作記錄了交掏。關(guān)于secondary alignment與supplementary alignment的區(qū)別我認(rèn)為就是分別對(duì)應(yīng)Multiple mapping reads與Chimeric alignment比對(duì)情況。詳見(jiàn)https://sourceforge.net/p/samtools/mailman/message/33235303/

• 某一加和結(jié)果數(shù)值只會(huì)是一組屬性特征值的和刃鳄，不會(huì)出現(xiàn)不同組合均能加和等于一個(gè)特定和的結(jié)果耀销；
• 例如下圖，具體flag值含義參考上面解釋
  113=1+16+32+64
369=1+16+32+64+256
117=1+16+32+128
417=1+32+128+256
  
  image.png

這個(gè)網(wǎng)頁(yè)https://broadinstitute.github.io/picard/explain-flags.html提供了根據(jù)輸入的結(jié)果flag來(lái)分解原特征組合的功能铲汪，挺方便的熊尉。

alignment record

第三列&第四列：read的位置

• RNAME: Reference sequence NAME of the alignment；與header部分的@SQ SN保持一致掌腰；
• POS: read比對(duì)到參考序列(染色體)的最左邊位置(基于染色體狰住，且為1-based)；
• 如果只有pairs reads的一個(gè)成功比對(duì)到參考序列齿梁，那么未比對(duì)的reads的
  RNAME與POS應(yīng)與比對(duì)上的reads的描述一致催植；
• 如果pair reads均未成功比對(duì)到參考序列上，則RNAME列均為*勺择，POS列均為0
  
  unmapped reads

簡(jiǎn)單來(lái)說(shuō)1-based就是序列第一個(gè)堿基序號(hào)為1创南，例如SAM, VCF, GFF and Wiggle formats等；0-based就是序列第一個(gè)堿基序號(hào)為0省核，例如 BAM, BCFv2, BED, and PSL formats等稿辙。

第五列：MAPQ

• MAPping Quality. It equals ?10log 10 Pr{mapping position is wrong}, rounded to the nearest
  integer.
• 一般來(lái)說(shuō)MAPQ值越大，表示MAPping Quality越高气忠；但255表示不可計(jì)算該read的MAPQ(我理解就是指unmapped的reads)

第六列：CIGAR

• 這一列主要來(lái)描述序列的具體比對(duì)情況邻储，以堿基為單位赋咽，主要有如下圖幾種比對(duì)類型
  
  image.png
  
  M表示read堿基比對(duì)到了參考序列的堿基上(有趣的是允許夾雜mismatch的情況，為SNP variant calling基礎(chǔ))
  I相比于參考序列吨娜，read里插入insert了新堿基
  D相比于參考序列脓匿，read里刪除了部分堿基
  N比較特殊一般用于mRNA-seq結(jié)果里，表示比對(duì)到了intron內(nèi)含子區(qū)域宦赠；
  S/H表示read序列兩端基本沒(méi)能比對(duì)到參考序列陪毡，而中間部分可以比對(duì)的情況；具體區(qū)別可參考http://www.aigenetic.com/index.php/2018/03/19/soft-clip-%E4%B8%8E-hard-clip%E7%9A%84%E5%90%AB%E4%B9%89%E6%8F%8F%E8%BF%B0/
  P針對(duì)是多條reads比對(duì)到同一參考序列區(qū)域勾扭，其中有一條reads存在insert的情況毡琉，具體見(jiàn)https://davetang.org/wiki/tiki-index.php?page=SAM#Padded_alignment
• Sum of lengths of the M/I/S/=/X operations shall equal the length of SEQ，其中=X不常見(jiàn)尺借。并且第四列read的起始位置也是從M/I/S/=/X開始計(jì)算的

其實(shí)總結(jié)來(lái)看，N與S/H分別表示了兩種特殊比對(duì)情況spliced alignment與clipped alignment精拟，前者是中間沒(méi)比對(duì)上燎斩，而兩端比對(duì)上了；后者使中間比對(duì)上蜂绎，而兩端沒(méi)比對(duì)上栅表。

clipped alignment -- 3S8M1D6M4S

spliced alignment -- 9M32N8M

• 舉個(gè)例子，如下圖表示3 bases aligned followed by 1 base deleted, 2 next ones aligned, 1 base inserted and the last one aligned.

  
  
  image.png

alignment record

第七师枣、八怪瓶、九列：描述mate read信息

• RNEXT是該read的mate read比對(duì)到的參考染色體，有三種情況践美；
  （1）=：read與mate read比對(duì)到同一染色體上洗贰；
  （2）若mate read比對(duì)到其它染色體上，則用相應(yīng)的染色體名稱即可（SN）
  （3）若mate read未比對(duì)成功陨倡，用*表示
• PNEXT是該read的mate read比對(duì)到的參考染色體的起始位置敛滋；
• TLEN列主要描述pair reads的"相隔距離"，如下圖兴革。
  （1）僅考慮M/I/=/X(excludes soft-clipped bases)情況；
  （2）pair reads的該列絕對(duì)值相同，只是左邊的reads為正值艺晴，右邊的reads為負(fù)值兽赁；
  （3）如果pair reads分別比對(duì)到不同染色體上，那么該值就是0
  
  image.png

特殊的TLEN：Note: these two definitions agree in most alignments, but differ in the case of overlaps where the first segment aligns beyond the start of the last segment.

image.png

第十擎勘、十一列：序列信息

• 第十列：原始序列組成咱揍，等于原來(lái)fastq的第二行；
• 第十一列：序列的Phred質(zhì)量值棚饵，等于原來(lái)fastq的第四行述召；詳見(jiàn)之前fastq的筆記朱转。

第12+n列：metadata（optimal）

• TAG:TYPE:VALUE格式：TAG表示標(biāo)簽名，一般是兩個(gè)大寫字母积暖；TYPE表示VALUE的數(shù)據(jù)類型藤为；VALUE表示該read的VALUE值
• Predefined standard tags可參考：https://samtools.github.io/hts-specs/SAMtags.pdf
• 以BWA mem比對(duì)結(jié)果產(chǎn)生的TAG結(jié)果為例，解釋如下圖夺刑。
  
  image.png

至此缅疟，關(guān)于sam格式的基本介紹大致如上，主要參考了http://samtools.github.io/hts-specs/SAMv1.pdf教程手冊(cè)遍愿，其中還有很多深入的知識(shí)存淫，值得以后深入探索~~

image.png

三、samtools轉(zhuǎn)換sam沼填、bam

• sam轉(zhuǎn)bam

samtools view -bS SRR1663608.sam > SRR1663608.bam

• bam轉(zhuǎn)sam

samtools view -h -o SRR1663608.sam SRR1663608.bam

• bam結(jié)果統(tǒng)計(jì)

#查看bam
samtools view -h SRR1663608.bam | more
#the number of records (alignments) 
samtools view -c SRR1663608.bam
#Displays basic alignment stats based on flag
samtools flagstat SRR1663608.bam

更多關(guān)于sam/bam格式的操作桅咆，可見(jiàn)之間生信技能樹Jimmy大神布置的一些練習(xí)題，我自己也做了坞笙，詳見(jiàn)Linux生信練習(xí)3--sam/bam

筆記圖片大多來(lái)自網(wǎng)上岩饼，侵刪~ 筆記中如有錯(cuò)誤之處，歡迎指出薛夜！