2.1.Expression Matrix

本节教学如何使用featureCounts统计落在每个基因的exons上的reads数量，以及如何将这一步的结果整理成表达矩阵。

1) Files Needed

• 本节需要大家下载一个新的docker镜像bioinfo_expmat-3.0.tar.gz，加载镜像的方式和之前相同:

docker load -i bioinfo_expmat-3.0.tar.gz  #load the image

docker run --name=bioinfo_featurecount -dt -h featurecount_docker --restart unless-stopped -v C:\Users\usrname\Desktop\bioinfo_tsinghua_share:/home/test/share bioinfo/expmat:3.0 /bin/bash

docker exec -it bioinfo_featurecount /bin/bash

cd /home/test

• 本节中我们会用到两个bam文件也可以直接从 **** Files needed 中的Files/ 路径下的相应文件夹中下载。

2) Count reads mapped to a gene

2a) Sequencing protocols

测序方法分为非链特异性测序方法和链特异性测序方法两种(Strand specific and nonspecific sequencing protocols)。非链特异性测序方法得到的reads没有方向性，无法判断reads是属于Gene A还是属于Gene B(Antisense RNA)的。链特异性测序方法得到的reads具有方向，可以根据reads方向判断reads是属于Gene A还是属于Gene B。

• 非链特异性测序方法的基本流程如下。

测序得到的reads1.fastq和reads2.fastq没有方向性，因此我们将mapping到Gene A的所有reads都归为Gene A的reads。

• 链特异性测序方法的基本流程如下。

链特异性测序方法根据reads1.fastq的方向性又分为forward和reverse两种，更多细节请参考课件 PPT。

2b) Count reads/counts

对不同的sequencing protocols，计算每个基因（如Gene A）上的mapped reads的方法也不同。

2c) How to tell which sequencing protocol being used

为了准确的计算出每个基因的mapped reads，我们需要知道测序方法的方向性。主要有以下途径。

• 资料参考：当明确知道测序方法，可以通过经验或查阅资料得知不同测序方法的方向性。

• 软件判别：当数据测序方法未知时，我们可以通过RseQC软件判断测序方法的方向性。RSeQC中的infer_experiment.py命令可以用来确定连特异性，执行该功能时需要用到需要确定的RNA-seq数据mapping后得到的.SAM/.bam文件作为输入以及参考基因组.bed(12列)。具体操作如下:

cd /home/test
/usr/local/bin/infer_experiment.py -r GTF/Arabidopsis_thaliana.TAIR10.34.bed -i bam/Shape01.bam

输出为

Reading reference gene model GTF/Arabidopsis_thaliana.TAIR10.34.bed ... Done
Loading SAM/BAM file ...  Total 200000 usable reads were sampled


This is PairEnd Data
Fraction of reads failed to determine: 0.0231
Fraction of reads explained by "1++,1--,2+-,2-+": 0.4792
Fraction of reads explained by "1+-,1-+,2++,2--": 0.4977

从结果中我们可以看到，"1++,1--,2+-,2-+"与"1+-,1-+,2++,2--"的比例几乎相同，有很大的把握认定这个数据是由非链特异性建库得到的。

结果的解释详情可参考RSeQC的infer_experiment.py部分。

3) featureCounts

3a) Count

cd /home/test

/home/software/subread-2.0.3-source/bin/featureCounts \
-s 0 -p -t exon -g gene_id \
-a GTF/Arabidopsis_thaliana.TAIR10.34.gtf \
-o result/Shape01.featurecounts.exon.txt bam/Shape01.bam

featureCounts参数说明：

• -s: 是指输入的bam文件的方向性。0代表非链特异测序；1代表reads1.fastq的方向性分为forward；2代表reverse。

• -a: GTF/GFF基因组注释文件

• -p: 表明输入的bam文件是paired end数据处理得到的

• -F: 指定-a注释文件的格式，默认是GTF

• -g: 从注释文件中提取Meta-features信息用于read count，默认是gene_id，取值参考gtf第9列注释信息。

• -t: 从注释文件中提取某部分信息用于read count，其是默认将exon作为一个feature，取值范围是gtf文件中的第3列的值，可选CDS等。

• -o: 输出文件前缀

• -T: 线程数

对于paired end数据，给定-p参数后，featureCounts早期版本中默认把一个read pair，即测序文库中的一个RNA fragment count一次。在我们这里用到的subread-2.0.3中，默认按reads的数量来count，即一个read pair会被count两次；如果希望对fragment进行count，还需要额外加上一个--countReadPairs参数。 这两种统计方法都是合理的(我们前面给出的read counting的例子对应新版本的默认行为)，但是大家在使用前还是应当注意看一下当前使用的版本的帮助文档(featureCounts -h)。

3b) Generate count matrix for multiple samples

• 这里我们只提供了一个样本，实际情况中大家可能会需要将多个样本产生的count合并成一个表达矩阵（expression matrix）。

• 表达矩阵一般行为基因，列为样本。featureCount的输出前两行为注释信息，第一列为gene id,第7列为count值。大家可以用R语言中读取每一个样本的表达值，再合并成表达矩阵。

4) Homework

• 1) 请阐述 RNA-seq 中归一化基因表达值的几种基本计算方法。

• 2) 根据下述图片描述，填出对应选项:

• 3) 通过软件计算，判断给出文件shape02数据是来自哪一种sequencing protocols （strand nonspecific, strand specific - forward, strand specific - reverse)，并选择合适的参数计算shape02的read count matrix，给出AT1G09530基因(PIF3基因)上的counts数目。

• 4) tumor-transcriptome-demo.tar.gz提供了结肠癌(COAD)，直肠癌(READ)和食道癌(ESCA)三种癌症各50个样本的bam文件用featureCount计算产生的结果。请大家编写脚本将这些文件中的counts合并到一个矩阵中(行为基因，列为样本), 计算logCPM的Z-score，并用 heatmap 展示，提供代码和heatmap。根据heatmap可视化的结果，你认为这三种癌症中哪两种癌症的转录组是最相似的?

• 为了大家在本地机器上计算方便，我们只从featureCount的结果中随机选择了2000个基因。在实际中我们会把所有基因都用上(对于人类如果使用gencode的注释大约是60000个左右)，或只考虑蛋白编码基因(人类有20000个左右)

• 可使用tar xvzf tumor-transcriptome-demo.tar.gz解压文件

• 在R语言中可以用下面的方式计算CPM:

CPM.matrix <- t(1000000*t(counts.matrix)/colSums(counts.matrix))
log10.CPM.matrix <- log10(CPM.matrix+1) # 1 为pseudocount, 避免count为0时对数未定义的情况 

• 除此之外利用R包edgeR的cpm函数计算也是可以的:

# 如果没有安装edgeR，可通过BiocManager::install("edgeR")安装
library(edgeR)
# count.matrix 行为基因，列为样本，数值为counts
y <- DGEList(counts = count.matrix) # 定义edgeR用于存储基因表达信息的DGEList对象
CPM.matrix <- edgeR::cpm(y,log=F) # 计算CPM
log10.CPM.matrix <- log10(CPM.matrix+1) # 1 为pseudocount, 避免count为0时对数未定义的情况 

• 一个基因在某样本z-score为该基因在该样本表达量减去该基因平均表达量，除以该基因表达量的标准差

z.scores <- (log10.CPM.matrix - rowMeans(log10.CPM.matrix))/apply(log10.CPM.matrix,1,sd)
# apply(log10.CPM.matrix,1,sd)表示计算每行(1表示行,2表示列)的标准差(sd函数)
# rowMeans(log10.CPM.matrix)和apply(log10.CPM.matrix,1,mean)效果是一样的

• 作图建议使用pheatmap package,可参考 I.3.2 Plot with R 的例子。可视化时可以将>2的z-score clip到2，<-2的z-score类似的clip到-2，避免为了展示出个别outlier的数值，导致绝大部分样本看起来颜色差异不明显的情况。