2.1.Expression Matrix

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

1) Files Needed

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--"的比例几乎相同,有很大的把握认定这个数据是由非链特异性建库得到的。

结果的解释详情可参考RSeQCinfer_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的数值,导致绝大部分样本看起来颜色差异不明显的情况。

Last updated