[转载][bwa]短序比对工具简介bowtie&vs&BWA&
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）
人 homo sapiens
小鼠 mus musculus
大鼠 rattus norvegicus
果蝇 Drosophila melanogaster
线虫 Caenorhabditis elegans
之后先后出现了两个重要的比对软件， 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版本。但是在的讨论中，他与BWA的作者Heng
Li打起了口水战，他们都声称自己的软件才是准确率最高的。甚至两人还给出了截然相反的两组比较数据（下图）
Wei Shi的比较结果
Heng Li的比较结果
由于公说公有理，婆说婆有理，加上subread还需要时间来考验，所以现在还无从判断谁的更具有优势。
由于选择变多了，人们往往会变得无所适从。就我个人经验而言（其实是对研究机构前人脚本学习和接受），对于ChIP-seq,
RNA-seq，多使用bowtie2，因为它快速，下游结合cufflinks等结果验证率很高。对于SNP， Indels,
methylation分析，使用BWA，下游结合GATK可能会好一点。
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。高通量测序数据的比对小结----bwa、bowtie、bowtie2、samtools
高通量测序数据的比对小结----bwa、bowtie、bowtie2、samtools
1、bwa和samtools的安装BWA的安装1.下载BWA1wget http://sourceforge.net/projects/bio-bwa/files/bwa-0.7.10.tar.bz2/download2.解压1tar jxvf bwa-0.7.10.tar.bz23.编译1cd bwa-0.7.102makeSamtool安裝1.下载Samtool1wget http://sourceforge.net/projects/samtools/files/samtools/0.1.19/tar.bz2/download2.解压Samtool1tar jxvf samtools-0.1.19.tar.bz23.安装1cd samtools-0.1.192make之外Samtool还有一些脚本程序在misc这个文件夹里面最后将bwa和samtools的工作目录加入PATH1export PATH=$PATH:bwa的绝对路径:samtool的绝对路径~/.bashrc2source ~/.bashrc这样就可以直接使用bwa和Samtool了。samtools一般是统计比对信息的，一些处理工作等。bwa的使用：BWA SYNOPSIS语法：1）bwa index ref.fa2）bwa aln ref.fa short_read.fq aln_sa.sai3）bwa samse ref.fa aln_sa.sai short_read.fq aln-se.sam3）bwa sampe ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq aln-pe.sam【bwa bwasw ref.fa long_read.fq aln.sambwa mem ref.fa reads.fq aln-se.sambwa mem ref.fa read1.fq read2.fq aln-pe.sam】bwa的使用需要两种输入文件：Reference genome data（fasta格式 .fa, .fasta, .fna）Short reads data (fastaq格式 .fastaq, .fq)step 1: 建立 Index根据reference genome data(e.g. reference.fa) 建立 Index File1bwa index reference.fabwa index 指令更多的用法及 options，通过以下的命令来查看1bwa indexstep 2: 寻找 SA coordinates如果是pair-end 数据（leftRead.fastq和rightRead.fastq）两个文件分别处理；single-end数据则直接处理。1bwa aln reference.fa leftRead.fastq leftRead.sai2bwa aln reference.fa rightRead.fastq rightRead.sai3bwa aln reference.fa singleRead.fastq singleRead.sai如果希望多线程运行，在其中加入 -t这个参数，另外-f这个参数可以指定结果输出文件，如:1bwa aln -c -t 3 -f leftreads.sai reference.fa leftreads.fastqstep 3:转换SA coordinates输出为sam如果是pair-end数据1bwa sampe -f pair-end.sam reference.fa leftRead.sai rightRead.sai leftRead.fastq rightread.fastq如果是single reads数据1bwa samse -f single.sam reference.fa single.sai single.fastq值此Reads的mapping工作已经完成。接下来用samtools后续处理。2、samtools使用方法名称:&samtools Sequence Alignment/Map (SAM)格式的应用程序语法:samtools view ‐hS -t ref_list.txt ‐o aln.bam aln.samsamtools sort aln.bam aln.sortedsamtools index aln.sorted.bam【samtools view aln.sorted.bam chr2:20,100,000‐20,200,000samtools merge out.bam in1.bam in2.bam in3.bamsamtools faidx ref.fastasamtools pileup ‐f ref.fasta aln.sorted.bamsamtools tview aln.sorted.bam ref.fasta】描述:Samtools是一系列处理BAM格式序列的应用。它从SAM(Sequence Alignment/Map)格式输入或者输出为SAM格式，可以进行排序，合并和建立索引，并且允许快速地检索任意区域的读段（reads）。命令和选项：1）samtoolsview&[‐huHS] [‐ in.refList] [‐o output] [‐f reqFlag] [‐F skipFlag] [‐q minMapQ] [‐l library] [‐
read‐Group] in.bam|in.sam [region1 [...]]，提取/打印所有的或者部分序列文件。如果没有指定区域，所有的序列都将会被打印。否则，仅仅打印覆盖特定区域的序列才会被打印出来。如果一个序列覆盖多个区域，那么它将会被多次打印出来。一个区域可以以如下例格式呈现：`chr2 (the whole chr2), `chr2:1000000 (region starting from 1,000,000bp) or `chr2:1,000,000‐2,000,000 (region between 1,000,000 and 2,000,000bp including the end points).（即染色体名：起始打印位点..终止打印位点）位置是以1开始记录第一个碱基的方式定位（1-based coordinate）。-b 以BAM格式输出，可以用于samtools的后续分析-u 以未压缩的BAM格式输出，可以节约时间，一般在管道执行时使用-h 在结果中包含头header-H 只输出头-S 输入文件为SAM格式，如果确实@SQ头，则需要-t选项-t FILE 这个文件是以TAB分隔的，每行必须包含参考序列名字（如染色体名），一行只包含单个参考序列名字，其它区域被忽略。这个文件定义排序是的参考序列顺序。如果你运行samtools faidx ref.fa，得到的索引结果文件可以作为这里的FILE文件（见6）。-o FILE 输出文件[stdout]，如果在该管道中可以使用command1|tee FILE|command2。-f INT 只输出在FLAG中含有整个INT的序列，INT可以使用十六进制的/^0x[0‐9A‐F]+/ [默认0]-F INT 跳过含有INT的序列-q INT 跳过MAPQ（定位的质量）比INT小的序列[默认0]-l STR 只输出STR库中的序列-r STR 只输出在STR读段组中的序列2）samtoolsimport&in.ref_list in.sam out.bam，从0.1.4起，它与samtools view ‐t in.ref_list ‐o out.bam in.sam意义相同。3）samtoolssort&[‐
] [‐m maxMem] in.bam out.prefix，根据左起位点对序列排序，将产生out.prefix.bam文件，当整个序列无法完全装载到内存（‐m选项控制）时，可能还会产生out.prefix.%d.bam这个临时文件。-n 设定排列方式为读段名称而非染色体位置-m 设定内存使用量[](默认500M)4）samtoolsmerge&[‐h inh.sam] [‐
] out.bam in1.bam in2.bam [...]，将多个排序后的序列文件合并。参考头文件列出所有的输入BAM文件和inh.sam中@SQ头，如果存在，必须全部指向同一系列的参考序列。参考头文件列表（除非被-h选项替代）和in1.bam中的“@”头将会被拷贝到out.bam中，其他的输入文件中的头将会被忽略。-h FILE 使用FILE中的行作为@头拷贝到out.bam中，代替从in1.bam中拷贝的头。（FILE实际上是BAM格式的。但是其中的任何序列信息都会被忽略）。-n 输入文件是以读段名称排序的，而非染色体定位（当sort中使用-n选项时，此处应该选-n）。5）samtoolsindex&aln.bam 对排序后的序列索引，用于快速的随机处理，此处将产生aln.bam.bai文件。6）samtools&faidx&ref.fasta [region1 [...]]，对FASTA格式的参考基因组序列建立索引或者提取已经索引的序列中的子序列（subsequence）。如果没有指定区域，它将会索引文件，并在磁盘上建立ref.fasta.fai，子序列将会被检索并且以FASTA格式打印在标准输出。输入文件可以压缩成RAZF格式。7）samtoolspileup&[‐f in.ref.fasta] [‐ in.ref_list] [‐l in.site_list] [‐iscgS2] [‐T theta] [‐N nHap] [‐
pairDiffRate] in.bam|in.sam，以pileup格式输出序列，在pileup格式中，每一行代表一个基因组位点，包括染色体名称、位点、参考碱基、读段质量和序列定位质量。匹配、错配、插入和缺失、链、mapping质量以及读段的开始和结束位点都在读段碱基一列。在这一列中，点代表与参考序列的正向链匹配，逗号表示与参考序列的反向链匹配，ACTGN代表正向链上的错配，而actgn代表反向链上的错配，+[0‐9]+[ACGTNacgtn]+这一模式表示在参考序列之间存在插入序列，模式中，数字代表插入长度，其后为插入序列。‐[0‐9]+[ACGTNacgtn]+这一模式代表与参考序列相比存在缺失，后面以*代表被删除的碱基。在读段碱基列中，符号^标记一个读段片段的起始，它是读段中被N/S/H CIGAR操作符分隔开临近的子序列（subsequence）。^后面紧跟的ASCII码字符减去33给出mapping质量。$符号标记读段片段（read segment）的结尾。（注意-c选项不同的输出）如果应用-c选项，一致碱基、Phred-scaled一致性质量、SNP质量（指与参考序列相同的Phred-scaled一致性的可能性）和包含该位点的读段mapping质量均方根（root mean square RMS）将会被插入到参考碱基和读段序列之间（即会多出几列）。插入缺失将占有额外一列。每个插入缺失行包含染色体名称、位点、*、基因型、一致性质量、SNP质量、mapping质量均方根，#支持第一种等位序列的读段，#支持第二种等位序列的读段，#与前两种等位形式不同且包含插入缺失的读段。-s 最后一行打印mapping质量。这一选项使得输出结果易于分析，但较耗费空间-S 输入文件时SAM格式-i 只输出含有插入缺失的pileup行-f FILE FASTA格式的参考序列，如果缺少索引文件将产生FILE.fai-M INT设定mapping质量上限（cap）为INT-t FILE以import 命令描述的格式列出参考名称和序列长度列表。如果出现这一个选项，samtools将假设输入文件为SAM格式；否则假设为BAM格式。-l(L) FILE指定pileup输出的位点，第一列为染色体，第二列为1-based位点。其他的列将会被忽略，推荐与-s选项同时使用，因为默认格式我们可能无法知道mapping质量-c 使用MAQ一致性模型检测一致序列，在使用-g和-c选项时，仅-T、-N、-l和-r选项可用-g 产生GLFv3格式的基因型相似性，这一选项抑制-c、-i和-s选项-T FLOAT MAQ一致性检测模式的theta参数[0.85]（错误依赖系数）-N INT样本中的单倍型数量[默认2，要求=2]-r FLOAT 每一对单倍体之间期望的片段差异值[0.001]-I(i) INT 序列中插入缺失的Phred概率[40]8）samtoolstview&in.sorted.bam [ref.fasta]，文本定位查看器，在查看器中，使用？取得帮助，按下g从一个区域开始位置检查定位（alignment），格式如chr10:10,000,000。9）samtools&fixmate&in.nameSrt.bam out.bam，为以名称排序的定位alignment填入配对坐标，ISIZE（inferred insert size猜测的插入序列大小）和配对相应的标签（flag）10）samtoolsrmdup&input.srt.bam out.bam，移除潜在的PCR重复：如果多个读段含有相同的外部位点（external coordinates），只保留那些高mapping质量的碱基对。这一命令只适用于FR开始的并且要求正确设置ISIZE。11）samtoolsrmdupse&input.srt.bam out.bam，除去单端读段（single-ended reads）。这一命令将把所有的读段当作单端读段，即使它们实际上是配对的。12）samtools&fillmd&[‐e] aln.bam ref.fasta，产生MD标签，如果已经存在并且产生的MD标签不同，这个命令将给予警告。-e 如果读段碱基与参考序列相同，将其转换成=，indel检测目前还不支持=作为碱基标志。SAM格式：SAM是以TAB分隔的，除去以@开始的header行，每一定位行（alignment）包含以下项目：|Col | Field | Description |+‐‐‐‐+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+| 1 | QNAME | Query (pair) NAME 读段对应模板名称|| 2 | FLAG | bitwise FLAG 按位的标签（意义见后面）|| 3 | RNAME | Reference sequence NAME 参考序列名称|| 4 | POS | 1‐ased leftmost POSition/coordinate of clipped sequence 1-based左起位置/整齐的序列定位|| 5 | MAPQ | MAPping Quality (Phred‐scaled) 读段定位质量|| 6 | CIAGR | extended CIGAR string 扩展的CIGAR字符串（详见SAM格式详解）|| 7 | MRNM | Mate Reference sequence NaMe (`= if same as RNAME) 配对的参考序列名称，=表示相同|| 8 | MPOS | 1‐ased Mate POSistion 参考序列中1-based配对位置|| 9 | ISIZE | Inferred insert SIZE 猜测的插入序列大小||10 | SEQ | query SEQuence on the same strand as the reference 与参考序列处于同一链的读段序列||11 | QUAL | query QUALity (ASCII‐33 gives the Phred base quality) 度短序列的碱基质量ASCII码-33为碱基质量||12 | OPT | variable OPTional fields in the format TAG:VTYPE:VALUE 变量选项，格式TAG:VTYPE:VALUE |+‐‐‐‐+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+下表列出了FLAG列的含义：+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+| Flag | Description |+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+|0x0001 (1)| the read is paired in sequencing |读段序列是成对的|0x0002 (2)| the read is mapped in a proper pair |读段定位在适当位置|0x0004 (4)| the query sequence itself is unmapped |读段序列自身没有定位|0x0008 (8)| the mate is unmapped |与其配对的读段为定位|0x0010 (16)| strand of the query (1 for reverse) |读段对应链|0x0020 (32)| strand of the mate |配对链|0x0040 (64)| the read is the first read in a pair |读段是读段对的第一个|0x)| the read is the second read in a pair |读段是读段对的第二个|0x)| the alignment is not primary |定位不是最优选|0x)| the read fails platform/vendor quality checks |读段质量未生成|0x)| the read is either a PCR or an optical duplicate |读段是PCR或者光学重复+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+限制：1、 在bam_import.c, bam_endian.h, bam.c和bam_aux.c中的非定位单词2、 CIGAR操作符P不能正确处理3、 合并时，输入文件需要有相同数目的参考基因序列。设备要求可以降低，另外，合并时不会自动重构头字典（header dictionary）。用户必须提供正确的头，Picard做合并表现更好。4、 samtools的rmdup无法处理单端数据，不能移除染色体之间的重复。Picard表现更好。3、bowtie的使用现在举个例子来探讨bowtie的使用方法：现在有GENOME.fa、高通量测序数据Reads.fa，我们希望将Reads.fa比对到基因组GENOME.fa上。（一）、对Reference文件(GENOME.fa)建库1bowtie-build GENOME.fa GENOME.fa建库步骤可能需要1h甚至几个小时，建议在后台执行：nohup bowtie-build GENOME.fa GENOME.fa（二）、将Reads.fa比对到GENOME.fa上，只能比对到正链，且匹配到基因组不多于20个不同位置，允许有1个错配1bowtie -f-a-m20-v1--al Reads_aligned --un Reads_unaligned --norc GENOME.fa Reads.fa Reads.bwt 2 log（三）、bowtie输出结果的说明12sample001_x75 + Chr1 12453 ATCGGCCAATTACGGACTTAA IIIIIIIIIIIIIIIIIIIII 4 9:GT1 2 3 4 5 6 7 8第2列为+时：最左端的数字9表示query从5端数起，第10个碱基为T，而对应的reference为G;第2列为-时：最左端的数字9表示query先作反向互补，然后从3端数起，第10个碱基为T，而对应的reference为G;4、bowtie21）下载（wget http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.3/bowtie2-2.2.3-linux-x86_64.zip/download）解压 unzip bowtie2-。。。。2）Set the BT2_HOME environment variable to point to the new Bowtie 2 directory containing the bowtie2, bowtie2-build and bowtie2-inspect binaries. 【将 其加入环境变量（确保bowtie2,bowtie2-build or bowtie2-inspect能使用】3）建索引：$ bowtie2-build genome.fasta index&$ bowtie2 -p 8 -x index -1 reads1.fq -2 reads2.fq -S out.sam以一个例子详细说明：【【【【Bowtie 2 comes with some example files to get you started. The example files are not scient we use theLambda phage reference genome simply because its short, and the reads were generated by a computer program, not a sequencer. However, these files will let you start running Bowtie 2 and downstream tools right away.First follow the manual instructions to obtain Bowtie 2. Set the BT2_HOME environment variable to point to the new Bowtie 2 directory containing the bowtie2, bowtie2-build and bowtie2-inspect binaries. This is important, as the BT2_HOME variable is used in the commands below to refer to that directory.Indexing a reference genomeTo create an index for the Lambda phage reference genome included with Bowtie 2, create a new temporary directory (it doesnt matter where), change into that directory, and run:$BT2_HOME/bowtie2-build $BT2_HOME/example/reference/lambda_virus.fa lambda_virusThe command should print many lines of output then quit. When the command completes, the current directory will contain four new files that all start with lambda_virus and end with .1.bt2, .2.bt2, .3.bt2, .4.bt2, .rev.1.bt2, and .rev.2.bt2. These files constitute the index - youre done!You can use bowtie2-build to create an index for a set of FASTA files obtained from any source, including sites such asUCSC, NCBI, and Ensembl. When indexing multiple FASTA files, specify all the files using commas to separate file names. For more details on how to create an index with bowtie2-build, see the manual section on index building. You may also want to bypass this process by obtaining a pre-built index. See using a pre-built index below for an example.Aligning example readsStay in the directory created in the previous step, which now contains the lambda_virus index files. Next, run:$BT2_HOME/bowtie2 -x lambda_virus -U $BT2_HOME/example/reads/reads_1.fq -S eg1.samThis runs the Bowtie 2 aligner, which aligns a set of unpaired reads to the Lambda phage reference genome using the index generated in the previous step. The alignment results in SAM format are written to the file eg1.sam, and a short alignment summary is written to the console. (Actually, the summary is written to the standard error or stderr filehandle, which is typically printed to the console.)To see the first few lines of the SAM output, run:head eg1.samYou will see something like this:@HD VN:1.0 SO:unsorted@SQ SN:gi|9626243|ref|NC_| LN:48502@PG ID:bowtie2 PN:bowtie2 VN:2.0.1r1 0 gi|9626243|ref|NC_| M * 0 0 TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACNGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTNTAAAGGCGGCGAGCGCGGCTTTTCCG +@6:27(F5)9)B:%B+A-%5A?2$HCB0B+0=D7E/.03#!.F77@6B==?C7;))%;,3-$.A06+-1/@@?,26=?*@0;$:;??G+:#+(A?9+10!8!?()?7C AS:i:-5 XN:i:0 XM:i:3 XO:i:0 XG:i:0 NM:i:3 MD:Z:59G13G21G26 YT:Z:UUr2 0 gi|9626243|ref|NC_| M * 0 0 NTTNTGATGCGGGCTTGTGGAGTTCAGCCGATCTGACTTATGTCATTACCTATGAAATGTGAGGACGCTATGCCTGTACCAAATCCTACAATGCCGGTGAAAGGTGCCGGGATCACCCTGTGGGTTTATAAGGGGATCGGTGACCCCTACGCGAATCCGCTTTCAGACGTTGACTGGTCGCGTCTGGCAAAAGTTAAAGACCTGACGCCCGGCGAACTGACCGCTGAGNCCTATGACGACAGCTATCTCGATGATGAAGATGCAGACTGGACTGC (#!!+!$%+(+)%)%!+!(++)###!!(%(+%$%*%%#$%#%#!)*(#)(($$%+#%*)*#*%*)(%+!%%*$%#+)$+))*+!*)!*!(*#+(%)*(!$*!!%$$!!(*$#!$%%#)$#+%*+)!*)+(#!)!%*#*)*))($+*%%)!*)!((%+%$###+((!*(($*!*(+)%#$+(**$$+*!#%))(+(!%+ AS:i:-14 XN:i:0 XM:i:8 XO:i:0 XG:i:0 NM:i:8 MD:Z:0A0C0G0A108C23G9T81T46 YT:Z:UUr3 16 gi|9626243|ref|NC_| M * 0 0 GGGCGCGTTACTGGGATGATCGTGAAAAGGCCCGTCTTGCGCTTGAAGCCGCCCGAAAGAAGGCTGAGCAGCAGACTCAAGAGGAGAAAAATGCGCAGCAGCGGAGCGATACCGAAGCGTCACGGCTGAAATATACCGAAGAGGCGCAGAAGGCTNACGAACGGCTGCAGACGCCGCTGCAGAAATATACCGCCCGTCAGGAAGAACTGANCAAGGCACNGAAAGACGGGAAAATCCTGCAGGCGGATTACAACACGCTGATGGCGGCGGCGAAAAAGGATTATGAAGCGACGCTGTAAAAGCCGAAACAGTCCAGCGTGAAGGTGTCTGCGGGCGAT 7F$%6=$:9B@/F=?!D?@0(:A*)7/9C6#16:C(.CC;#.;;2$4D:?B!689?(0(G7+0=@37F)GG=?958.D2E04CE,*AD%G0.%$+A:H;?872:88?E6((CF)6DF#.)=BD-=CB080E5BH77:@70#4%A5=6.2/1;9-H6)=$/0;5E:8G!@::1?2DC7C*;@*#.1C0.DH/20,!C-#,6@%+D(AG-).?.00@)/F8?B!170,):?A7#1(A@0E#A.*DC.E)AH+.,5,252?:G,FD0B8D-6$65DD!A/3B*31?6 AS:i:-22 XN:i:0 XM:i:8 XO:i:0 XG:i:0 NM:i:8 MD:Z:80C4C16A52T23G30A8T76A41 YT:Z:UUr4 0 gi|9626243|ref|NC_| M * 0 0 GGGCCAATGCGCTTACTGATGCGGAATTACGCCGTAAGGCCGCAGATGAGCTTGTCCATATGACTGCGAGAATTAACNGTGGTGAGGCGATCCCTGAACCAGTAAAACAACTTCCTGTCATGGGCGGTAGACCTCTAAATCGTGCACAGGCTCTGGCGAAGATCGCAGAAATCAAAGCTAAGT(=8B)GD04*G%4F,1A.C7=F$,+#6!))43C,5/5+)?-/0/D3=-,2/+.1?@-;)00!3!7BH$G)HG+ADC#-9F)77$?.0)@5;4,!0-#C!15CF8HB+B==H7,/)C5)5*+(F5A%D,EA(G9E07/E?4%;#92)5+@7:A.(BG@BG86@.G AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:77C106 YT:Z:UUr5 0 gi|9626243|ref|NC_| M * 0 0 GTCAGGAAAGTGGTAAAACTGCAACTCAATTACTGCAATGCCCTCGTAATTAAGTGAATTTACAATATCGTCCTGTTCGGAGGGAAGAACGCGGGATGTTCATTCTTCATCACTTTTAATTGATGTATATGCTCTCTT 9%D)A03E1-*7=),:F/0!6,D9:H,9D%:0B(%E,(8EFG$E89B$27G8F*2+4,-!,0D5()=(FGG:5;3*@/.0F-G#5#3-(FDFEG?)5.!)AGADB3?6(@H(:B6!;6G,.?% AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:138 YT:Z:UUr6 16 gi|9626243|ref|NC_| M2D119M * 0 0 TCGATTTGCAAATACCGGAACATCTCGGTAACTGCATATTCTGCATTAAAAAATCAACGCAAAAAATCGGACGCCTGCAAAGATGAGGAGGGATTGCAGCGTGTTTTTAATGAGGTCATCACGGGATNCCATGTGCGTGACGGNCATCGGGAAACGCCAAAGGAGATTATGTACCGAGGAAGAATGTCGCT 1H#G;H$E*E#*)2%66?=9/9=;4)4/@%+5#@#$4A*!D==8#1*A9BA=:(1+#C.#(3#H=9E)AC*5,AC#E536*2?)H14?9B=7(3H/B:+A:8%1-+#(E%$$H5%*CF8!G5B+A4F$7(:?0$?G+$)B-?20F=D!38BH,%=85@+ AS:i:-13 XN:i:0 XM:i:2 XO:i:1 XG:i:2 NM:i:4 MD:Z:72^TT55C15A47 YT:Z:UUr7 16 gi|9626243|ref|NC_| M * 0 0 TCAGCCGGACGCGGGCGCTGCAGCCGTACTCGGGGATGACCGGTTACAACGGCATTATCGCCCGTCTGCAACAGGCTGCCAGCGATCCGATGGTGGACAGCATTCTGCTCGATATGGACANGCCCGGCGGGATGGTGGCGGGG -/@*7A0)2,AAH@%B)*5*23B/,)90.B@%=FE,E063C9?,:26$-0:,.,18494.;FFA;76+5$C:$!A*,B,)@85D%C*:)30@85;?.B$05=@95DCDH53!8G:F:B7/A.E:434 AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:98G21C22 YT:Z:UUThe first few lines (beginning with @) are SAM header lines, and the rest of the lines are SAM alignments, one line per read or mate. See the Bowtie 2 manual section on SAM output and the SAM specification for details about how to interpret the SAM file format.Paired-end exampleTo align paired-end reads included with Bowtie 2, stay in the same directory and run:$BT2_HOME/bowtie2 -x lambda_virus -1 $BT2_HOME/example/reads/reads_1.fq -2 $BT2_HOME/example/reads/reads_2.fq -S eg2.samThis aligns a set of paired-end reads to the reference genome, with results written to the file eg2.sam.Local alignment exampleTo use local alignment to align some longer reads included with Bowtie 2, stay in the same directory and run:$BT2_HOME/bowtie2 --local -x lambda_virus -U $BT2_HOME/example/reads/longreads.fq -S eg3.samThis aligns the long reads to the reference genome using local alignment, with results written to the file eg3.sam.Using SAMtools/BCFtools downstreamSAMtools is a collection of tools for manipulating and analyzing SAM and BAM alignment files. BCFtools is a collection of tools for calling variants and manipulating VCF and BCF files, and it is typically distributed with SAMtools. Using these tools together allows you to get from alignments in SAM format to variant calls in VCF format. This example assumes that samtools andbcftools are installed and that the directories containing these binaries are in your PATH environment variable.Run the paired-end example:$BT2_HOME/bowtie2 -x $BT2_HOME/example/index/lambda_virus -1 $BT2_HOME/example/reads/reads_1.fq -2 $BT2_HOME/example/reads/reads_2.fq -S eg2.samUse samtools view to convert the SAM file into a BAM file. BAM is a the binary format corresponding to the SAM text format. Run:samtools view -bS eg2.sam eg2.bamUse samtools sort to convert the BAM file to a sorted BAM file.samtools sort eg2.bam eg2.sortedWe now have a sorted BAM file called eg2.sorted.bam. Sorted BAM is a useful format because the alignments are (a) compressed, which is convenient for long-term storage, and (b) sorted, which is conveneint for variant discovery. To generate variant calls in VCF format, run:samtools mpileup -uf $BT2_HOME/example/reference/lambda_virus.fa eg2.sorted.bam | bcftools view -bvcg - eg2.raw.bcfThen to view the variants, run:bcftools view eg2.raw.bcf】】】】
发表评论：
TA的最新馆藏