理论基础
在上篇的比对中，我们需要纠结是否真的需要比对，如果你只需要知道已知基因的表达情况，那么可以选择alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）。到了这一篇的基因（转录本）定量，需要考虑的因素就更加多了，以至于我不知道如何说清才能理清逻辑。
定量分为三个水平
- 基因水平(gene-level)
- 转录本水平(transcript-level)
- 外显子使用水平(exon-usage-level)。
在基因水平上，常用的软件为HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。以常用的HTSeq-count为例，这些工具要解决的问题就是根据read和基因位置的overlap判断这个read到底是谁家的孩子。值得注意的是不同工具对multimapping reads处理方式也是不同的，例如HTSeq-count就直接当它们不存在。而Qualimpa则是一人一份，平均分配。
image
对每个基因计数之后得到的count matrix再后续的分析中，要注意标准化的问题。如果你要比较同一个样本(within-sample)不同基因之间的表达情况，你就需要考虑到转录本长度，因为转录本越长，那么检测的片段也会更多，直接比较等于让小孩和大人进行赛跑。如果你是比较不同样本（across sample）同一个基因的表达情况，虽然不必在意转录本长度，但是你要考虑到测序深度（sequence depth)，毕竟测序深度越高，检测到的概率越大。除了这两个因素外，你还需要考虑GC%所导致的偏差，以及测序仪器的系统偏差。目前对read count标准化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之间的差异与换算方法已经有文章进行整理和吐槽了。但是，有一些下游分析的软件会要求是输入的count matrix是原始数据，未经标准化，比如说DESeq2，这个时候你需要注意你上一步所用软件会不会进行标准化。
在转录本水平上，一般常用工具为Cufflinks和它的继任者StringTie， eXpress。这些软件要处理的难题就时转录本亚型（isoforms）之间通常是有重叠的，当二代测序读长低于转录本长度时，如何进行区分？这些工具大多采用的都是expectation maximization（EM）。好在我们有三代测序。
上述软件都是alignment-based，目前许多alignment-free软件，如kallisto, silfish, salmon，能够省去比对这一步，直接得到read count，在运行效率上更高。不过最近一篇文献[1]指出这类方法在估计丰度时存在样本特异性和读长偏差。
在外显子使用水平上，其实和基因水平的统计类似。但是值得注意的是为了更好的计数，我们需要提供无重叠的外显子区域的gtf文件[2]。用于分析差异外显子使用的DEXSeq提供了一个Python脚本（dexseq_prepare_annotation.py）执行这个任务。
小结
计数分为三个水平： gene-level, transcript-level, exon-usage-level
标准化方法： FPKM RPKM TMM TPM
输出表达矩阵
在RNA-Seq分析中，每一个基因就是一个feature（特征？），而基因被认为是它的所有外显子的和集。在可变剪切分析中，可以单独把每个外显子当作一个feature。而在ChIP-Seq分析中，feature则是预先定义的结合域。但是确定一个read到底属于哪一个feature有时会非常棘手。因此HTSeq提供了三种模式，示意图见前一幅图
- the union of all the sets S(i) for mode union. This mode is recommended for most use cases.
- the intersection of all the sets S(i) for mode intersection-strict.
- the intersection of all non-empty sets S(i) for mode intersection-nonempty.
基本用法非常的简单：
htseq-count [options] <alignment_file> <gtf_file>
# 安装
conda install htseq
# 利用htseq-count对sort之后的bam文件进行reads计数
htseq-count -r pos -f bam /home/cenhui2018/QWJ/sequence_data/20191030_NGS_DATA/19R577_paired.hisat2_sorted.bam /home/cenhui2018/QWJ/sequence_data/genome/Anotation/gencode.vM23.annotation.gtf > 19R577_paired.count 2>/home/cenhui2018/QWJ/sequence_data/20191030_NGS_DATA/19R577_htseq.log &
运行的时间会比较久，所以可以去了解不同参数的用法了，其中比较常用的为：
- -f bam/sam： 指定输入文件格式，默认SAM
- -r name/pos: 你需要利用samtool sort对数据根据read name或者位置进行排序，默认是name
- -s yes/no/reverse: 数据是否来自于strand-specific assay。DNA是双链的，所以需要判断到底来自于哪条链。如果选择了no， 那么每一条read都会跟正义链和反义链进行比较。默认的yes对于双端测序表示第一个read都在同一个链上，第二个read则在另一条链上。
- -a 最低质量， 剔除低于阈值的read
- -m 模式 union（默认）, intersection-strict and intersection-nonempty。一般而言就用默认的，作者也是这样认为的。
- -i id attribute: 在GTF文件的最后一栏里，会有这个基因的多个命名方式（如下）， RNA-Seq数据分析常用的是gene_id， 当然你可以写一个脚本替换成其他命名方式。
对count数据进行查看
# 利用wc指令我们可以计算count文件的行数
> wc -l *.count
运行结果：
统计count文件的行数
看下每个文件的格式，
# 查看前4行,第一列ensembl_gene_id,第二列read_count计数
> $ head -n 4 *.count
# 查看后四行
> $ tail -n 4 *.count
运行结果：
查看count文件的首行及末行
合并表达矩阵并进行注释（R中进行）
上一步得到的2个单独的矩阵文件，现在要把这2个文件合并为行为基因名，列为样本名，中间为count的矩阵文件。
-
从上面看出需要至少做两步工作才能更好理解和往下进行分析
第一，需要把2个文件合并；
第二，需要把ensembl_gene_id转换为gene_symbol;（这一步不进行也行，后面还需要）
运行R设置工作路径并查看当前目录下的文件
> setwd("/home/cenhui2018/QWJ/sequence_data/20191030_NGS_DATA/")
> list.files()
运行结果
(1) 载入数据，添加列名
# 防止R自动把字符串string的列辨认成factor
> options(stringAsFactor = FALSE)
# 从19R576_paired.count文件中读取数据，并添加列名，之后生成新文件命名为control_W55
> control_W55 <- read.table("19R576_paired.count",sep = "\t",col.names=c("gene_id","control_W55"))
# 查看前6行
> head(control_W55)
# 类似的方法处理19R577_paired.count，并存储到
> test_K54 <- read.table("19R577_paired.count",sep = "\t",col.names=c("gene_id","test_K54"))
运行结果：
(2)数据整合
merge进行整合
gencode的注释文件中的gene_id（如ENSMUSG00000105298.13_3）在EBI是不能搜索到的，所以用gsub功能只保留ENSMUSG00000105298这部分
处理之前先看一下,也就是最后5行是我们不需要的，可以删除
查看文件的最后五行
整合文件：
# 将control_w55和test_K54这两个文件通过"gene_id"进行merge，生成新的merge好的文件raw_count
> raw_count <- merge(control_W55,test_K54, by = "gene_id")
# 查看merge后的文件的前6行
> head(raw_count)
# 查看后六行
> tail(raw_count)
运行结果：
合并文件并查看前六行
查看合并后的文件的后六行
删除前5行
这里要注意，因为读入之后顺序变了，删除的时候看下删除的是哪些行
> raw_count_filt <- raw_count[-1:-5,]
> head(raw_count_filt)
删除之后的结果：
因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。
# 第一步通过gsub工具将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
>ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id)
# 将ENSEMBL重新添加到raw_count_filt1矩阵
>row.names(raw_count_filt) <- ENSEMBL
# 查看>head(raw_count_filt)
执行结果：
可以将最后得到的raw_count_filt文件通过save()函数保存成.Rdata文件，然后复制放在自己的桌面，然后用自己的Rstudio打开
(3) 对基因进行注释-获取gene_symbol
以下两种方式可以进行
第一：去这里或这里的网页版，输入列表即可输出，不再赘述
第二：用bioMart对ensembl_id转换成gene_symbol
bioMart包是一个连接bioMart数据库的R语言接口，能通过这个软件包自由连接到bioMart数据库
# 首先检查BioManager包是否安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# 如果BioManager包已经安装，则通过BioManager包安装biomaRt包
BiocManager::install("biomaRt")
# 因为我的文件是放在桌面上的，所以先设置R的工作目录为桌面路径
> setwd("C://Users/My/Desktop/")
# 加载放在桌面上的.Rdata文件
> load("raw_count.Rdata")
# 加载所需要的R包
> library('biomaRt')
> library("curl")
# 用bioMart对差异表达基因进行注释
> mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
> my_ensembl_gene_id<-row.names(raw_count_filt)
> mms_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name',"description"),filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
执行完之后需要等待一段时间提交检索，下方会有一个进度条
基因注释结束后，查看注释结果
基因注释结束后的结果
合并数据成一个文件
注意：合并的话两个数据必须有共同的列名，我们先看一下
可见，两个文件没有共同的列名，所以要先给'raw_count_filt'添加一个‘ensembl_gene_id’的列名
方法如下:
# 将raw_count_filt的行名赋值到一个ensembl_gene_id的变量
> ensembl_gene_id<-rownames(raw_count_filt)
# 将行名单独变成一列加到raw_count_filt文件中
> raw_count_filt <- cbind(ensembl_gene_id,raw_count_filt)
# 添加列名
> colnames(raw_count_filt)[1] <- c("ensembl_gene_id")
# 将raw_count_filt文件和mms_symbols文件通过ensembl_gene_id而merge在一起，生成新的文件
> diff_name<-merge(raw_count_filt,mms_symbols,by="ensembl_gene_id")
# 查看新文件前6行
> head(diff_name)
代码执行结果：
执行结果
现在成功地生成了一个带有基因名称的count表达矩阵
参考链接：
https://blog.csdn.net/lztttao/article/details/82086346
https://www.jianshu.com/p/3a0e1e3e41d0
网友评论