链接地址：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 […]
Recent Comments