2020-01-14 bwa比对：入门

2020-01-14 本文已影响0人王子威PtaYoth

目前的算法为 bwa mem，先前的算法为bwa aln

bwa的工作原理

所有的比对工具均基于相同的原则：
1. 从参考基因组建立一个索引
2. 将FASTA和FASTQ文件中的序列同索引进行比对

索引构建包括对参考基因组进行预处理，以便程序可以有效地对其进行检索。不同程序将建立不同类型的索引，可能会产生多个名称或扩展名怪异的文件。 因此最好将参考基因组放在单独的目录中。
为参考基因组构建索引可能比较耗时，某些工具可支持下载构建好的索引。

如何使用bwa

首先为参考基因组建立索引

mkdir -p refs
# -p指同时创建可能需要的父目录

#下载埃博拉病毒FASTA格式的基因组
efetch -db=nuccore -format=fasta -id=AF086833 > ~/refs/AF086833.fa

#建立索引
bwa index ~/refs/AF086833.fa

建议用一个固定的变量名储存名称，然后在后续代码中声明变量名即可
例如：

# 指定变量ACC的值，这里指定为AF086833
ACC=AF086833

# 参考基因组保存在本地，$+变量名调用变量即可
REF=refs/$ACC.fa

# 下载序列
efetch -db=nuccore -format=fasta -id=$ACC > $REF.fa

# samtools建立索引，用于IGV可视化（这一步不太懂）
# 整合基因组浏览器(IGV)是一种高性能的可视化工具，用来交互式地探索大型综合基因组数据。
#samtools faidx  描述: index/extract FASTA
samtools faidx $REF

# 建立参考基因组
bwa index $REF

bwa查看bwa函数所有可用参数，包括bwa index 、bwa mem等

除了建好索引的参考基因组，还需要比对序列

SRA Accession: PRJNA257197
SRA（Short Read Archive）属于INSDC，保存高通量测序实验数据。
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA257197/为刚果埃博拉病毒的测序数据，1个run代表对一个标本测序，获取runinfo.csv

esearch -db sra -query PRJNA257197  | efetch -format runinfo > runinfo.csv
# esearch -db声明数据库，-query声明检索式
# efetch -format声明数据格式 runinfo为SraRunInfo XML格式

run的检索号：SRR1972739
https://www.ncbi.nlm.nih.gov/sra/?term=SRR1972739#

fastq-dump -X 10000 --split-files SRR1972739
#fastq-dump大概是一个sra文件的解压工具，但我对sra文件格式还不了解。
#--split-files: 将双端测序分为两份,放在不同的文件,但是对于一方有而一方没有的reads直接丢弃
#如果不清楚是单端测序还是双端测序，一律--split-3

同样，善用变量名可以简洁脚本

R1=SRR1972739_1.fastq
R2=SRR1972739_2.fastq

#随时echo查看变量名
echo $R1

对双端测序其中一个read pair跑bwa mem，输出一个SAM（序列比对图）文件，它是当下最新的生物信息学数据格式之一，目前已成为存储和表示所有高通量测序结果的标准方法。

bwa mem $REF $R1 > output.sam

SAM文件内部长这样：

@SQ SN:gi|10141003|gb|AF086833.2| LN:18959
@PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa mem /Users/ialbert/refs/ebola/1976.fa SRR1972739_1.fastq SRR1972739_2.fastq
SRR1972739.1 83 gi|10141003|gb|AF086833.2| 15684 60 69M32S = 15600 -153 TTTAGATTTAACAAGATACCGAGAAAATGAATTGATTTATGACAATAATCCTCTAAAAGGAGGACTCAAATGAGATATTGCAATTGAGTCCTCCTTTTAGA DDDDDEEEEEDEEFFFEDHHHHIJJJJJJJJJJJJJJJJJIGIJJJJJJJJJJJJJJJJJJJIGIGJJJJJJJJJJJJJJIJJJJJJIFHHHHFFFFFCCC NM:i:2 MD:Z:27C16G24 AS:i:59 XS:i:0 SA:Z:gi|10141003|gb|AF086833.2|,15735,+,33M68S,60,0;

SAM文件包含有关样品及其比对的所有已知信息。至此，我们不再查看FastQ文件，因为SAM格式包含FastQ测量中也存在的几乎所有信息。

因为这是一个双端测序数据，因此需要在双端模式下进行比对

bwa mem $REF $R1 $R2 > bwa.sam
#输出的sam文件还包含了双端比对的相关信息

其中列出了一些scoring参数：

获取bwa mem的相关参数

bwa mem获取所有参数
-x ont2d 该参数将scoring设置为Oxford Nanopore MinION测序数据比对的推荐值
-x pacbio 该参数将scoring设置为Pacbio测序数据比对的推荐值

学习资料：
《解读SRA数据库规律一文就够》
https://mp.weixin.qq.com/s/1BTerwyy1vD425bFMPc6RQ
《都8102年了，还用fastq-dump，快换fasterq-dump吧》
https://www.jianshu.com/p/5c97a34cc1ad
《Fastq-dump: 一个神奇的软件》（很有用，对参数提供了指导）
https://www.jianshu.com/p/a8d70b66794c