2019年5月30日，晚上，心情变好，好几天没更新了，看到男朋友在学一款软件，我也近朱者赤，来继续注释Day2-2中NGS分析流程中的一个重要软件——BWA

NGS基础

NGS分析注释

BWA

对应于NGS分析流程的这两步：

BWA是目前常用的将测序回来的reads比对到参考基因组上的软件，简单来说，参考基因组相当于一整块已知的地图，而reads是被切碎的，与已知地图存在少量区别的地图碎片，而BWA要做的就是找到两者之间的相似信息，并将地图碎片拼到已知地图上去。
关于BWA的使用贴来“庐州月光”（他可能喜欢听许嵩的歌）的这篇简书https://www.jianshu.com/p/1552cc6ac3be，明天再参考多几篇博文和实验室常用脚本，进行详细说明。

1. 首先BWA包含了以下3种算法
BWA-backtrack
BWA-SW
BWA-MEM
BWA-backtrack适合比对长度不超过100bp的序列；BWA-SW和BWA-MEM适合于长度为70-1M bp的序列；其中BWA-MEM是最新开发的算法，对于高质量的测序数据，其比对的速度更快，精确度更高，对于70-100bp的reads, BWA-MEM算法在比对长度为70-100bp的序列时，效果比BWA-backtrack 算法的效果更好。总而言之，通常情况下，选择BWA-MEM算法就好。

2. bwa 的源代码存储在github上，链接如下
https://github.com/lh3/bwa
安装过程如下：

git clone https://github.com/lh3/bwa.git
cd bwa
make

安装好之后，会出现一个名为bwa的可执行文件，输入下面命令可以查看帮助信息

./bwa

Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
Contact: Heng Li <<a href="mailto:lh3@sanger.ac.uk" _href="mailto:lh3@sanger.ac.uk">lh3@sanger.ac.uk</a>>
Usage:   bwa <command> [options]
Command: index         index sequences in the FASTA formatmem           BWA-MEM algorithmfastmap       identify super-maximal exact matchespemerge       merge overlapping paired ends (EXPERIMENTAL)aln           gapped/ungapped alignmentsamse         generate alignment (single ended)sampe         generate alignment (paired ended)bwasw         BWA-SW for long queriesshm           manage indices in shared memoryfa2pac        convert FASTA to PAC formatpac2bwt       generate BWT from PACpac2bwtgen    alternative algorithm for generating BWTbwtupdate     update .bwt to the new formatbwt2sa        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.

3. 比对之前首先要对参考基因组建立索引，命令如下：

bwa index ref.fasta

下面是对小鼠基因组构建索引的示例

bwa index mm10.fasta
├── mm10.fasta
├── mm10.fasta.amb
├── mm10.fasta.ann
├── mm10.fasta.bwt
├── mm10.fasta.pac
└── mm10.fastq.sa

/mnt/**/bwa index -a bwtsw ref.fasta

前面是bwa的绝对路径，index -a 可以指示构建索引的算法，包括两种，is和bwtsw，这里使用的是后者。对于大于2G的参考基因组文件需使用bwtsw算法，且使用该算法时，参考基因组大小必须大于10M。
除此之外，还可加入另外的参数：-p，可用于指示输出文件前缀，

即：建立索引：bwa index [-p prefix] [-a algoType] <ref.fa>

在网上查阅其他文章时，中间还有一步：https://www.plob.org/article/7009.html 这篇文章中的步骤与我们常用步骤略有不同，借鉴一部分即可。
寻找输入reads文件的SA坐标

## pair end
bwa aln hg19.fa read1.fq.gz -l 30 -k 2 -t 4 -I > read1.fq.gz.sai
bwa aln hg19.fa read2.fq.gz -l 30 -k 2 -t 4 -I > read2.fq.gz.sai
##single end：
bwa aln hg19.fa read.fq.gz -l 30 -k 2 -t 4 -I > read.fq.gz.sai

主要参数说明：
-o int：允许出现的最大gap数。
-e int：每个gap允许的最大长度。
-d int：不允许在3’端出现大于多少bp的deletion。
-i int：不允许在reads两端出现大于多少bp的indel。
-l int：Read前多少个碱基作为seed，如果设置的seed大于read长度，将无法继续，最好设置在25-35，与-k 2 配合使用。
-k int：在seed中的最大编辑距离，使用默认2，与-l配合使用。
-t int：要使用的线程数。
-R int：此参数只应用于pair end中，当没有出现大于此值的最佳比对结果时，将会降低标准再次进行比对。增加这个值可以提高配对比对的准确率，但是同时会消耗更长的时间，默认是32。
-I int：表示输入的文件格式为Illumina 1.3+数据格式。
-B int：设置标记序列。从5’端开始多少个碱基作为标记序列，当-B为正值时，在比对之前会将每个read的标记序列剪切，并将此标记序列表示在BC SAM 标签里，对于pair end数据，两端的标记序列会被连接。
-b ：指定输入格式为bam格式。
bwa aln hg19.fa read.bam > read.fq.gz.sai

4. 建立好参考基因组之后，就可以进行比对了。不同的比对算法，命令不同。

BWA-backtrack 算法
对应的子命令为aln/samse/sample
单端数据用法如下：

bwa aln ref.fa reads.fq > aln_sa.sai
bwa samse ref.fa aln_sa.sai reads.fq > aln-se.sam

双端数据用法如下：

bwa aln ref.fa read1.fq > aln1_sa.sai
bwa aln ref.fa read2.fq > aln2_sa.sai
bwa sampe ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam

BWA-SW 算法
对应的子命令为bwasw, 基本用法如下

bwa bwasw ref.fa reads.fq > aln-se.sam
bwa bwasw ref.fa read1.fq read2.fq > aln-pe.sam

BWA-MEM
算法对应的子命令为mem, 基本用法如下:

bwa mem ref.fa reads.fq > aln-se.sam
bwa mem ref.fa read1.fq read2.fq > aln-pe.sam

##Align reads with BWA-MEM algorithm/mnt/**/bwa mem -t 6 -M -R "@RG\tID:${lane_id}\tLB:${sample}\tPL:Illumina\tPU:${sample}\tSM:${sample}" \ref.fasta ${read1} ${read2} \> /**/00.BAM/${sample}.bwa.sam \2> /**/00.BAM/${sample}.bwa.log

采用常用的mem算法进行比对，-t 用于指示线程数；-M将较短的split hits标记为secondary，与picard兼容；-R设置reads标头，“\t”分割。-R STR Complete read group header line. ’\t’ can be used in STR and will be converted to a TAB in the output SAM. The read group ID will be attached to every read in the output. An example is ’@RG\tID:foo\tSM:bar’.
ID str：输入reads集ID号；LB：read集文库名；PL：测序平台（illunima或solid）；PU：测序平台下级单位名称（run的名称）；SM：样本名称。

${lane_id}和${sample}是上一步读取.fq文件时生成的变量，在下一篇博客介绍。
其基本语法如下：

bwa mem -R '@RG\tID:foo\tSM:bar\tLB:library1' <ref.fa> <read1.fa> <read2.fa> > aln-pe.sam

对于超长读长的reads,比如PacBio和Nanopore测序仪产生的reads, 用法如下:

bwa mem -x pacbio ref.fa reads.fq > aln.sam
bwa mem -x ont2d ref.fa reads.fq > aln.sam

上述的代码中, ref.fa指的是参考基因组索引的名字，对于上面提供的小鼠的例子而言，参考基因组索引的名字为mm10.fasta, 注意是不包含后缀的。默认用法都非常简单，有时还需要对参数进行调整，以mem子命令为例，常用的参数包括以下几个: -t指定线程数，默认为1，增加线程数，会减少运行时间；-p 忽略第二个输入序列，默认情况下，输入一个序列文件认为是单端测序，输入两个序列文件则是双端测序，加上这个参数后，会忽略第二个输入序列文件，把第一个文件当做单端测序的数据进行比对; -Y参数对数据进行soft clipping, 当错配或者gap数过多比对不上时，会对序列进行切除，这里的切除并只是在比对时去掉这部分序列，最终输出结果中序列还是存在的，所以称为soft clipping。

https://www.plob.org/article/7009.html可能还有的其他步骤：
（1）对sam文件进行进行重新排序（reorder）
由BWA生成的sam文件时按字典式排序法进行的排序（lexicographically）进行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在进行callsnp的时候是按照染色体组型（karyotypic）进行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要对原始sam文件进行reorder。可以使用picard-tools中的ReorderSam完成。

java -jar picard-tools-1.96/ReorderSam.jar
I=hg19.sam
O=hg19.reorder_00.sam
REFERENCE=hg19.fa

注意：
1 这一步的头文件可以人工加上，同时要确保头文件中有的序号在下面序列中也有对应的。虽然在GATK网站上的说明chrM可以在最前也可以在最后，但是当把chrM放在最后时可能会出错。
2 在进行排序之前，要先构建参考序列的索引。

samtools faidx hg19.fa

最后生成的索引文件：hg19.fa.fai。
3 如果在上一步想把大文件切分成小文件的时候，头文件可以自己手工加上，之后运行这一步就好了。

(2) 将sam文件转换成bam文件（bam是二进制文件，运算速度快）
这一步可使用samtools view完成。

samtools view -bS hg19.reorder_00.sam -o hg19.sam_01.bam

5. 排序并生成 .bam
这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。可以使用picard-tools中SortSam完成。
lb

 ## sort bam filejava -Djava.io.tmpdir=/**/tmp -jar /mnt/*/picard-tools-1.124/picard.jar SortSam \INPUT=/mnt/*/00.BAM/${sample}.bwa.sam \OUTPUT=/mnt/*/00.BAM/${sample}.bwa.sort.bam \SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT \>> /mnt/*/00.BAM/${sample}.bwa.log 2>&1 && \rm -v /mnt/*/00.BAM/${sample}.bwa.samecho "[`date`]: Start marking duplicates ${sample} ... "

-Djava.io.tmpdir: 获取操作系统缓存的临时目录
基本语法：

java -jar picard-tools-1.96/SortSam.jar
INPUT=hg19.sam_01.bam  这里也可以先转化为.bam文件，再sort
OUTPUT=hg19.sam.sort_02.bam
SORT_ORDER=coordinate

如何理解这里的：

  >> /mnt/*/00.BAM/${sample}.bwa.log 2>&1 && \rm -v /mnt/*/00.BAM/${sample}.bwa.sam

log 2>&1 && rm -v .sam

关于2>&1的含义：把标准错误输出指向log
本来1----->屏幕（1指向屏幕）
执行>log后， 1----->log (1指向log)
执行2>&1后， 2----->1 (2指向1，而1指向log,因此2也指向了log)
&&
&&与& 区别：两者都表示“与”运算，但是&&运算符第一个表达式不成立的话，后面的表达式不运算，直接返回。而&对所有表达式都得判断。
|| 与| 区别：两者都表示“或”运算，但是||运算符第一个表达式成立的话，后面的表达式不运算，直接返回。而|对所有表达式都得判断。

写得很乱，请原谅一个小白迷茫的絮絮叨叨，若以后有更好的理解，希望能整理成一个简单的说明文档

生信小白学习日记Day4Day5——NGS基础 NGS分析注释（BWA软件）相关推荐

生信小白学习日记Day3——NGS基础 NGS分析注解（质量分析软件）
2019年5月27日,天气舒适,忙碌一天之后开始今天的生信学习.今天就昨天Day2-2的一些标记加以查询说明,仅供参考. NGS基础 NGS分析注解 1. 质量分析软件昨天提到,拿到数据后可以通过一 ...
生信小白学习日记Day2-2——NGS基础 NGS分析
2019年5月26日下午,无意中看到hanli0902的关于NGS分析的博文https://blog.csdn.net/hanli1992/article/details/82790386有很多需要学 ...
生信小白学习日记Day2——NGS基础 illumina高通量测序原理
2019年5月26日,周日,小雨说明:阅读生信宝典和查阅文章的总结,原文请关注公众号生信宝典,参考的博文都附有链接,仅供参考. 生信宝典 NGS基础--高通量测序原理本文介绍了测序文库构建原理.链 ...
生信小白学习日记Day7——WGS分析流程（picard）
2019年6月2日,周日,天气晴,pass 上午.开始学习NGS分析,继BWA比对和SAM文件排序转BAM后的流程. NGS分析 step5 Mark Duplications 参考这篇:GATK使用 ...
生信小白的福音——免费在线分析扩增子数据SILVAngs
文章开头必须感谢一下宏基因组公众号和微信群的各位朋友,平时给予我的温暖和关怀,让我有了写文章的冲动(基情满满). 今天突然听到有个刚刚入坑的同学跟我说,做了60个扩增子不知道怎么分析.What?不会分 ...
210学习日记(18)_ARM基础知识
210学习日记(18) --ARM基础知识注意: 以下大部分类容都来自网上现成的(直接拷贝过来的,然后经整理)!!!! 问1:ARM处理器工作模式有几种?各种工作模式下分别有什么特点? 答1:ARM ...
适合Java零基础小白学习的Java零基础教程
很多Java零基础小白,在刚刚快入门的时候玩命的学习,玩命的记住Java原理,天天早上五点起床背Java的一些英文词汇,然后遇见一些未知的困难,让自己打到癫狂状态,逐渐迷失自我放弃Java,为了解决这 ...
生信小白入门必看网站！常用数据库分享
新手上路,如何快速了解自己课题,含有DNA.RNA序列等信息的核酸数据库肯定是需要了解滴.今天分享一部分大家会比较常用到的,可结合自己的研究背景挑选使用,主要是应用在组学研究or基因功能研究中. 生信 ...
生信漫谈送你超级好用的多序列拼接软件
作为科研人员(生物专业从业者,医学生,临床医生等),接触的各种测序非常多,那么怎么快速的分析你测序的结果呢,测了正向,测了反向序列,是不是测通了,只有拼接在一起才知道,科研真是费脑 https://m ...

生信小白学习日记Day4Day5——NGS基础 NGS分析注释（BWA软件）

NGS基础

NGS分析注释

BWA

生信小白学习日记Day4Day5——NGS基础 NGS分析注释（BWA软件）相关推荐

最新文章

热门文章