真是好久不見啦！在本該撰寫專題內容的整個八月和九月前半部分，小編經歷了離職、旅遊以及西湖到復旦的開學典禮、寢室大掃除搬遷、生病養病、選課搶課等等諸多磨難終於在中秋節來臨之前又完成了一章的內容（「一把辛酸淚」TnT），目測即將進入持續產出和攝入的狀態。
這講內容主要涉及常用常見的一些測序相關知識：測序數據比對結果的解讀和量化分析，其實也是分析NGS數據的基礎步驟了。
由於長時間沒有更新此專題，先回顧一下這個專題已經講了什麼吧（從五月開始到現在了嗯~ o(￣▽￣)o）：
1. 常見的linux操作：指令、文本（正則表達式）、git協作
Bioinformatic Data Skills 學習專題(4) git
Bioinformatic Data Skills 學習專題(3) 文本操作
Bioinformatic Data Skills 學習專題(2)：一些基礎的shell腳本和伺服器指令
2. 常見生物數據，genomic range
Bioinformatic Data Skills 學習專題(6) Genomic Range之三
Bioinformatic Data Skills 學習專題(6) Genomic Range之二
Bioinformatic Data Skills 學習專題(6) Genomic Range之一
Bioinformatic Data Skills 學習專題(5) Bioinformatics Data
SAM/BAM的數據格式組成：header/section/flags/CIGAR/Mapping quality
SAM/BAM文件處理的命令行工具: samtools
用IGV可視化
用Pysam自己做SAM/BAM處理工具（難度較大，日後再講）
這些內容集中回答了如下幾個問題……Q1: 比對後的文件格式：多了什麼？從上一章我們講到的原始測序數據用軟體進比對和回帖（主流比對軟體如bowtie, hisat，以及最近流行的超快轉錄組比對軟體salmon和kallisto），接下來就要解釋比對結果，它是以BAM/SAM文件格式呈現的，主要包括：
1. 頭信息（Header），以@開頭@SQ, reference sequences，通常是染色體名，主要欄位為SN(sequence name )，後面跟著染色體號；LN是則是sequnce length
@RG, 樣本名和組名(read group and sample metadata)， ID必須是唯一的，作者建議the name of the sequencing run and lane，可以用於後期區分batch effect；同時SM號碼可以用於存儲樣本的metadata信息(individual, treatment group, replicate)；PL 這個域可以存儲測序平臺信息，比如ILLUMINA，PACBIO等等。
@PG，存儲比對軟體的信息，ID是必須的；VN為版本號；CL則是command line的縮寫。通常比對軟體會自動加上這些信息
第一行的比對信息不以@開頭
2. 每個reads的具體信息（Alignmet Section）1 QNAME Query template/pair NAME
2 FLAG bitwise FLAG
3 RNAME Reference sequence NAME
4 POS 1-based leftmost POSition/coordinate of clipped sequence 位置信息
5 MAPQ MAPping Quality (Phred-scaled) 比對鹼基質量，詳見上一節
6 CIGAR extended CIGAR string, matching bases, insertions/deletions, clipping
7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME) MPOS 1-based Mate POSistion,下一個信息？
8 TLEN inferred Template LENgth (insert size)， template length for paired-end read
9 SEQ query SEQuence on the same strand as the reference
10 QUAL query QUALity (ASCII-33 gives the Phred base quality)
我們需要的信息主要從下面三個欄位獲得：
CIGAR string: matches/mismatches, insertions, deletions, soft or hard clipped, and so on.根據描述，就是看看reads是否有前面這些特徵
Mapping qualities: Q = -10 log10P(incorrect mapping position). 鹼基質量，於30就可以了。
相關代碼看SAM文件的header信息
samtools view -H celegans.sam
@SQ SN:I LN:15072434
@SQ SN:II LN:15279421
@SQ SN:III LN:13783801
看bitwise flag信息
$ samtools flags 147
0x93 147 PAIRED,PROPER_PAIR,REVERSE,READ2
$ samtools flags 0x93
0x93 147 PAIRED,PROPER_PAIR,REVERSE,READ2
Q2: 了解了數據格式。處理這些比對好的文件？換一種說法就是，怎麼分析SAM文件以得到所需要的特定reads（正確的比對結果）？1. 格式轉換可以參考samtools的詳細說明書(詳見http://www.htslib.org/)sam到bam文件格式轉化： samtools view-b celegans.sam>celegans_copy.bam包含header信息的轉換： samtools view-h celegans.bam>celegans_copy.sam對map好的文件進行排序： samtools sort celegans_unsorted.bam celegans_sorted由於體積小，一般sort好的bam文件被用於後續的分析，這樣我們用 samtools view的速度就會快很多，一般情況下，sort好的文件分析都需要以fai結尾的index文件： samtools index celegans_sorted.bam
2. reads挑選實例1：使用 samtools view找特定區域上的比對情況作者再github項目上分享了1000 Genomes Project Consortium的數據來實戰，具體可以到https://github.com/vsbuffalo/bds-files/tree/master/chapter-11-alignment 下載
samtools index NA12891_CEU_sample.bam
# 看一下 chromosome 1, 215,906,469-215,906,652上的比對reads情況
samtools view NA12891_CEU_sample.bam 1:215906469-215906652 | head -n 3
SRR003212.5855757 147 1 215906433 60 33S43M = 215906402 [...]
SRR003206.18432241 163 1 215906434 60 43M8S = 215906468 [...]
SRR014595.5642583 16 1 215906435 37 8S43M * 0 [...]
# 我們可以把獲得的信息重定向到別的文本
$ samtools view -b NA12891_CEU_sample.bam 1:215906469-215906652 >
USH2A_sample_alns.bam
# 如果你想獲得多個區域的比對情況，就需要用這些區域的BED文件，同時加上參數 -L
$ samtools view -L USH2A_exons.bed NA12891_CEU_sample.bam | head -n 3
SRR003214.11652876 163 1 215796180 60 76M = 215796224 92 [...]
SRR010927.6484300 163 1 215796188 60 51M = 215796213 76 [...]
SRR005667.2049283 163 1 215796190 60 51M = 215796340 201 [...]
3. reads挑選實例2：使用 samtools view-f篩選特定的比對結果比如我們想看看那些reads沒有比對上哪~ samtools flags unmap告訴我們：0x4 4 UNMAP，4就是沒有map上的bitflag，接下來就要用view -f的參數篩選出這些reads:
$ samtools view -f 4 NA12891_CEU_sample.bam | head -n 3
SRR003208.1496374 69 1 215623168 0 35M16S = 215623168 0 [...]
SRR002141.16953736 181 1 215623191 0 40M11S = 215623191 0 [...]
SRR002143.2512308 181 1 215623216 0 * = 215623216 0 [...]
參數 view-F的含義剛好相反，就是map好的reads信息了
samtools view -F 4 NA12891_CEU_sample.bam | head -n 3
SRR005672.8895 99 1 215622850 60 51M = 215623041 227 [...]
SRR010927.10846964 163 1 215622860 60 35M16S = 215622892 83 [...]
SRR005674.4317449 99 1 215622863 37 51M = 215622987 175 [...]
可以注意到這些bitflag有很多別的數字，想看69含義，如下：
$ samtools flags 69
0x45 69 PAIRED,UNMAP,READ1
同時我們也可以發過來用 samtools flags
$ samtools flags READ1,PROPER_PAIR
0x42 66 PROPER_PAIR,READ1
最後上一個-f -F的組合使用，需要paired且mapped的reads（有沒有像邏輯學的感覺）
0x1 1 PAIRED
$ samtools flags unmap,proper_pair
0x6 6 PROPER_PAIR,UNMAP
$ samtools view -F 6 -f 1 NA12891_CEU_sample.bam | head -n 3
SRR003208.1496374 137 1 215623168 0 35M16S = 215623168 [...]
ERR002297.5178166 177 1 215623174 0 36M = 215582813 [...]
SRR002141.16953736 121 1 215623191 0 7S44M = 215623191 [...]`
不能盲目使用，需要檢驗：
$ samtools view -F 6 NA12891_CEU_sample.bam | wc -l # total mapped and paired
233628
$ samtools view -F 7 NA12891_CEU_sample.bam | wc -l # total mapped, paired,
201101 # proper paired
$ samtools view -F 6 -f 1 NA12891_CEU_sample.bam | wc -l # total mapped, paired,
32527 # and not proper paired
$ echo "201101 + 32527" | bc
233628
Q3: 用samtools過濾得到正確的比對結果，然後？直接看fasta文件和bam/sam文件的關係可以用samtools的功能實現：
samtools tview -p 1:215906469-215906652 NA12891_CEU_sample.bam \
human_g1k_v37.fasta
到了這一步，我們已經能夠保證獲得的結果都是通過質量控制的，所以技術上的問題已經解決，接下來就要檢驗數據在生物學意義上是不是靠譜的。比如看看reads的分布是不是合理，一般都需要結合已有的生物學背景知識來進行檢驗，比如ChIPseq的話就需要看看input是不是均勻分布，IP的結果是不是真的蛋白富集在特定位點，甲基化數據就要看一些特異區域比如CpG島是不是高甲基化的，重複序列是不是低甲基化等等。通常我們需要把這些結果可視化再做檢驗，所以數據量化必不可少，大致過程如下：
step1. 量化步驟和腳本step2. IGV可視化的一些關鍵點我再六月初的時候寫了一篇關於所需要基因組的IGV參考格式製作，大家可以參考該文章: 何快速有效地將新拿到的參考基因組在IGV裡可視化
講到這裡，本章的主要內容差不多結束了，首先介紹了比對數據SAM文件格式的信息組成，然後講了處理和分析SAM文件的命令行工具samtools。那可不可以不用工具自己分析呢？大拿們當然早已推出各種包了，在此我僅僅摘錄一些書中部分（這裡摘錄的目的其實是今後打算推出單獨章節用於此類工具介紹），有興趣的可以先結合本書學習。
補充部分可能會下期或者最後補講：定製化的文件分析（Pysam）和自己編程創建SAM文件