小生这厢有礼了(BioFaceBook Personal Blog)

getting started with Bowtie 2: Lambda phage example-从这里开始使用Bowtie2：λ噬菌体的例子

Bowtie2自带了一些入门级的示例文件，这些示例文件并不具有科学含义，我们用λ噬菌体的参考基因组只是因为它很短，并且例子里面的reads是由一个电脑程序生成的而不是测序的结果。但是，这些文件能让你立即开始运行Bowtie2和下游的程序。

首先按照获取Bowtie2的指导下载它。设置Bowtie2环境变量BT2_HOME，把它指向含有bowtie2, bowtie2-build和bowtie2-inspect二进制文件的新Bowtie2的文件夹。这一步很重要，因为在下面的命令当中，变量BT2_HOME被用来代表那个文件夹的位置。

Indexing a reference genome-对参考基因组建立索引

为了对Bowtie2内置的λ噬菌体的参考基因组建索引，先新建一个临时文件夹（建在哪里无所谓），进入那个文件夹，然后运行：

$BT2_HOME/bowtie2-build $BT2_HOME/example/reference/lambda_virus.fa lambda_virus

这条命令在结束前应该会打印很多行输出。当其运行完毕时，当前文件夹会产生4个新的文件，它们的文件名都以lambda_virus开始，分别以.1.bt2, .2.bt2, .3.bt2, .4.bt2, .rev.1.bt2和.rev.2.bt2结束。这些文件构成了索引——你完成了！

你可以使用bowtie2-build对一组任意来源的FASTA文件构建索引，包括像UCSC，NCBI，和Ensembl这些站点。当对多个FASTA文件建立索引时，你要在指定所有的文件，并用逗号隔开。更多关于如何用bowtie2-build建立索引的信息，请查看使用手册的建立索引部分。你可能也会直接获取一个已经建好的索引，从而绕过这一步。使用已建好的索引给出了例子。

Aligning example reads-比对示例reads

在上一步创建的文件夹中，现在含有lambda_virus的索引文件。接下来，运行：

$BT2_HOME/bowtie2 -x lambda_virus -U $BT2_HOME/example/reads/reads_1.fq -S eg1.sam

这个命令会运行Bowtie2的比对软件，它会使用上一步建立的索引，把一组非双端测序的reads比对到λ噬菌体的参考基因组上。这步比对的结果是SAM格式的，输出文件是eg1.sam，同时比对的总结会被输出到终端控制台。（事实上，总结是被写进了“standard error” 或 “stderr”，即标准错误句柄里面，通常它会被输出到终端。）

要查看SAM结果的前几行，运行：

head eg1.sam

你会看到类似于下面的东西：

@HD VN:1.0 SO:unsorted @SQ SN:gi|9626243|ref|NC_001416.1| LN:48502 @PG ID:bowtie2 PN:bowtie2 VN:2.0.1 r1 […]

The summary of tophat-cufflinks protocol is like that:

step1: generate a tophat_out folder with bam files tophat -G genes.gtf <index> sample1_1.fq sample1_2.fq tophat -G genes.gtf <index> sample2_1.fq sample2_2.fq step2: generate new .gtf files (assemble isoform) cufflinks sample1/accepted_hits.bam cufflinks sample2/accepted_hits.bam step3: prepare a text file named assemblies.txt with following gtf files cat << EOF > assemblies.txt […]

链接地址：http://blog.sciencenet.cn/blog-635619-884213.html

人的基因组一共有两万多个基因，但是这些基因不是每时每刻都在表达，在不同发育时期和不同组织中，基因的表达是不同的，一个检测这些表达的有效的方法就是RNA-seq，它结合了下一代测序的技术来对细胞整个的mRNA进行测序，从而确定每一个基因的表达量和表达区段，主要用在分析不同条件下细胞内基因表达差异和分析基因表达的不同可变剪接上。

RNAseq分析大致分下面几个步骤，首先要把测到的序列map到基因组上，然后根据map到的区段对细胞构建转录本，然后比较几种细胞的转录本并且合并，最后衡量差异和可变间接和其他的分析。

Mapping

所有的序列分析的第一步，都是把测到的序列map到基因组上，这样就能知道序列原来是在基因组的什么地方。mapping一般基于两种快速索引算法，一种是哈希，MOSAIK，SOAP，SHRiMP用的就是这种算法，在对参照基因组建好哈希表之后，可以在常数次的运算里查找到给定序列的位置，虽然高效，但是由于基因组有些区段重复性很高，所以查找次数虽是常数，但有时会变得非常大，降低效率；还有一种叫Burrows-Wheeler变换，BWA，Bowtie 和SOAP2都是用它，Burrows-Wheeler变换的设计比哈希更加巧妙，它最开始是一种文本压缩算法，文本重复性越高，它的压缩比就越大，这正好克服了基因组重复性高的问题，而且对于一个精确的序列查找，最多在给定序列的长度的次数里就能找到匹配，所以说基于Burrows-Wheeler变换的软件在mapping里用得更加广泛。可是RNAseq的map还有另外一个问题，那就是要允许可变剪接的存在，因为一条RNA不一定是一个外显子表达出来的，也有可能是几个外显子结合在了一起，原来基因里的内含子被空了出来，这些内含子的长度从五十到十万个碱基不等，如果直接用DNAseq的方法的话去在基因组里寻找，有些正好在两个exon连接处的序列就会有错配，而且有些在进化过程中遗漏下来的假基因是没有intron的，这样就导致有些序列会被map到假基因上去，使假基因的表达变得很高，所以，传统的bwa和bowtie在RNAseq里都不是最好的选择。

更加适合RNA mapping的软件需要克服上面的两个问题，Tophat，subread, STAR, GSNAP, RUM, MapSplice都是为RNA测序而开发的，我只用过Tophat，它的新版Tophat2，在map的过程中分三个步骤，如果基因注释文件存在的话，它会先用注释文件的转录组来map，然后再对剩下的序列用bowtie进行普通的map，最后再用bowtie里用过的所有的序列做剪接map，所以跟其他的软件比起来会有比较高的正确率。在运行tophat之前，要对参考基因组作index，samtools可以轻松搞定。

samtools faidx hg38.fasta

如果用的是tophat的话，还需要用bowtie2做index.

bowtie2-build genome.fa genome

hg38.fasta是人类的参考基因组，对参考基因组做index是为了提高mapping时查找的效率.

(对懒人来说，在这里可以直接下载到Bowtie打包做好的index文件，这样就可以省略掉做index的步骤了 )

然后用Tophat开始mapping：

path/tophat

-p 8

-G $path_ref/Homo.GTF

-o tophat_output

$path_ref/hg38

single_end.fastq;

-p指定用几个线程来工作，-G指定注释文件的位置，-o指定输出文件的路径和文件名，最后两个参数分别告诉参考基因组的位置和要map的fastq文件。在参考基因组里不用加.fa的后缀，因为程序还要去寻找其他后缀的index文件。

这步的运行时间是根据fastq文件的大小和设定的线程数来决定的，一般单端8个线程需要每G一个小时，双端各4G，线程数设为16的话需要五六小时，运行完之后会在fastq文件的目录里产生一个-o命令指定的文件夹，这个文件夹里有几个bam文件和bed文件，还有一个summary，在下一步需要用到的是accepted_hits.bam这个文件。

注释文件的后缀是GTF，它包含了所有已知的基因的外显子在基因组中的位置。（hg38的注释文件可以在这里下载）所以对于已经map在基因组上的序列，我们可以直接根据它的位置从注释文件里查找它是不是属于一个外显子，或者是一个转录本。对于要不要在这一步提供注释文件，各有各的看法，我用单端测序的序列用两种方法做了实验，发现他们有些差别：

有注释的时候：

没注释的时候：

发现没注释的时候有更多的序列被map上去，但是重复的map也变多了，但是既然Tophat2三步map的第一步是根据注释文件来map，我还是觉得在运行tophat的时候用注释文件比较好。

Mapping的最后一步是去除map到基因组中多于一处的序列，如果出现好几个序列都map在完全相同的一个区段，那么就应该只保留一个这样的序列，所以，只保留匹配最高的那一个。而且这样的序列占很大一部分，这步也很简单，samtools里的rmdup可以轻松解决：

samtools rmdup -s input.bam output.bam

-s小写是告诉samtools，bam文件是单端测序的结果，不指定-s的话默认是双端。

2。构建转录本

Mapping完了以后，cufflinks就可以把map到基因组里的序列组装成一个转录组了，这个转录组理论上包含了所有当时细胞里的所有mRNA，组装好的转录组包含了可能的剪切信息和所有转录的表达量，这个表达量是根据map到基因组的序列的总数和每个转录片断的长度进行归一化的，听起来比较难懂，它是对于在转录片断里的每一千个碱基对，在每一百万个成功map的序列中，map在这一千个碱基对上的序列的比例，fragments per kilobase of transcript per million mapped […]

(This is now obsolete. The use of the Blast2GO Command Line is recommended for this task.)

This page describes how to create a local installation of the B2G MySQL database and is intended for users with knowledge in MySQL databases and general server administration tasks. If you are not able to follow the steps given […]