目前的算法为 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
网友评论