生物資訊資料格式：sam,bam格式
文章目錄
- 資料獲取
- 格式說明
- 例項演練
- 統計共有多少條reads（pair-end-reads這裡算一條）參與了比對參考基因組
- 統計共有多少種比對的型別（即第二列數值有多少種）及其分佈
- 篩選出比對失敗的reads，看序列特徵
- 比對失敗的reads區別成單端失敗和雙端失敗情況，並且拿到序列ID
- 篩選出比對質量值大於30的情況（看第五列）
- 篩選出比對成功，但是並不是完全匹配的序列
- 篩選出insert size長度大於1250bp的pair-end reads
- 統計參考基因組上面各條染色體的成功比對reads數量
- 篩選出原始fq序列裡面有N的比對情況
- 篩選出原始fq序列裡面有N，但是比對的時候卻是完全匹配的情況
- sam檔案裡面的標頭檔案行數
- sam檔案裡面每一行的tags個數一樣嗎？
- sam檔案裡每一行的tags個數分別是多少？
- sam檔案裡記錄的參考基因組染色體長度分別是
- 找到比對情況有insertion情況的
- 找到比對情況有deletion情況的
- 取出位於參考基因組某區域的比對記錄，比如5013到50130區域
- 把sam檔案按照染色體及起始座標排序
- 找到102M3D11M的比對情況，計算其reads片段長度
資料獲取
首先安裝bowtie短序列比對軟體
bam和sam檔案主要是bwa、bowtie、tophat等序列比對工具產生的
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
unzip bowtie2-2.3.4.3-linux-x86_64.zip
ln -s ~/local/app/bowtie2-2.3.4.3-linux-x86_64/bowtie2 ~/bin
生成sam,bam檔案
cd ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads/
bowtie2 -x .. bam 是sam二進位制格式
格式說明
sam(Sequence Alignment Mapping) 序列比對對映
sam分為兩部分，註釋資訊（header section）和比對結果部分（alignment section）
註釋資訊可有可無，都是以@開頭，用不同的tag表示不同的資訊，主要有
- @HD，說明符合標準的版本、對比序列的排列順序；
- @SQ，參考序列說明；例如
@SQ SN:chrY LN:57227415 - @RG，比對上的序列（read）說明；
- @PG，使用的程式說明；例如
@PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa mem /home/sunchengquan/data/database/hg38.fa IonXpress_001.fq - @CO，任意的說明資訊。
比對結果部分（alignment section），每一行表示一個片段（segment）的比對資訊，包括11個必須的欄位（mandatory fields）和一個可選的欄位，欄位之間用tag分割。必須的欄位有11個，順序固定，不可用時，根據欄位定義，可以為’0‘或者’*‘
-
QNAME,序列的名字（Read的名字）
-
FLAG, 概括出一個合適的標記，各個數字分別代表
1 read有多種測序資料，一般理解為有雙端測序資料 read paired
2 比對結果是一個pair-end比對，一正一負完美的比對上，read mapped in proper pair
4 這條read沒有比對上，read unmapped
8 這個序列是pair中的一個但是沒有比對上 mate unmapped
16 序列與參考序列反向互補 read reverse strand
32 這個序列在pair-end中的的mate序列與參考序列反響互補 mate reverse strand
64 序列是 mate 1 first in pair
128 序列是 mate 2 second in pair
256 第二次比對 not primary alignment
… -
RNAME,參考序列的名字（染色體）
-
POS,在參考序列上的位置（染色體上的位置）
-
MAPQ, mapping qulity 越高則位點越獨特
bowtie2有時並不能完全確定一個短的序列來自參考序列的哪個位置，特別是對那些比較簡單的序列。但是bowtie2會給出一個值來顯示這個段序列來自某個位點的概率值，這個值就是mapping qulity。Mapping qulity的計算方法是：Q=-10log10p，Q是一個非負值，p是這個序列不來自這個位點的估計值。
假如說一條序列在某個參考序列上找到了兩個位點，但是其中一個位點的Q明顯大於另一個位點的Q值，這條序列來源於前一個位點的可能性就比較大。Q值的差距越大，這獨特性越高。 -
CIGAR,代表比對結果的CIGAR字串，如37M1D2M1I，這段字元的意思是37個匹配，1個參考序列上的刪除，2個匹配，1個參考序列上的插入。M代表的是alignment match(可以是錯配)
standard cigar:
M match
I insertion
D deletionextended cigar
N gap
S substitution
H hard clipping
P padding
= sequence match
X sequence mismatch -
RNEXT, mate 序列所在參考序列的名稱; 下一個片段比對上的參考序列的編號，沒有另外的片段，這裡是’*‘，同一個片段，用’=‘；
-
PNEXT, mate 序列在參考序列上的位置;下一個片段比對上的位置，如果不可用，此處為0；
-
TLEN，估計出的插入片段的長度，當mate 序列位於本序列上游時該值為負值。Template的長度，最左邊得為正，最右邊的為負，中間的不用定義正負，不分割槽段（single-segment)的比對上，或者不可用時，此處為0；
-
SEQ，read的序列;序列片段的序列資訊，如果不儲存此類資訊，此處為’*‘，注意CIGAR中M/I/S/=/X對應數字的和要等於序列長度；
-
QUAL，ASCII碼格式的序列質量;序列的質量資訊，格式同FASTQ一樣。
-
可選的欄位（field）
NM:i 經過編輯的序列
MD:Z 代表序列和參考序列錯配的字串
AS:i 匹配的得分
insert size示意圖：
insertsize = mate位置B1（第8列）+mate長度（B2-B1）- 比對起始位置A1（4列）
FLAG含義查詢網站： http://broadinstitute.github.io/picard/explain-flags.html ， 假如說標記為以上列舉出的數目，就可以直接推斷出匹配的情況。假如說標記不是以上列舉出的數字，比如說 83=（64+16+2+1），就是這幾種情況值和。
例項演練
統計共有多少條reads（pair-end-reads這裡算一條）參與了比對參考基因組
[sunchengquan 09:59:59 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |wc
20000 391929 7049181
mate1 和 mate2 算一條read
$ cat tmp.sam |grep -v '^@' |cut -f1 |uniq |wc -l
10000
統計共有多少種比對的型別（即第二列數值有多少種）及其分佈
[sunchengquan 10:10:04 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |cut -f2 |sort |uniq -c|sort -nrk1,1
4650 83
4650 163
4516 99
4516 147
213 77
213 141
165 69
165 137
153 73
153 133
136 89
136 165
125 153
125 101
24 65
24 129
16 177
16 113
2 81
2 161
篩選出比對失敗的reads，看序列特徵
[sunchengquan 10:19:05 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print}'|wc -l
1005
[sunchengquan 10:20:28 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print}'|cut -f10 |head -3
GCGTNGTACNNGNNGATGNTTACGACTAANACTTNNNNNTNCNGNNT
ANTGGNTNAGCGGNCGNAGNGNTCNGNAGCNANNAANCTGNTATGNTNNNTNNATCNNNCNNNNGTGNNAANNNCCAGNNNNGGCNNNNCATACANNNTNNNCNTNNNTACANNATANNCGNNATTNNACNGTNNCGNGNNCCGNGTNNNNN
NNTNNNNNAGNTNAANANANGTNTTGNGAANCTNNANACNNCTTTTNNNNTATGNTNNTNAGCCTAG
[sunchengquan 10:20:37 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$10}'|head -3
GCGTNGTACNNGNNGATGNTTACGACTAANACTTNNNNNTNCNGNNT
ANTGGNTNAGCGGNCGNAGNGNTCNGNAGCNANNAANCTGNTATGNTNNNTNNATCNNNCNNNNGTGNNAANNNCCAGNNNNGGCNNNNCATACANNNTNNNCNTNNNTACANNATANNCGNNATTNNACNGTNNCGNGNNCCGNGTNNNNN
NNTNNNNNAGNTNAANANANGTNTTGNGAANCTNNANACNNCTTTTNNNNTATGNTNNTNAGCCTAG
比對失敗的reads區別成單端失敗和雙端失敗情況，並且拿到序列ID
[sunchengquan 10:23:58 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort |uniq -c|grep -w 2 |head -2
2 r1019
2 r1105
[sunchengquan 10:24:20 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort |uniq -c|grep -wv 2 |head -2
1 r1007
1 r1010
篩選出比對質量值大於30的情況（看第五列）
[sunchengquan 10:24:42 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($5>30)print$1,$5}'|head -3
r1 42
r1 42
r2 42
篩選出比對成功，但是並不是完全匹配的序列
[sunchengquan 10:26:19 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6!="*")print$6}'|grep '[IDNSHPX]'|head -4
72M2D119M
126M3D61M
70M1D71M
53M5D134M
篩選出insert size長度大於1250bp的pair-end reads
[sunchengquan 10:48:55 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($9>1250||$9<-1250)print}'|less -S
統計參考基因組上面各條染色體的成功比對reads數量
[sunchengquan 22:25:45 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cut -f3 tmp.sam |sort |uniq -c
426 *
19574 gi|9626243|ref|NC_001416.1|
測試資料只有一個染色體
篩選出原始fq序列裡面有N的比對情況
[sunchengquan 11:35:18 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($10 ~/N/) print $0}' tmp.sam |less -S
篩選出原始fq序列裡面有N，但是比對的時候卻是完全匹配的情況
[sunchengquan 11:42:42 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($10 ~/N/) print}' tmp.sam |awk '{if($6 !~ /[IDNSHPX]/) print}' |less -S
[sunchengquan 11:42:52 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($10 ~/N/) print}' tmp.sam |awk '{if($6 !~ /[IDNSHPX*]/) print}' |less -S
sam檔案裡面的標頭檔案行數
[sunchengquan 12:37:28 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep '^@' tmp.sam
@HD VN:1.0 SO:unsorted
@SQ SN:gi|9626243|ref|NC_001416.1| LN:48502
@PG ID:bowtie2 PN:bowtie2 VN:2.3.4.3 CL:"/home/sunchengquan/local/app/bowtie2-2.3.4.3-linux-x86_64/bowtie2-align-s --wrapper basic-0 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq"
sam檔案裡面每一行的tags個數一樣嗎？
[sunchengquan 13:43:16 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam |cut -f 12-1000|head -6
AS:i:-5 XN:i:0 XM:i:3 XO:i:0 XG:i:0 NM:i:3 MD:Z:59G13G21G26 YS:i:-4 YT:Z:CP
AS:i:-4 XN:i:0 XM:i:4 XO:i:0 XG:i:0 NM:i:4 MD:Z:65T3A19T107T0 YS:i:-5 YT:Z:CP
AS:i:-14 XN:i:0 XM:i:8 XO:i:0 XG:i:0 NM:i:8 MD:Z:0A0C0G0A108C23G9T81T46 YS:i:-5 YT:Z:CP
AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:1C182C3 YS:i:-14 YT:Z:CP
AS:i:-22 XN:i:0 XM:i:8 XO:i:0 XG:i:0 NM:i:8 MD:Z:80C4C16A52T23G30A8T76A41 YS:i:-3 YT:Z:CP
AS:i:-3 XN:i:0 XM:i:3 XO:i:0 XG:i:0 NM:i:3 MD:Z:23C1T39G1 YS:i:-22 YT:Z:CP
sam檔案裡每一行的tags個數分別是多少？
[sunchengquan 13:50:54 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam |cut -f 12-1000|awk '{print NF}' |sort |uniq -c
457 1
548 2
579 8
18416 9
[sunchengquan 13:59:09 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam |cut -f 12-1000|awk '{split($0,a, " ");print length(a)}' |sort|uniq -c
457 1
548 2
579 8
18416 9
sam檔案裡記錄的參考基因組染色體長度分別是
[sunchengquan 14:02:12 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep '^@SQ' tmp.sam |cut -f3|cut -d: -f2
48502
找到比對情況有insertion情況的
[sunchengquan 14:05:50 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6 ~/I/)print}' tmp.sam |less -S
至少有3個I
[sunchengquan 14:12:07 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6 ~/(I.*){3,}/)print}' tmp.sam |less -S
找到比對情況有deletion情況的
[sunchengquan 14:12:27 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6 ~/D/)print}' tmp.sam |less -S
取出位於參考基因組某區域的比對記錄，比如5013到50130區域
[sunchengquan 14:13:04 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($4>5013 && $4<50130) print}' tmp.sam |less -S
把sam檔案按照染色體及起始座標排序
[sunchengquan 14:30:02 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam | awk '{print $4":" $0}' |sort -t: -nrk1,1|less -S
找到102M3D11M的比對情況，計算其reads片段長度
[sunchengquan 14:40:31 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6=="102M3D11M")print $1,"read length:"length($10)}' tmp.sam
r284 read length:113