基于eXpress对转录组和基因组进行量化

General workflow

eXpress是一个通用的丰度估计工具，它可以应用于任意靶序列和高通量测序reads。 靶序列可以是任意基因组区域，例如RNA-seq中的转录本。因此，一般的流程应该是这样的：

1. 选择你要分析的数据

2. 产生靶序列的集合

3.将目的片段比对到靶序列上

4. eXpress需要的参数包括目的片段，然后进行靶序列的丰度估计

5. 额外的下游分析

图示:

这个教程涉及如下工具:

• Gene Expression Omnibus (GEO) - 取得数据
• The Short Read Archive (SRA) -取得数据
• SRA Toolkit - 从压缩的测序数据中抽取FASTQ文件
• UCSC Genome Browser - 取得注释信息
• Bowtie - 进行比对
• Bowtie 2 -进行比对
• eXpress - 靶序列丰度估计

其它有用的工具（不限于RNA-seq）还包括 here.

---

例子: 没有参考基因组序列也没有注释信息的情况

有的时候，你将研究没有参考基因组序列的物种，或者参考序列质量较差。这通常意味着你也没有转录组序列。下面要进行的步骤经常是从头组装转录组。接下来，我们将使用Bowtie2进行片段比对。

取得数据

为了取得一个数据，我们将使用GEO访问号。如果你没有一个GEO访问号，而是仅仅想要浏览数据，你可以跟随这个tutorial。为了演示的目的，我们将要研究牦牛转录组。

为什么选择牦牛呢？因为牦牛是天生没有气味的。事实上，是牦牛毛无味。为了下载数据，就直接去GEO吧，然后在"GEO accession" 输入访问号“GSE33300”，点击“GO”。

这将将你带到主实验页面。可以看到有6个来自不同器官的不同的样品。我们先看一下"GSM823609 brain"。点击这个实验，并点击ftp链接下载SRA文件。点击目录可以看到SRR361433.sra. 这是一个paired end的RNA-seq数据，我们将使用如下命令抽取数据

$ fastq-dump --split-3 SRR361433.sra

结果产生两个文件，SRR361433_1.fastq和SRR361433_2.fastq. 注意到使用 --split-3. 只有当你下载的数据是paired end的情况下才需要用这个参数. 

通过从头组装进行注释

在这个演示中，我们花点时间进行一个完全从头组装分析，而不使用基因组信息，使用的工具是Trinity.
使用如下参数运行Trinity:

$ Trinity.pl --seqType fq --JM 200G --left SRR361433_1.fastq --right SRR361433_2.fastq --CPU 2

在几个小时内，我们将在trinity_out_dir中得到几个文件, 包括注释文件Trinity.fasta.这将是新的注释文件. 从这里开始，如果你有参考基因组的情况下，分析将大致相同（下面的例子也是一样的）.

从这里，可以下载组装文件.

比对

建立索引

在你进行任何比对之前，你首先需要建立靶序列的一系列索引文件.

$ cd trinity_out_dir$ bowtie2-build --offrate 1 Trinity.fasta Trinity

这将在trinity_out_dir中建立索引，base name是Trinity，bowtie2需要的参数只需要写到base name结束即可，不需要后面的部分. 这个索引将允许bowtie2快速将reads比对到靶序列。Offrate 1可以加快比对速度，代价是需要的硬盘空间增大.

进行比对

使用一行命令即可运行bowtie2,

$ bowtie2 -a -X 800 -p 4 -x trinity_out_dir/Trinity \    -1 SRR361433_1.fastq -2 SRR361433_2.fastq | samtools view -Sb - > hits.bam

• -a - 想让bowtie2报告所有的比对可用此参数，适用于转录组分析，不适用于比对到基因组的分析(eXpress将处理多比对的问题，非常慢！)
• -X 800 - 将设置片段长度（fragment length）不超过800（实际的应用一般设置到300以内就可以了）. 这个设置将完全将RNA-seq的测序片段纳入分析的范畴
• -p 4 - 使用4个CPUs进行比对。你应该尽可能多的使用多个CPUs以加快速度.
• -x ... - Bowtie 2索引
• -1 ... -2 ... - RNA-seq实验的左reads和右reads

几十分钟到几小时后(取决于你使用多少CPUs), 你应该看到如下控制台信息:

[samopen] SAM header is present: 165714 sequences.32959992 reads; of these:  32959992 (100.00%) were paired; of these:    4385612 (13.31%) aligned concordantly 0 times    22571641 (68.48%) aligned concordantly exactly 1 time    6002739 (18.21%) aligned concordantly >1 times    ----    4385612 pairs aligned concordantly 0 times; of these:      352368 (8.03%) aligned discordantly 1 time    ----    4033244 pairs aligned 0 times concordantly or discordantly; of these:      8066488 mates make up the pairs; of these:        5942057 (73.66%) aligned 0 times        1356189 (16.81%) aligned exactly 1 time        768242 (9.52%) aligned >1 times90.99% overall alignment rate

运行eXpress

我们现在运行eXpress. 简单的说，你应该准备下面的数据:

• multi-FASTA format的靶参考序列（注释）
• 比对得到的BAM格式的文件

运行eXpress非常简单,

$ express -o xprs_out trinity_out_dir/Trinity.fa hits.bam

参数 -o 建立一个新的文件夹并将结果存储到xprs_out. 下面一行的参数是一个multi-FASTA文件，包含靶序列. 最后，hits.bam是比对过的reads文件.

你将看到如下类似的输出:

Attempting to read 'hits.bam' in BAM format...Parsing BAM header...Loading target sequences and measuring bias background... Processing input fragment alignments...Synchronized parameter tables. Fragments Processed (hits.bam): 1000000          Number of Bundles: 145443Synchronized parameter tables.Fragments Processed (hits.bam): 2000000          Number of Bundles: 144074

最终，你将得到以下文件：

• results.xprs
• params.xprs

在个人电脑中，只是用两核心CPUs, eXpress在30分钟内得到28,613,101个片段的量化计算.

分析results.xprs

results.xprs是一个tab分隔的文件，包括如下列:

• bundle_id - The bundle that the particular target belongs to
• target_id - The name of the target from the multi-FASTA annotation
• length - The true length of the target sequence
• effective length - The length adjusted for biases
• total counts - The number of fragments that mapped to this target
• unique counts - The number of unique fragments that mapped to this target
• estimated counts - The number of counts estimated by the online EM algorithm
• effective counts - The number of fragments adjusted for biases: est_counts * length / effective_length
• alpha - Beta-Binomial alpha parameter
• beta - Beta-binomial beta parameter
• FPKM - Fragments Per Kilobase per Million Mapped. Proportional to (estimated counts) / (effective length)
• FPKM low
• FPKM high
• solveable flag - Whether or not the likelihood has a unique solution

你可以将这个结果文件导入Excel等工具. 你将可能对FPKM进行排序感兴趣。

结果文件可从这里下载。

---

第二个例子：使用已知注释分析大的测序数据

这里，我们将看到一个测序深度非常大的数据. 我们将从Peng et. al.（发表在Nature Biotechnology上）这里看到数据。 一共有 ~520,000,000个reads. 我们将只关注polyA的数据集, 包含402,000,000个reads.

取得数据

如果你已经有一个感兴趣的数据集, 直接将访问号在GEO中输入即可. 这里，作者在文章中提供了SRA的访问号.如果你点击它，或者在SRA搜索该访问号(SRA043767.1) ，你将不会找得到它。将 '.1'去除再搜它，就能找到了. 

这里有很多数据集，但是我们只对polyA的数据感兴趣. 有时候，你可能想通过阅读文件描述决定下载哪些文件，但是这里不行. 为了找到这个信息，我们阅读"Supplemental Information"后，发现了描述数据的页面 (Supplementary Table 1. Overview of sequencing samples).

根据这个表格，我们对HUMwktTBRAAPE_* 和HUMwktTBRBAPE感兴趣. 我们首先点击HUMwktTBRAAPE. 当从SRA下载完后,需要首先确定它是single-end还是paired-end，因为抽取FASTQ的方式不同. 如果我们看到"Spot description"信息, 我们能看出reads是paired-end. 点击size后面的链接, 然后下载SRR324687.sra和SRR325616.sra. 针对HUMwktTBRBAPE也是类似的，相应的SRA文件是SRR324688.sra.

SRA格式的reads是以压缩的方式存储的,但是大多数工具无法直接处理它对应的平展格式FASTQ. 幸运的是，我们可以很容易地抽取FASTQ文件通过使用 fastq-dump程序（来自SRA Toolkit）.

$ fastq-dump --split-3 SRR324687.sra $ fastq-dump --split-3 SRR324688.sra $ fastq-dump --split-3 SRR325616.sra

注意到使用 --split-3. 如果你有paired end的数据，你必须这么用. 等到fastq-dump运行结果后, 你应该有6个.fastq结尾的文件了.

使用已知信息产生靶序列

在你做任何定量工作之前，你应该确定靶序列对于你的分析是不是充足的。对于RNA-seq而言，通常需要GTF格式的转录组注释。因为eXpress需要一般性的靶序列, 必须使用GTF文件为每个转录本生成一个FASTA记录. 

如果你分析的是一个常见的物种，例如人，实际上已经有很多注释了. 如果你分析一个研究较少的物种（例如牦牛），你可能需要自己做从头转录组组装了.

毫无疑问，最简单的建立靶序列的文件的方式是使用已知注释. 如果你研究的物种在UCSC中存在，那更简单了. 对于一个不在UCSC中注释的物种，可以先参考一下附件，这里提供了获得基因组注释的步骤:

1. 访问http://genome.ucsc.edu
2. 点击主菜单的"Tables" 
3. 选择基因组
   • 例如选择 Mammal -> Human -> hg19
4. 选择注释
   • 例如选择 'Ensembl'. 点击: "Ensembl Genes"
5. 改变输出格式"output format"到"sequence"
6. 输入一个具有描述性的名字
   
   • 例如: ensembl-hg19.fa.gz
7. 压缩以节省下载时间
8. 最后的界面看起来应该类似于:
   
   
9. 确保一切都正确, 然后点击"get output". 在下一个界面中, 点击"genomic",然后点击submit. 如果你没看到这些的话，说明你做错了什么! 返回重试!
10. 重要: 在下一个界面中, 确保"Introns"没被选上. 其它就选择默认的就可以了. 点击"get sequence", 你的下载就该开始了. 你的界面应该是类似于此的.
    
    

使用解压缩软件解压下载的文件，然后简单看下.

$ gunzip ensembl-hg19.fa.gz$ head ensembl-hg19.fa>hg19_ensGene_ENST00000237247 range=chr1:66999066-67210057 5'pad=0 3'pad=0 strand=+ repeatMasking=noneAGTTTGATTCCAGAGCCCCACTCGGCGGACGGAATAGACCTCAGCAGCGGCGTGGTGAGGACTTAGCTGGGACCTGGAATCGTATCCTCCTGTGTTTTTTCAGACTCCTTGGAAATTAAGGAATGCAATTCTGCCACCATGATGGAAGGATTGAAAAAACGTACAAGGAAGGCCTTTGGAATACGGAAGAAAGAAAAGGACACTGATTCTACAGGTTCACCAGATAGAGATGGAATTCAGGGGAAAAAAAAGACCCAGAAGACTCAGCTTCTCCTCACCTCTTGCTTCTGGCTCAGAGCCCTCTCGTTAACTCTGTCTCAGAAGAAAAGCAATGGGGCACCAAATGGATTTTATGCGGAAATTGATTGGGAAAGATATAACTCACCTGAGCTGGATGAAGAAGGCTACAGCATCAGACCCGAGGAACCCGGCTCTACCAAAGGAAAGCAC

基于这个multi-FASTA靶序列文件, 你现在可以建立一个bowtie索引了.

比对

建立索引

使用Bowtie,

$ bowtie-build --offrate 1 ensembl-hg19.fa ensembl-hg19

bowtie-build建立了注释的索引，这个索引对于bowtie快速地到靶序列的比对是必须的. Offrate 1可以加快比对速度，代价是需要的硬盘空间增大.

这个索引大概需要25分钟才能建立起来. Ensembl注释通常是比较的注释. 

需要注意的是，你不需要每次都建一次索引，建立一个注释的索引只需要一次，以后直接引用它就可以了. 

进行比对

现在，我们进行实际的比对. 

$ bowtie -aS --chunkmbs 128 -v 3 -X 800 -p 40 /home/lmcb/bowtie_indices/ensembl-hg19/ensembl-hg19 \     -1 SRR324687_1.fastq,SRR324688_1.fastq,SRR325616_1.fastq \    -2 SRR324687_2.fastq,SRR324688_2.fastq,SRR325616_2.fastq \    | samtools view -Sb - > hits.bam

这里使用的参数有些不同，因为使用的是Bowtie:

• -a - 想让bowtie2报告所有的比对可用此参数，适用于转录组分析，不适用于比对到基因组的分析(eXpress将处理多比对的问题，非常慢！)
• -S - 以SAM格式输出
• -v 3 - 允许多达3次错配
• -X 800 - 将设置片段长度（fragment length）不超过800（实际的应用一般设置到300以内就可以了）. 这个设置将完全将RNA-seq的测序片段纳入分析的范畴
• -p 4 - Use四个CPUs进行比对. 越多CPUs越好.
• -1 ... -2 ... - RNA-seq实验 Paired end的FASTQ的左reads和右reads

然后，我们使用管道(|)将输出的SAM直接导入samtools以直接生成BAM文件，可以节省磁盘空间，得到hits.bam.当然，你可以使用附录里提供的脚步简化工作量.

[samopen] SAM header is present: 181648 sequences.# reads processed: 421836549# reads with at least one reported alignment: 184591763 (43.76%)# reads that failed to align: 237244786 (56.24%)Reported 242921711 paired-end alignments to 1 output stream(s)

使用40个核心的CPUs, 这个比对过程在大约11个小时可以完成. 你可能注意到了比对率有些低. 然而，这并不意外，因为测定的主体是中国人群，可能和参考基因组有些不同.

运行eXpress进行丰度估计

因为我们之前运行过的例子使用的是相同的名字，所以eXpress的命令是相同的.

$ express -o xprs_out /home/lmcb/bowtie_indices/ensembl-hg19/ensembl-hg19.fa hits.bam

参数 -o 新建了文件夹病将结果存储在xprs_out.下面一行的参数是一个multi-FASTA文件，包含靶序列. 最后，hits.bam是比对过的reads文件.

你将看到如下类似的输出:

Attempting to read 'hits.bam' in BAM format...Parsing BAM header...Loading target sequences and measuring bias background... Processing input fragment alignments...Synchronized parameter tables.Fragments Processed (hits.bam): 1000000      Number of Bundles: 135798Fragments Processed (hits.bam): 2000000      Number of Bundles: 124769Fragments Processed (hits.bam): 3000000      Number of Bundles: 118656Fragments Processed (hits.bam): 4000000      Number of Bundles: 114476....

最终，你将得到以下文件：

• results.xprs
• params.xprs

在个人电脑中，只是用两核心CPUs, eXpress在120分钟内可运行完.

Analyzing results.xprs

results.xprs是一个tab分隔的文件，包括如下列:

• bundle_id - The bundle that the particular target belongs to
• target_id - The name of the target from the multi-FASTA annotation
• length - The true length of the target sequence
• effective length - The length adjusted for biases
• total counts - The number of fragments that mapped to this target
• unique counts - The number of unique fragments that mapped to this target
• estimated counts - The number of counts estimated by the online EM algorithm
• effective counts - The number of fragments adjusted for biases: est_counts * length / effective_length
• alpha - Beta-Binomial alpha parameter
• beta - Beta-binomial beta parameter
• FPKM - Fragments Per Kilobase per Million Mapped. Proportional to (estimated counts) / (effective length)
• FPKM low
• FPKM high
• solveable flag - Whether or not the likelihood has a unique solution

你可以将这个结果文件导入Excel等工具. 你将可能对FPKM进行排序感兴趣。

结果文件可从这里下载。

---

附录: 有用的脚本

#!/bin/bash# Author: Harold Pimentel if [ $# -eq 0 ]; then    echo "Usage: bowtie2.sh DIR1 [DIR2 DIR3 ... DIRN]" 1>&2    exit 1fi # TODO: Configure the variables in this block!BWT2=`which bowtie2`IDX=/home/lmcb/bowtie2_indices/ensembl-mm9/ensembl-mm9N_THR=40 for DIR in "$@"do    echo "Aligning reads in $DIR..."     if [ ! -e $DIR ]; then        echo "ERROR: Couldn't find directory '${DIR}'" 1>&2        echo "\tskipping...." 1>&2        continue    fi     LEFT=`ls ${DIR}/*_1.fastq | sed ':a;N;$!ba;s/\n/,/g'`    RIGHT=`ls ${DIR}/*_2.fastq | sed ':a;N;$!ba;s/\n/,/g'`     if [ -z "$LEFT"  ] || [ -z "$RIGHT" ]; then        echo "Didn't find matching paired-end reads in $DIR"        echo "Trying single-end..."         SINGLE=`ls ${DIR}/*_1.fastq | sed ':a;N;$!ba;s/\n/,/g'`        if [ -z $SINGLE ]; then            echo "ERROR: Couldn't find any single-end reads either... skipping $DIR"            continue        fi    fi     if [ -z "$SINGLE" ]; then        CMD="$BWT2 -a -X 800 -p $N_THR -x $IDX -1 $LEFT -2 $RIGHT | samtools view -Sb - > ${DIR}/hits.bam"    else        CMD="$BWT2 -a -X 800 -p $N_THR -x $IDX -U $SINGLE | samtools view -Sb - > ${DIR}/hits.bam"    fi     echo "Executing..."    echo $CMD     # Before you run it, you might want to comment the next line so that you can do a "dry run"    $CMD 2>&1 | tee ${DIR}/bwt2.log     unset SINGLEdone

虽然，不能保证的这个脚本在你那儿也能运行成功!, 但是这个脚本对于我能运行成功，并且希望你能使用. 正如脚本快要结束前的注释, 你可能想要注释掉$CMD行,去观察程序的运行的详细过程，并且确保一切如你所愿的运行，然后不需要的话可以再注释上.
如果你将它拷贝并粘贴到一个叫bowtie2.sh的文件中并且使它可执行的话, 你就需要使用如下命令:

$ chmod 755 bowtie2.sh$ ./bowtie2.sh exp1 exp2 exp3

其中 exp1, exp2, exp3 是包含 *.fastq reads的文件夹. 这些reads应该满足假设的左reads *_1.fastq以及右reads的 *_2.fastq的命名习惯.

附录：基于定制注释数据生成multi-FASTA文件

使用UCSC基因组浏览器

使用这个方法，如果:

• 在UCSC没有你感兴趣的注释数据
• 甚至连参考基因组也没有

直接访问UCSC Genome Broswer并再次点击"Tables" . 点击"add custom tracks" ，然后再点击"Choose file"以上传注释. 点击submit. 

上传成功后, 你将进入一个新的界面，然后点击 "go to table browser."

使用从人类数据集中抽取序列的方法来抽取你自己的注释对应的序列.

在本地电脑上使用自定义参考序列

如果你已经安装了TopHat, 其中就包括了一个工具 gtf_to_fasta. 这个工具需要有以下文件,

• multi-FASTA格式的基因组序列
• GFF/GTF格式的注释

使用这个工具也是蛮简单的:

$ gtf_to_fasta ensembl-hg19.gtf hg19.fa ensembl-hg19.fa

其中不同的参数一次代表 [annotation.gtf] [genome.fa] [target out]. 最后得到的靶序列应该看来是这样的:

$ head ensembl-hg19.fa >99061 chr10:134210672-134231367 ENST00000305233ggcggcggcggcggcggcggcggcggcggggcgggggcggcCTGGGACGCGGCGGGAGCATGGAGCCGCGCGCCGGCTGCCGGCTGCCGGTGCGGGTGGAGCAGGTCGTCAACGGCGCGCTGGTGGTCACGGTGAGCTGCGGCGAGCGGAGCTTCGCGGGGATCCTGCTGGACTGCACGAAAAAGTCCGGCCTCTTTGGCCTACCCCCGTTGGCTCCGCTGCCCCAGGTCGATGAGTCCCCTGTCAACGACAGCCATGGCCGGGCTCCCGAGGAGGGGGATGCAGAGGTGATGCAGCTGGGGTCCAGCTCCCCCCCTCCTGCCCGCGGGGTTCAGCCCCCCGAGACCACCCGCCCCGAGCCACCCCCGCCCCTCGTGCCGCCGCTGCCCGCCGGAAGCCTGCCCCCGTACCCTCCCTACTTCGAAGGCGCCCCCTTCCCTCACCCGCTGTGGCTCCGGGACACGTACAAGCTGTGGGTGCCCCAGCCGCCGCCCAGGACCATCAAGCGCACGCGGCGGCGTCTGTCCCGCAACCGCGACC

如果你倾向于使用转录本的名字而不是包含位置信息的索引，你可以使用awk写的命令处理这个文件:

$ awk '{if ($1 ~ /^>/) print ">"$3; else print $0}' ensembl-hg19.fa > ensembl-hg19_fixed_names.fa$ head ensembl-hg19_fixed_names.fa >ENST00000305233ggcggcggcggcggcggcggcggcggcggggcgggggcggcCTGGGACGCGGCGGGAGCATGGAGCCGCGCGCCGGCTGCCGGCTGCCGGTGCGGGTGGAGCAGGTCGTCAACGGCGCGCTGGTGGTCACGGTGAGCTGCGGCGAGCGGAGCTTCGCGGGGATCCTGCTGGACTGCACGAAAAAGTCCGGCCTCTTTGGCCTACCCCCGTTGGCTCCGCTGCCCCAGGTCGATGAGTCCCCTGTCAACGACAGCCATGGCCGGGCTCCCGAGGAGGGGGATGCAGAGGTGATGCAGCTGGGGTCCAGCTCCCCCCCTCCTGCCCGCGGGGTTCAGCCCCCCGAGACCACCCGCCCCGAGCCACCCCCGCCCCTCGTGCCGCCGCTGCCCGCCGGAAGCCTGCCCCCGTACCCTCCCTACTTCGAAGGCGCCCCCTTCCCTCACCCGCTGTGGCTCCGGGACACGTACAAGCTGTGGGTGCCCCAGCCGCCGCCCAGGACCATCAAGCGCACGCGGCGGCGTCTGTCCCGCAACCGCGACC

All better!