bwa的使用方法
bwa的使用需要两中输入文件：step 1: 建立 Index
根据reference genome data(e.g. reference.fa) 建立
新浪博客
bwa的使用方法
bwa的使用需要两中输入文件：step 1: 建立 Index
根据reference genome data(e.g. reference.fa) 建立
Index File
bwa index 指令更多的用法及 options，通过以下的命令来查看
step 2: 寻找 SA coordinates
如果是pair-end 数据（leftRead.fastq和rightRead.fastq）两个文件分别处理
如果希望多线程运行，在其中加入 -t这个参数，另外-f这个参数可以指定结果输出文件，如:
step 3:转换SA coordinates输出为sam
如果是pair-end数据
如果是single reads数据
其他：
fai是对ref基因组文件建的索引，方便软件快速随机读取基因组序列
sai是将fastq比对后出来的文件，用于最后输出比对结果sam文件的
官方文档
http://www.bbioo.com/lifesciences/40-113315-1.html
http://bio-bwa.sourceforge.net/bwa.shtml
[BWA]bwa的使用
P Danecek,
http://www.chenlianfu.com/?p=2103
bwa，即Burrows-Wheeler-Alignment Tool。
BWA 包含 3 种算法：
BWA-bactrack: 用于进行 Illumina reads 的比对。reads 的长度最大为 100bp。 BWA-SW: 用于比对 long-read ，支持的长度为 70bp-1Mbp；同时支持剪接性比对。 BWA-MEM: 最新的算法。和 BWA-SW 的适用性一致，但是更加快速和准确；同时与 BWA-bactrack 相比，在对 70-100bp reads 的比对上有更优的性能。
常用例子和参数如下：
$ bwa index ref.fa -p genome -p STR 输出的数据库前缀。 默认和输入的文件名一致，则输出的数据库在其输入文件所在的文件夹，并以该文件名为前缀。 -a [is|bwtsw] 输入构建 index 的算法。 is 算法简单快速，是默认的选项，但是不能用于基因组大于 2GB 的数据库； bwtsw 适合于大基因组，比如人类基因组。
BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点。这对 long reads 的比对时比较重要的结果。但是却会和 Picard 的 markDuplicates 程序部兼容。
常用例子和参数如下：
$ bwa mem genome reads.fq > aln-se.sam $ bwa mem genome read1.fq read2.fq > aln-pe.sam -t INT 使用的线程数 -p 若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(输入的文件若有2个，则忽略之)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。 -M 加入此参数用于将 shorter split hits 标记为次优，有利于兼容 Picard。
常用例子和参数如下：
$ bwa bwasw genome long_read.fq > aln.sam $ bwa bwasw genome read1.fq read2.fq > aln-pe.sam -t INT 使用的线程数
常用例子：
$ bwa aln genome read1.fq > aln_sa1.sai $ bwa aln genome read2.fq > aln_sa2.sai $ bwa samse genome aln_sa1.sai read1.fq > aln_se.sam $ bwa sampe genome aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam
此条目是由chenlianfu发表在未分类分类目录的。将固定链接加入收藏夹。
[bwa]短序比对工具简介bowtie
http://pgfe.umassmed.edu/ou/archives/3197
有趣的是，大部分的short read比对工具都是由中国人写出来的。因此可以说华大基因（BGI， Beijing Genomics Institute, Chinese Academy of Science）是中国NGS测序技术的摇篮。
速度上较有优势的short read（短序）比对工具最早出现的是SOAP（表1）。它很好地解决了一个问题，那就是如何在小内存（4G）的机器上将短序比对至人类基因组这样的大数据上去。我们都知道，人类基因组的大小为3.2G（表2），光把这样大的数据读入内存都是一件不太容易的事情。所以SOAP对NGS的贡献是值得我们记住的。SOAP在设计之初是针对single-end reads， 所以对paired-end的支持不被大家看好。它的成功也逐步被后起之秀所掩盖。
表1：时代表
表2：基因组大小（数据来源Ensenbl， release 72）
之后先后出现了两个重要的比对软件， MAQ/BWA以及Bowtie。MAQ/BWA是Heng Li发表的。Li是从华大基因走出来的人，后来去了Wellcome Trust
Sanger Institute, 现在在哈佛Broad Institute。MAQ的引用率非常高，并成为了Li的成名作。之后写作的BWA以准确率高而闻名，是SNP分析的首选比对软件。
而Bowtie借着其算法上的优势，在运算速度上一举成名。如果对速度的要求高于准确率的时候，bowtie就成了不二选择。bowtie被广泛地应用于ChIP-seq, RNA-seq的分析当中。
NovoAlign是一款商业软件，但是如果只是科研用途的话，可以直接从其网站上下载到编译好的程序（只支持unix/linux/mac）。它也有MPI版本。但是因为在单机上运行效率问题以及商业化的原因，它的应用并不象BWA和Bowtie那样广泛。
Subread是最新出现的比对软件。作者Wei Shi在文章中声称subread在速度和准确率上都较之前的主流软件有优势。并且它还有R/Bioconductor版本。但是在SEQanswers的讨论中，他与BWA的作者Heng Li打起了口水战，他们都声称自己的软件才是准确率最高的。甚至两人还给出了截然相反的两组比较数据（下图）
由于选择变多了，人们往往会变得无所适从。就我个人经验而言（其实是对研究机构前人脚本学习和接受），对于ChIP-seq, RNA-seq，多使用bowtie2，因为它快速，下游结合cufflinks等结果验证率很高。对于SNP， Indels, methylation分析，使用BWA，下游结合GATK可能会好一点
【vcftools】变异
http://blog.sciencenet.cn/blog-479743-762926.html
The variant call format and VCFtools（http://www.ncbi.nlm.nih.gov/pubmed/21653522）.
http://blog.sciencenet.cn/blog-479743-762926.html
The variant call format and VCFtools（http://www.ncbi.nlm.nih.gov/pubmed/21653522）.
[vcftools]基于全基因组snp数据如何进行主成分分析
1）全基因组snp数据格式为 .vcf
2）利用vcftools软件进行格式转换：vcftools --vcf tmp.vcf --plink --out tmp
此时会生成两个文件：tmp.ped 和 tmp.map
3）利用plink软件进行数据格式转换：./plink --noweb --file tmp --make-bed --out tmp
注意，输入文件和输出文件都不需要文件名的后缀，此时生成3个文件：tmp.bed，tmp.bim 和 tmp.fam
4）利用gcta软件进行pca构建
4.1 ./gcta --bfile tmp --make-grm --autosome --out tmp
此时生成一个文件：tmp.grm.gz
4.2 ./gcta --grm tmp --pca 3 --out pcatmp
此时生成两个文件：pcatmp.eigenval 和 pcatmp.eigenvec
5）将生成的pcatmp.eigenvec用文本编辑器打开，在最上面加入一行：1 2 pc1 pc2 pc3（之间以空格隔开），保存
6）打开R软件
6.1 输入文件：a <- read.table('D:/pcatmp.eigenvec', header=TRUE)
6.2 绘散点图：plot(a$pc1,a$pc2, pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10) , main='pca',xlab='pc1',ylab='pc2')
6.3 添加图例： legend('bottomleft', c('CL','IN','GZ','DA','PP','YN','DX','JY','NP','SL'), pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10))
文件 > 另存为 > Jpeg or Tiff
That's all, Game over. 再次向基因组-health (213256700)予以致谢！
bwa index 指令更多的用法及 options，通过以下的命令来查看
step 2: 寻找 SA coordinates
如果是pair-end 数据（leftRead.fastq和rightRead.fastq）两个文件分别处理
如果希望多线程运行，在其中加入 -t这个参数，另外-f这个参数可以指定结果输出文件，如:
step 3:转换SA coordinates输出为sam
如果是pair-end数据
如果是single reads数据
其他：
fai是对ref基因组文件建的索引，方便软件快速随机读取基因组序列
sai是将fastq比对后出来的文件，用于最后输出比对结果sam文件的
官方文档
http://www.bbioo.com/lifesciences/40-113315-1.html
http://bio-bwa.sourceforge.net/bwa.shtml
[BWA]bwa的使用
P Danecek, http://www.chenlianfu.com/?p=2103
1. 简介
最近需要使用软件GAM-NGS软件，但是该软件貌似不支持bowtie2的结果。可能需要用到BWA，于是参考BWA的网站及其Manual。bwa，即Burrows-Wheeler-Alignment Tool。
BWA 包含 3 种算法：
BWA-bactrack: 用于进行 Illumina reads 的比对。reads 的长度最大为 100bp。 BWA-SW: 用于比对 long-read ，支持的长度为 70bp-1Mbp；同时支持剪接性比对。 BWA-MEM: 最新的算法。和 BWA-SW 的适用性一致，但是更加快速和准确；同时与 BWA-bactrack 相比，在对 70-100bp reads 的比对上有更优的性能。
2. BWA 的下载和安装
$ wget http://jaist.dl.sourceforge.net/project/bio-bwa/bwa-0.7.9a.tar.bz2 $ tar jxf bwa-0.7.9a.tar.bz2 -C /opt/biosoft/ $ cd /opt/biosoft/bwa-0.7.9a/ $ make $ echo 'PATH=$PATH:/opt/biosoft/bwa-0.7.9a' >> ~/.bashrc $ source ~/.bashrc3. BWA用法
2.1 index
在进行 reads 的比对前，需要对 fasta 文件构建 FM-index。常用例子和参数如下：
$ bwa index ref.fa -p genome -p STR 输出的数据库前缀。 默认和输入的文件名一致，则输出的数据库在其输入文件所在的文件夹，并以该文件名为前缀。 -a [is|bwtsw] 输入构建 index 的算法。 is 算法简单快速，是默认的选项，但是不能用于基因组大于 2GB 的数据库； bwtsw 适合于大基因组，比如人类基因组。
2.2 mem
该算法先使用 MEM(maximal exact matches) 进行 seeding alignments，再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点。这对 long reads 的比对时比较重要的结果。但是却会和 Picard 的 markDuplicates 程序部兼容。
常用例子和参数如下：
$ bwa mem genome reads.fq > aln-se.sam $ bwa mem genome read1.fq read2.fq > aln-pe.sam -t INT 使用的线程数 -p 若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(输入的文件若有2个，则忽略之)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。 -M 加入此参数用于将 shorter split hits 标记为次优，有利于兼容 Picard。
2.3 bwasw
对输入的第1个文件的所有序列进行比对。如果输如有 2 个文件，则进行 paired-end 比对，此模式仅对 Illumina 的 short-insert 数据进行比对。在 Paired-end 模式下，BWA-SW依然输出剪接性比对结果，但是这些结果会标记为 not properly paired; 同时如果有多个匹配位点，则不会写入 mate 的匹配位置。常用例子和参数如下：
$ bwa bwasw genome long_read.fq > aln.sam $ bwa bwasw genome read1.fq read2.fq > aln-pe.sam -t INT 使用的线程数
2.4 backtrack
经典的 bwa 先使用 aln 命令将单独的 reads 比对到参考序列，再使用 samse 或 sampe 生成 sam 文件。常用例子：
$ bwa aln genome read1.fq > aln_sa1.sai $ bwa aln genome read2.fq > aln_sa2.sai $ bwa samse genome aln_sa1.sai read1.fq > aln_se.sam $ bwa sampe genome aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam
此条目是由chenlianfu发表在未分类分类目录的。将固定链接加入收藏夹。
[bwa]短序比对工具简介bowtie
http://pgfe.umassmed.edu/ou/archives/3197
有趣的是，大部分的short read比对工具都是由中国人写出来的。因此可以说华大基因（BGI， Beijing Genomics Institute, Chinese Academy of Science）是中国NGS测序技术的摇篮。
速度上较有优势的short read（短序）比对工具最早出现的是SOAP（表1）。它很好地解决了一个问题，那就是如何在小内存（4G）的机器上将短序比对至人类基因组这样的大数据上去。我们都知道，人类基因组的大小为3.2G（表2），光把这样大的数据读入内存都是一件不太容易的事情。所以SOAP对NGS的贡献是值得我们记住的。SOAP在设计之初是针对single-end reads， 所以对paired-end的支持不被大家看好。它的成功也逐步被后起之秀所掩盖。
表1：时代表
| 软件 |
发表年代 |
PMID |
| SOAP |
2008 |
18227114 |
| Maq |
2008 |
18714091 |
| BWA |
2009 |
19451168 |
| Bowtie |
2009 |
19261174 |
| NovoAlign |
2009 |
|
| Subread |
2013 |
23558742 |
| 物种 |
碱基数 |
基因数 |
| 人 homo
sapiens |
3.3G |
～21K |
| 小鼠 mus
musculus |
3.5G |
～23K |
| 大鼠 rattus
norvegicus |
2.6G |
～23K |
| 果蝇 Drosophila melanogaster |
169M |
～14K |
| 线虫 Caenorhabditis elegans |
103M |
～21K |
而Bowtie借着其算法上的优势，在运算速度上一举成名。如果对速度的要求高于准确率的时候，bowtie就成了不二选择。bowtie被广泛地应用于ChIP-seq, RNA-seq的分析当中。
NovoAlign是一款商业软件，但是如果只是科研用途的话，可以直接从其网站上下载到编译好的程序（只支持unix/linux/mac）。它也有MPI版本。但是因为在单机上运行效率问题以及商业化的原因，它的应用并不象BWA和Bowtie那样广泛。
Subread是最新出现的比对软件。作者Wei Shi在文章中声称subread在速度和准确率上都较之前的主流软件有优势。并且它还有R/Bioconductor版本。但是在SEQanswers的讨论中，他与BWA的作者Heng Li打起了口水战，他们都声称自己的软件才是准确率最高的。甚至两人还给出了截然相反的两组比较数据（下图）
由于选择变多了，人们往往会变得无所适从。就我个人经验而言（其实是对研究机构前人脚本学习和接受），对于ChIP-seq, RNA-seq，多使用bowtie2，因为它快速，下游结合cufflinks等结果验证率很高。对于SNP， Indels, methylation分析，使用BWA，下游结合GATK可能会好一点
【vcftools】变异
http://blog.sciencenet.cn/blog-479743-762926.html
The variant call format and VCFtools（http://www.ncbi.nlm.nih.gov/pubmed/21653522）.
http://blog.sciencenet.cn/blog-479743-762926.html
The variant call format and VCFtools（http://www.ncbi.nlm.nih.gov/pubmed/21653522）.
[vcftools]基于全基因组snp数据如何进行主成分分析
1）全基因组snp数据格式为 .vcf
2）利用vcftools软件进行格式转换：vcftools --vcf tmp.vcf --plink --out tmp
此时会生成两个文件：tmp.ped 和 tmp.map
3）利用plink软件进行数据格式转换：./plink --noweb --file tmp --make-bed --out tmp
注意，输入文件和输出文件都不需要文件名的后缀，此时生成3个文件：tmp.bed，tmp.bim 和 tmp.fam
4）利用gcta软件进行pca构建
4.1 ./gcta --bfile tmp --make-grm --autosome --out tmp
此时生成一个文件：tmp.grm.gz
4.2 ./gcta --grm tmp --pca 3 --out pcatmp
此时生成两个文件：pcatmp.eigenval 和 pcatmp.eigenvec
5）将生成的pcatmp.eigenvec用文本编辑器打开，在最上面加入一行：1 2 pc1 pc2 pc3（之间以空格隔开），保存
6）打开R软件
6.1 输入文件：a <- read.table('D:/pcatmp.eigenvec', header=TRUE)
6.2 绘散点图：plot(a$pc1,a$pc2, pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10) , main='pca',xlab='pc1',ylab='pc2')
6.3 添加图例： legend('bottomleft', c('CL','IN','GZ','DA','PP','YN','DX','JY','NP','SL'), pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10))
文件 > 另存为 > Jpeg or Tiff
That's all, Game over. 再次向基因组-health (213256700)予以致谢！