RNA-seq基本介绍及实战——表达定量

实战环节

上一期对所拥有的数据做了个序列比对，现在进入表达定量的环节

表达定量

在进行表达定量的处理之前，需要对原始的比对文件进行处理，这里就有以下步骤：

• 使用picard / samtools
• 将sam格式转换为bam格式
• 对bam文件进行排序
• 去除比对得分较低的序列
• 如果需要可以去除重复reads

在这里将会以三种方法进行表达定量的操作，分别是STAR+RSEM进行表达定量，另外一个就是使用Kallisto进行表达定量操作。最后一个就是使用featureCounts软件进行操作。

STAR+RSEM

这个方法分为两个步骤

• 准定定量分析所需文件
• 利用STAR结果进行定量分析

在进行这个方法之前，需要对RSEM这个软件进行安装

RSEM这个软件的安装方法同样也很简单：

## 下载 RSEM
wget -c https://github.com/deweylab/RSEM/archive/v1.3.1.tar.gz
cd RSEM-1.3.1
## 安装 RSEM
make
make install

接下来构建准备文件，在主目录下创建arab_RSEM文件夹，随后输入以下命令：

rsem-prepare-reference --gtf 00ref/Araport11_GFF3_genes_transposons.201606.gtf 00ref/TAIR10_Chr.all.fasta arab_RSEM/arab_rsem

该命令有三个参数：

• –gtf：标注GTF文件的位置在哪
• 00ref/TAIR10_Chr.all.fasta：该参数是表示参考基因组的位置在哪
• arab_RSEM/arab_rsem：表示输出文件的位置

运行后的输出结果为：

rsem-extract-reference-transcripts arab_RSEM/arab_rsem 0 00ref/Araport11_GFF3_genes_transposons.201606.gtf None 0 00ref/TAIR10_Chr.all.fasta
Parsed 200000 lines
Parsed 400000 lines
Parsed 600000 lines
Parsed 800000 lines
Parsing gtf File is done!
00ref/TAIR10_Chr.all.fasta is processed!
58699 transcripts are extracted.
Extracting sequences is done!
Group File is generated!
Transcript Information File is generated!
Chromosome List File is generated!
Extracted Sequences File is generated!

rsem-preref arab_RSEM/arab_rsem.transcripts.fa 1 arab_RSEM/arab_rsem
Refs.makeRefs finished!
Refs.saveRefs finished!
arab_RSEM/arab_rsem.idx.fa is generated!
arab_RSEM/arab_rsem.n2g.idx.fa is generated!

查看arab_RSEM文件夹可看到生成一批文件：

这些文件表明RSEM的准备工作已经完成。

在准备好文件后，下一步就是调用RSEM进行表达量的计算，在shell中输入以下命令：

rsem-calculate-expression --paired-end --no-bam-output --alignments -p 5 -q 03align_out/sample1_Aligned.toTranscriptome.out.bam arab_RSEM/arab_rsem 04rsem_out/sample1_rsem

接下来解释一下各种参数的含义：

• paired-end：表示序列是双端测序，如果用到双端测序则需要加上这个参数
• no-bam-output：只进行定量的表达计算，不进行bam的输出
• -p：表示运行所用的进程数目
• -q：输入转录本文件
• arab_RSEM/arab_rsem：表示上一步索引文件的位置
• 04rsem_out/sample1_rsem：输出文件的目录以及前缀

在shell中的输出如下：

rsem-parse-alignments arab_RSEM/arab_rsem 04rsem_out/sample1_rsem.temp/sample1_rsem 04rsem_out/sample1_rsem.stat/sample1_rsem 03align_out/sample1_Aligned.toTranscriptome.out.bam 3 -tag XM -q

rsem-build-read-index 32 1 1 04rsem_out/sample1_rsem.temp/sample1_rsem_alignable_1.fq 04rsem_out/sample1_rsem.temp/sample1_rsem_alignable_2.fq

rsem-run-em arab_RSEM/arab_rsem 3 04rsem_out/sample1_rsem 04rsem_out/sample1_rsem.temp/sample1_rsem 04rsem_out/sample1_rsem.stat/sample1_rsem -p 5 -q
Time Used for EM.cpp : 0 h 00 m 10 s

rm -rf 04rsem_out/sample1_rsem.temp

其输出如下图所示:

这里有两个结果，一个是基于基因的结果，另外一个是基于转录本的结果

其中sample1_rsem.genes.results如图所示：

而sample1_rsem.isoforms.results如图所示：

这个就是STAR+RSEM进行表达定量的一个流程

Kallisto分析

这个方法也是分为两个步骤

• 利用转录本参考序列文件构建索引
• 进行无比对定量分析

在运行这个命令之前，需要对kallisto进行安装，同样的，这个软件的安装方法也很简单：

# conda安装
conda install kallisto
# github安装
git clone https://github.com/pachterlab/kallisto.git

cd kallisto

mkdir build
cd build

cmake .. # 这里没有sudo权限的装不上 建议改为：
cmake .. -DCMAKE_INSTALL_PREFIX:PATH=$HOME/kallisto/bin # 或者自己的其他目录

make # 或者make install

安装完毕之后我们运行命令以构建索引：

madir arab_kallisto
cd arab_kallisto
kallisto index -i arab_kallisto ../arab_RSEM/arab_rsem.transcripts.fa

这时候的shell会输出一段过程：

[build] loading fasta file ../arab_RSEM/arab_rsem.transcripts.fa
[build] k-mer length: 31
[build] warning: clipped off poly-A tail (longer than 10)
        from 26 target sequences
[build] warning: replaced 697 non-ACGUT characters in the input sequence
        with pseudorandom nucleotides
KmerStream::KmerStream(): Start computing k-mer cardinality estimations (1/2)
KmerStream::KmerStream(): Start computing k-mer cardinality estimations (1/2)
KmerStream::KmerStream(): Finished
CompactedDBG::build(): Estimated number of k-mers occurring at least once: 54535532
CompactedDBG::build(): Estimated number of minimizer occurring at least once: 13421357
CompactedDBG::filter(): Processed 97555899 k-mers in 58699 reads
CompactedDBG::filter(): Found 54357654 unique k-mers
CompactedDBG::filter(): Number of blocks in Bloom filter is 372803
CompactedDBG::construct(): Extract approximate unitigs (1/2)
CompactedDBG::construct(): Extract approximate unitigs (2/2)
CompactedDBG::construct(): Closed all input files

CompactedDBG::construct(): Splitting unitigs (1/2)

CompactedDBG::construct(): Splitting unitigs (2/2)
CompactedDBG::construct(): Before split: 343803 unitigs
CompactedDBG::construct(): After split (1/1): 343803 unitigs
CompactedDBG::construct(): Unitigs split: 951
CompactedDBG::construct(): Unitigs deleted: 0

CompactedDBG::construct(): Joining unitigs
CompactedDBG::construct(): After join: 319356 unitigs
CompactedDBG::construct(): Joined 24543 unitigs
[build] building MPHF
[build] creating equivalence classes ... 
[build] target de Bruijn graph has k-mer length 31 and minimizer length 23
[build] target de Bruijn graph has 319356 contigs and contains 54404035 k-mers 

后面会生成一个arab_kallisto文件，这个文件就是我们所需要的索引文件

下一步就是运行分析的命令：

kallisto quant -i arab_kallisto/arab_kallisto -o 05kallisto_out/sample1 02clean_data/sample1_paired_clean_R1.fastq 02clean_data/sample1_paired_clean_R2.fastq

接下来解释下面参数的意思：

• quant：quant模式，用来数有多少个COUNT
• -i：后面输入索引文件的位置
• -o：后面是表示输出文件的位置
• 后面的两个参数表示的是清洗后pair的数据

运行过程如下所示：

[quant] fragment length distribution will be estimated from the data
[index] k-mer length: 31
[index] number of targets: 58,699
[index] number of k-mers: 54,404,035
[quant] running in paired-end mode
[quant] will process pair 1: 02clean_data/sample1_paired_clean_R1.fastq
                             02clean_data/sample1_paired_clean_R2.fastq
[quant] finding pseudoalignments for the reads ... done
[quant] processed 196,007 reads, 194,599 reads pseudoaligned
[quant] estimated average fragment length: 155.256
[   em] quantifying the abundances ... done
[   em] the Expectation-Maximization algorithm ran for 766 rounds

输出文件如图所示：

打开abundance.tsv查看内容：

这里就是利用Kallisto进行表达定量的过程。

FeatureCounts

STAR+featureCounts去进行表达定量的方法相较于前两个方法有一个很大的特点，那就是快，其安装方法也非常简单：

# conda安装
conda install subread

在安装之后，定位在03align_out文件夹，使用以下命令运行FeatureCounts

featureCounts -p -a ../00ref/Araport11_GFF3_genes_transposons.201606.gtf -o out_counts.txt -T 6 -t exon -g gene_id sample*_Aligned.sortedByCoord.out.bam

接下来解释一下以下参数的具体含义：

• -p：表示序列是双端测序
• -a：后面明确要求输入gtf文件
• -o：输出运行log文件
• -T：软件运行时所需要调用的线程数目
• -t、-g：通过什么方法进行定量，输出什么结果
• 最后一个是定量时所需要的输入文件

运行后会弹出这个框：

其目录中会产生多两个文件：

其内部写入的信息如下：

表达矩阵

上面用了3个方法进行了表达定量的分析，做完表达定量后，我们还需要对表达定量的结果转换为表达矩阵。

04rsem_out中，在shell中输入以下命令：

rsem-generate-data-matrix *_rsem.genes.results > output.matrix

打开输出的output.matrix查看一下：

可以看到这里的表达矩阵已经做出来了。

下一讲将会聚焦到如何分析差异基因这一步