map與alignment

新一代測序資料的處理過程中，map/alignment總是同時出現，那麼二者的區別是什麼呢？大致上可以這樣理解：map這個動作表示的是：把短片段對應到參考基因組的某個位置上。這個動作只回答“短片段來自長片段的什麼位置？”這個問題；而align這個動作早在二代測序之前很久就有了，這個動作表示兩條（或多條）序列（不管是核酸序列還是蛋白序列）的比較，回答的問題是“這兩條序列之間有什麼異同？”。顯然在處理重測序資料的過程中，要先map，再align，即先把read定位到參考基因組的某個位置上，再比較read和參考基因組序列的異同。另外還可以認為：凡是要align，那一定是要先map的；但是map過後，卻不一定要align（比如我只關注CNV或表達量等資訊，而不關心具體序列的時候，就只需要map確定下read來自參考序列的哪個位置，而不需要align獲取具體的序列間異同資訊）。可見在SAM檔案中，map和align這兩個動作都要有，這可能就是這兩個動作總是同時出現的原因。理解了這兩個名字的含義和區別，就可以理解為什麼是“mapping” position，因為確定position資訊只需要mapping。map和align兩個詞，似乎都可以用中文“比對”一詞來表示。下文中如果沒有必要區分二者，將全都譯作“比對”。

bwa

命令演算法：

bwa <command> [options]

bwa index 建立索引
bwa mem BWA_MEM algorithm
bwa bwasw BWA_SW for long queries
bwa aln gapped/ungapped alignment
bwa samse generate alignment(single ended)
bwa sampe generate alignment(paired ended)

可通過bwa index/bwa mem/bwa aln檢視相應命令的幫助

比對流程

-建立一個已知參考基因組的索引

索引的建立是對參考基因組進行重新格式化，每個比對軟體的索引不同。基因組越大，構建索引的時間越長。很多比對工具在其官方網站上有提供已經預建的索引檔案，可直接下載使用。

-FASTA/FASTQ檔案中的reads比對到索引上

示例

參照原文用ebola病毒測序資料和1976年Mayinga參考基因組。
mkdir -p ～/refs/ebola #建立資料夾用於存放參考基因組檔案及將建立的索引存放於此
下載ebola病毒參考基因組檔案，命名為1976.fa
efetch -db=nuccore -format=fasta -id=AF086833 > ~/refs/ebola/1976.fa
建立索引
bwa index ~/refs/ebola/1976.fa
還好前段時間看了下shell程式設計，才看懂文中的變數的用法，基礎差的學生就得想辦法惡補。

REF=~/refs/ebola/1976.fa #將後面的路徑賦值給變數REF
bwa index $REF #shell中用$加變數名用來引用變數值
ls ～/refs/ebola #檢視建立索引後目錄內檔案的變化

1976.fa 1976.fa.amb 1976.fa.ann 1976.fa.bwt 1976.fa.pac 1976.fa.sa

獲取SRR1972739的ebola病毒的測序資料

esearch -db sra -query PRJNA257197  | efetch -format runinfo > runinfo.csv
fastq-dump -X 10000 --split-files SRR2571197#通過fastq-dump命令將sra資料還原為fastq檔案，雙端測序檔案加上引數--split-files

命令結束後會在當前目錄下產生兩個測序資料
SRR1972739_1.fastq和SRR1972739_2.fastq

R1=SRR1972739_1.fastq #定義定量R1
R2=SRR1972739_2.fastq #定義變數R2

執行比對
bwa mem $REF $R1 $R2 > bwa.sam

bwa mem

bowtie

Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. It is particularly good at aligning reads of about 50 up to 100s or 1,000s of characters, and particularly good at aligning to relatively long (e.g. mammalian) genomes. Bowtie 2 indexes the genome with an FM Index to keep its memory footprint small: for the human genome, its memory footprint is typically around 3.2 GB. Bowtie 2 supports gapped, local, and paired-end alignment modes.

演算法可參考天天向上的部落格

工作流程

-建立索引
-將FASTA/FASTQ格式的reads與索引進行比對

bowtie2 -h

Usage:
bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

<bt2-idx> Index filename prefix (minus trailing .X.bt2).
NOTE: Bowtie 1 and Bowtie 2 indexes are not compatible.
<m1> Files with #1 mates, paired with files in <m2>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<m2> Files with #2 mates, paired with files in <m1>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<r> Files with unpaired reads.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<sam> File for SAM output (default: stdout)

<m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
specified many times. E.g. '-U file1.fq,file2.fq -U file3.fq'.

建立索引檔案
bowtie2-build -h

Usage: bowtie2-build [options]* <reference_in> <bt2_index_base>
reference_in comma-separated list of files with ref sequences
bt2_index_base write bt2 data to files with this dir/basename
*** Bowtie 2 indexes work only with v2 (not v1). Likewise for v1 indexes. ***
Options:
-f reference files are Fasta (default)
-c reference sequences given on cmd line (as
<reference_in>)
--large-index force generated index to be 'large', even if ref
has fewer than 4 billion nucleotides
-a/--noauto disable automatic -p/--bmax/--dcv memory-fitting
-p/--packed use packed strings internally; slower, less memory
--bmax <int> max bucket sz for blockwise suffix-array builder
--bmaxdivn <int> max bucket sz as divisor of ref len (default: 4)
--dcv <int> diff-cover period for blockwise (default: 1024)
--nodc disable diff-cover (algorithm becomes quadratic)
-r/--noref don't build .3/.4 index files
-3/--justref just build .3/.4 index files
-o/--offrate <int> SA is sampled every 2^<int> BWT chars (default: 5)
-t/--ftabchars <int> # of chars consumed in initial lookup (default: 10)
--threads <int> # of threads
--seed <int> seed for random number generator
-q/--quiet verbose output (for debugging)
-h/--help print detailed description of tool and its options
--usage print this usage message
--version print version information and quit

bowtie2-build $REF $REF #建立索引

第一個路徑表示參考檔案，第二個表示索引的字首名稱，可以是anything。

ls $REF.*
執行比對
bowtie2 -x $REF -1 $R1 -2 $R2 > bowtie.sam

比對軟體的比較與選擇

SAM檔案

SAM代表Sequence Alignment/Map格式，是一種製表符分隔的文字格式，包含一個可選的頭部分（header section，有人稱之為“註釋部分”），和一個比對部分（alignment section）。如果包含頭部分，那麼頭部分必須置於比對部分之前。頭部分的行以@符號開頭，而比對部分的行不以@符號開頭。比對部分的每一行包含11個必選的欄位，用於說明重要的比對資訊，如比對位置（mapping position）等；另有可變數量的可選欄位，用於儲存其他資訊（flexible）或比對軟體特異的資訊。

header部分

@開頭，後面緊接著一個或兩個字母的頭記錄型別編碼。每行都是以製表符分隔的，除了@CO行，每個資料欄位都服從‘TAG:VALUE’這一格式，其中TAG是一個兩字元的字串，定義了VALUE的格式和內容。所以頭部分的行都符合這一正規表示式： /^@[A-Z][A-Z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/或者 /^@CO\t.*/

比對部分

每行表示一個片段（segment）的線性比對（linear alignment）。每行有11個必選欄位，資訊缺失欄位用0或者“*”表示。

FLAG欄位的理解

SAM 檔案是比對結果的展示, 格式詳細介紹見官方文件, 其中第二列為按位標記的 flag, 這是一種常用且高效的儲存多個布林特徵值的方法, 可通過按位取 '與' 和取 '或' 快速提取對應的值

舉個簡單的例子: 在 SAM 格式中, 當 flag 為 1, 也即對應的二進位制為 01 時, 表示該 read 有多個測序資料 (template having multiple segments in sequencing), 一般理解為有雙端測序資料 (另一條沒被過濾掉), 而 flag 為 2, 也即二進位制 10 時, 表示這條 read 的多個片斷都有比對結果 (each segment properly aligned according to the aligner), 通常理解為雙端 reads 都比對上了, 那麼就可以推斷出 flag 為 3 時, 也即二進位制的 11, 表示該 read 有另一端的 read 並且比對成功, 可以看到, 其實就是 01 加 10, 反之, 當我們看到二進位制的 11, 就可以想到這表示這兩位的狀態均為 1 (通常即為真), 再查詢定義就知道了具體的含義

這樣就通過一個位來反映一個 True/ False 的狀態值, 同時方便了人腦和電腦, 在 SAM flag 中, 目前共有 12 種不同狀態, 列舉如下:
（來源：http://not.farbox.com/post/sam_flag）

SAM Format提供了工具可快速查詢各個value的含義。
將其按需求自由組合, 再結合使用 samtools 軟體 view 命令的 -f (滿足) -F (不滿足) 和 -c (計數) 引數, 就可以實現各種過濾和統計了, 比如 samtools view -f 4 就會提取出未比對上的 reads, samtools view -c -F 4 就能統計出比對上的 reads 數, 也即比對率的分子部分

CIGAR欄位

第10及11列就是我們FASTQ檔案中的序列及phred+33對應的ASCII字元