拿到質(zhì)控后的數(shù)據(jù)就可以開始進行數(shù)據(jù)的比對了。（這里復習一下）需要確定好參考基因組儿倒，參考基因組選擇標準：質(zhì)量高晾嘶、完整拐揭、和某一親本相近（是親本之一的基因組當然更好了）
比對軟件
由于轉(zhuǎn)錄組存在可變剪切，所以和DNA重測序數(shù)據(jù)還是有區(qū)別的火脉。簡單來說就是DNA重測序數(shù)據(jù)包含基因組所有的堿基信息（內(nèi)含子外顯子都有）牵舵，但是RNA數(shù)據(jù)只有表達數(shù)據(jù)（外顯子，同時存在可變剪切）倦挂。所以比對時要注意這一點畸颅。DNA數(shù)據(jù)一般使用BWA，RNA除了BWA還可以使用Hisat2等方援。
BWA
下載安裝
wget http://superb-sea2.dl.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2
tar jxf bwa-0.7.12.tar.bz2
cd bwa-0.7.12 && make
#編譯完成后將軟件寫入環(huán)境變量中
echo 'PATH=$PATH:~/softwarepath/bwa-0.7.12/setup' >> ~/.bashrc && source ~/.bashrc
#測試是否安裝好
bwa
Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.12
Contact: Heng Li <lh3@sanger.ac.uk>
Usage: bwa <command> [options]
Command: index index sequences in the FASTA format
mem BWA-MEM algorithm
fastmap identify super-maximal exact matches
pemerge merge overlapping paired ends (EXPERIMENTAL)
aln gapped/ungapped alignment
samse generate alignment (single ended)
sampe generate alignment (paired ended)
bwasw BWA-SW for long queries
fa2pac convert FASTA to PAC format
pac2bwt generate BWT from PAC
pac2bwtgen alternative algorithm for generating BWT
bwtupdate update .bwt to the new format
bwt2sa generate SA from BWT and Occ
Note: To use BWA, you need to first index the genome with `bwa index'.
There are three alignment algorithms in BWA: `mem', `bwasw', and
`aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
first. Please `man ./bwa.1' for the manual.
使用
bwa包含三種算法：
BWA-backtrack: 比對Illumina 序列没炒，read length： ≤100bp。
BWA-SW: 比對 long read犯戏，也支持剪切比對送火。read length：70bp-1Mbp。
BWA-MEM: 比對 long read先匪，也支持剪切比對种吸。read length：70bp-1Mbp。（該算法更新呀非，準確性和速度都更優(yōu)）
一般使用步驟：
#ref.fa即為參考基因組坚俗，genome為建立索引庫的前綴，不加-p genome岸裙，默認建立索引以ref.fa為前綴
bwa index ref.fa -p genome
#因為BSA分析中的重測序數(shù)據(jù)大多是150 PE猖败，這里就直接介紹mem算法的用法，別的算法不再詳細介紹
bwa mem [-aCHMpP] [-t nThreads] [-k minSeedLen] [-w bandWidth] [-d zDropoff] \
[-r seedSplitRatio] [-c maxOcc] [-A matchScore] [-B mmPenalty] [-O gapOpenPen] \
[-E gapExtPen] [-L clipPen] [-U unpairPen] [-R RGline] [-v verboseLevel] \
db.prefix reads.fq [mates.fq]
常用參數(shù)注釋：
-t INT 線程數(shù)降允，默認是1恩闻。
-M 將 shorter split hits 標記為次優(yōu)，以兼容 Picard’s markDuplicates 軟件拟糕。
-p 若無此參數(shù)：輸入文件只有1個判呕，則進行單端比對；若輸入文件有2個送滞，則作為paired reads進行比對侠草。若加入此參數(shù)：則僅以第1個文件作為輸入(輸入的文件若有2個，則忽略之)犁嗅，該文件必須是read1.fq和read2.fa進行reads交叉的數(shù)據(jù)边涕。
-R STR 完整的read group的頭部，可以用 '\t' 作為分隔符，sam文件會識別功蜓。需要指定RG表頭是因為同一樣品可包括多個不同lane园爷、不同文庫的測序結(jié)果，此外不同樣品比對結(jié)果合并時式撼，也需要通過RG進行標記區(qū)分童社。
-T INT 當比對的分值比 INT 小時，不輸出該比對結(jié)果著隆，這個參數(shù)只影響輸出的結(jié)果扰楼，不影響比對的過程。
-a 將所有的比對結(jié)果都輸出美浦，包括 single-end 和 unpaired paired-end的 reads弦赖，但是這些比對的結(jié)果會被標記為次優(yōu)。
SE（單）：
bwa aln [options] ref.fa read.fq > aln_sa.sai
bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam
PE（雙）：
bwa aln [options] ref.fa read1.fq > aln_sa1.sai
bwa aln [options] ref.fa read2.fq > aln_sa2.sai
bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam
#實例（雙端）浦辨，sam文件即為我們的比對結(jié)果文件蹬竖。
##不指定樣品名
bwa mem ref.fa read1.fq read2.fq > out.sam
##指定樣品名，sample即為樣品名
bwa mem -M -R "@RG\tID:sample\tLB:sample\tSM:sample" ref.fa read1.fastq read2.fastq > out.sam 2> out.sam.log
當我們想要提高比對的速率的時候流酬，可以直接使用shell的基本命令split币厕，或者寫python腳本對測序數(shù)據(jù)進行分割再并行比對（每四行為一組），最后再將比對文件合并起來康吵。
注意：1劈榨、因為涉及到后面的變異檢測，所以當基因組存在染色體長度大于500M時晦嵌，一定要對染色體拆分成小于500M的序列再進行比對同辣，在變異檢測完畢后再進行合并（拆分基因組的腳本我會放在流程腳本里，這里不再詳細記錄）惭载；2旱函、一定不能是基因組拆分進行比對，因為數(shù)據(jù)是無序的描滔，且涉及到多重比對的問題０舴痢！含长！
比對結(jié)果格式解讀
@HD VN:1.5 SO:coordinate
@SQ SN:Chr01 LN:91897250
@RG ID:sample LB:sample SM:sample
@PG ID:bwa CL:~/softwarepath/bwa-0.7.12/setup/bwa mem -M -R "@RG\tID:sample\tLB:sample\tSM:sample" ref.fa read1.fastq read2.fastq > out.sam 2> out.sam.log PN:bwa VN:0.7.17-r1188
E150015117L1C019R0062820889 433 Chr01 186 0 72M78H ChrD02_S6 89876638 0 AAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCTA ?A+<80-?C5-@A5BC<7>A<@CC-2A9CCB67@2BCC/;8:@CCC=66@CCA<>=BCCCAB68B@C<?<=C NM:i:0 MD:Z:72 MC:Z:150M AS:i:72 XS:i:72 RG:Z:CY619 SA:Z:ChrD05_S5,59252485,-,51S99M,0,1;
E150015079L1C036R0342398185 161 Chr01 361 2 6S14M1I86M1I14M4D28M ChrD02_S1 14286 0 CCTAAAAAACCCTAAACCCTAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCAAAAAACCCTAAAAAAACCCTAAACCCTAAACCCTAAACCCTAAACCATAAACCCCTAAACCCTAAAAAACCCTAAACCCTAAACCCTAAACCCT B-AB<CBA/65@=5B<.C96C4:B<#AB2=8);18B7/;=??C86>)@AB<8C4<CCA2<C69@>=?2A=A76=6/867/C=4>A--=(2B=A1;@@/9.#@'4C>).@?92<C421AAC(B2@99-C6#=&'<<9@4C2'B.B$B/CC% NM:i:7 MD:Z:95C18^CCCT28 MC:Z:147M3S AS:i:113 XS:i:109 RG:Z:CY619
E150015079L1C003R0260125142 81 Chr01 368 0 38M1D112M ChrD03 10858 0 AAACCGTAAACCCTAAAGCCTAAATCCTAAACCCTAAAGCAAAAAACCCTAAAAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAAAGCTAAACCCTAAACGCTAAACCCTAAACC @2C;15CCCC26<BCAB$?BCBC>/<.@@CB9C.C?CC#CBBA6BC:C6ACCBCC>CC=@CCCA7A>CACAC6CCACC?1?CCC=B9CCCCC.2-CCC<;C4CCCC;5?CCCCC-CCCAC&(CCC;A*C5CAC7C'9BBCCC?7CCC>8B NM:i:8 MD:Z:5C11C6C13^C0C81C0C13C14 MC:Z:3M1D147M AS:i:108 XS:i:106 RG:Z:CY619
#@開頭的列都是表頭信息行券腔，不能少
HD：（排序類型）
頭部區(qū)第一行：VN是格式版本；SO表示比對排序的類型拘泞，有unknown（default）纷纫，unsorted，queryname和coordinate幾種陪腌。
SQ：染色體號\t染色體長度
RG：樣本信息
PG: 比對命令行辱魁，可以直接由此查到比對方式烟瞧，參考基因組，測序數(shù)據(jù)染簇，樣本名等信息
#比對內(nèi)容行每列內(nèi)容注釋：
第一列：QNAME参滴。比對的序列名稱，類似于E150015117L1C019R0062820889锻弓。
第二列：FLAG Bwise FLAG砾赔，比對類型：paring，strand弥咪，mate strand等过蹂。（這個值需要進行計算）
第三列：RENAME 比對上的參考序列名，類似于Chr01
第四列：POS 1-Based的比對起始物理位置
第五列：MAPQ 比對質(zhì)量
第六列：CIGAR Extended CIGAR string（操作符：MIDNPSH=X）比對結(jié)果信息聚至，包含數(shù)字和字符。
第七列：MRNM 相匹配的另外一條序列本橙，比對上的參考序列名扳躬，“=”：比對到一個位置；“*”：無匹配甚亭；染色體號：比對多個位置贷币。
第八列：MPOS 1-Based leftmost Mate Position 與該reads對應的mate pair reads的比對位置
第九列：ISIZE 插入片段長度
第十列：SEQ 和參考序列在同一個鏈上比對的序列（若比對結(jié)果在負義鏈上，則序列是其反向重復序列亏狰，反向互補序列）
第十一列：QUAL 比對序列的質(zhì)量（ASCII-33=Phred base quality）reads堿基質(zhì)量值
可選列：以TAG：TYPE：VALUE的形式提供額外的信息役纹。
比對文件格式處理
samtools
下載安裝
# 建議直接conda安裝
conda create -n samtools
conda activate samtools
conda install -c bioconda samtools -y
#測試是否安裝成功
samtools
Program: samtools (Tools for alignments in the SAM format)
Version: 1.3.1 (using htslib 1.3.1)
Usage: samtools <command> [options]
Commands:
-- Indexing
dict create a sequence dictionary file
faidx index/extract FASTA
index index alignment
-- Editing
calmd recalculate MD/NM tags and '=' bases
fixmate fix mate information
reheader replace BAM header
rmdup remove PCR duplicates
targetcut cut fosmid regions (for fosmid pool only)
addreplacerg adds or replaces RG tags
-- File operations
collate shuffle and group alignments by name
cat concatenate BAMs
merge merge sorted alignments
mpileup multi-way pileup
sort sort alignment file
split splits a file by read group
quickcheck quickly check if SAM/BAM/CRAM file appears intact
fastq converts a BAM to a FASTQ
fasta converts a BAM to a FASTA
-- Statistics
bedcov read depth per BED region
depth compute the depth
flagstat simple stats
idxstats BAM index stats
phase phase heterozygotes
stats generate stats (former bamcheck)
-- Viewing
flags explain BAM flags
tview text alignment viewer
view SAM<->BAM<->CRAM conversion
depad convert padded BAM to unpadded BAM
使用
這里主要介紹如何將sam文件處理成二進制bam文件，建立索引暇唾，排序促脉，去重，以及統(tǒng)計比對信息
#sam文件查看
less -S out.sam
#bam文件查看
samtools view -hS out.sort.rmdup.bam |less -S
##以上二者區(qū)別就是是否為二進制文件
#sam轉(zhuǎn)bam
samtools view -hb out.sam >out.bam
#拆分并行比對后會涉及到bam合并的問題,merge也可以合并
samtools cat 1.bam 2.bam ... -o out.bam
#bam排序,先建索引（生成*.bam.bai文件策州，當染色體長度過長時則會出現(xiàn)索引無法建立的情況）再排序
samtools index out.bam
samtools sort out.bam -o out.sort.bam
#去重瘸味，s:SE,S:PE,默認雙端
samtools rmdup -S out.sort.bam out.sort.rmdup.bam
#統(tǒng)計比對信息
samtools flagstat out.sort.rmdup.bam > flagstat.txt
#flagstat.txt文件內(nèi)容
7524 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
7523 + 0 mapped (99.99% : N/A)
7524 + 0 paired in sequencing
3760 + 0 read1
3764 + 0 read2
7510 + 0 properly paired (99.81% : N/A)
7522 + 0 with itself and mate mapped
1 + 0 singletons (0.01% : N/A)
12 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
第一行：QC pass的reads的數(shù)量為7514，未通過QC的reads數(shù)量為0够挂，意味著一共有7524*2條reads旁仿；
第四行：重復reads的數(shù)量泥彤，QC pass和failed锯厢，這里都是0條
*第五行：比對到參考基因組上的reads數(shù)量，這里是共7523條比對上了參考基因組昆咽，比對率為99.99%办悟；
第六行：paired reads數(shù)據(jù)數(shù)量尘奏，這里是7524對；
第七行：read1比對上的數(shù)量誉尖；
第八行：read2比對上的數(shù)量罪既；
*第九行：雙端都正確匹配到參考序列的reads數(shù)量，這里是7510，比對率為99.81%琢感；
第十行：雙端都匹配到了參考序列上的reads數(shù)量(但是并不一定比對到同一條染色體上)丢间；
*第十一行：雙端中只有一條與參考序列相匹配的數(shù)量；
*第十二行：雙端比對到不同染色體的數(shù)量驹针；
第十三行：雙端比對到不同染色體的且比對質(zhì)量值大于5的數(shù)量烘挫。
## *為主要數(shù)據(jù)關注行
#覆蓋率統(tǒng)計
/work1/Software/samtools-1.4/setup/bin/samtools depth -a out.sort.rmdup.bam > out_depth.txt 2>out_depth.statlog
#out_depth.txt文件內(nèi)容
#染色體號 物理位置 覆蓋堿基數(shù)
Chr01 1 0
Chr01 2 10
Chr01 3 0
Chr01 4 0
## a為有堿基覆蓋的位點總長度
a=`awk '$3!=0' out_depth.txt|wc -l `
## b是基因組總長度
samtools faidx ref.fa
b=`awk '{sum+=$2}END{print sum} ref.fai`
## a/b 即為覆蓋率了
往期回顧
BSA分析（三）——測序數(shù)據(jù)的質(zhì)控
參考資料
https://blog.csdn.net/genome_denovo/article/details/78712972
