ballgown包可以读入Stringtie 的转录组组装及定量数据，进行基因差异表达分析。

1. 数据读入

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("ballgown")

require(ballgown)
library(dplyr)
library(genefilter)

project_dir <- '/home/test/rna_seq_data' # 项目文件夹，每个样品的数据放在不同的子文件夹里。

# 含有.ctab结尾的五个文件
# 注意samplePattern，文件名的前缀
m_bg <- ballgown(dataDir = project_dir,samplePattern = "SRR")

# length(structure(m_bg)$exon)
# length(structure(m_bg)$intron)
# length(structure(m_bg)$trans)

# 样品分组
pheno_data <- data.frame(id=sampleNames(m_bg), group=c(rep("1",3),rep("2",2)))
pData(m_bg) = pheno_data

2. 数据筛选

texpr(m_bg)[1:3,1:5] #  转录本表达矩阵

# 筛选FPKM在个样本中方差大于1的转录本
# subset {ballgown} 需要加载genefilter包
m_bg_filt = subset(m_bg,"rowVars(texpr(m_bg)) >1",genomesubset=TRUE)

3. 分析差异表达基因

# Identify genes that show statistically significant differences between groups
results_transcripts = stattest(m_bg_filt,feature="transcript", 
                               covariate="group", 
                               getFC=TRUE, meas="FPKM")

results_genes = stattest(m_bg_filt, feature="gene",covariate="group",
                        getFC=TRUE,meas="FPKM")

colnames(results_genes)

# Add gene names and gene IDs
results_transcripts =
  data.frame(geneNames=ballgown::geneNames(m_bg_filt),
             geneIDs=ballgown::geneIDs(m_bg_filt), results_transcripts)

head(results_transcripts)
colnames(results_transcripts)

# sort
results_transcripts = arrange(results_transcripts,pval) # arrange {dplyr}
results_genes = arrange(results_genes,pval)

head(results_transcripts)
head(results_genes)

# write to files
write.csv(results_transcripts, "/home/test/rna_seq_data/diff_transcript_results.csv",
          row.names=FALSE)
write.csv(results_genes, "/home/test/rna_seq_data/diff_gene_results.csv",
          row.names=FALSE)

4. 筛选差异表达基因

# choose
diff_trans_filt <- subset(results_transcripts,results_transcripts$qval<0.05)
diff_gene_filt <- subset(results_genes,results_genes$qval<0.05)

diff_trans_filt[1:3,1:5]

dim(diff_trans_filt)
dim(diff_gene_filt)

# p-value和q-value是统计学检验变量，衡量“假阳性概率”，应用到基因检测结果中，可衡量“某个基因差异 # 表达的假阳性概率”，代表差异显著性，小于0.05代表结果有差异。
# 如果p-value或q-value/越低，那么“该基因差异结果”是假阳性的概率就越低，可靠性就越高。
# q-value相比于p-value更加严格，当差异基因结果较少时，可退而求其次根据p-value筛选。
# 当然，用q值筛选可能会过滤掉少部分真的有差异的基因，所以，q值是个双刃剑。但，相比绝大部分基因的假 # 阳性，以及真阳性被滤掉的小概率，这部分的真阳性的丢失也不是很重要了。

5. 作图查看

tropical= c('darkorange', 'dodgerblue','hotpink',
            'limegreen', 'yellow')
palette(tropical)

fpkm = texpr(m_bg_filt,meas="FPKM")
fpkm = log2(fpkm+1)

# 不同组之间的fpkm差异
boxplot(fpkm,col=as.numeric(pData(m_bg)$group),
        las=2,ylab='log2(FPKM+1)')
# 转录本作图
plotTranscripts(ballgown::geneIDs(m_bg_filt)[200], m_bg_filt,
                main=c('Gene PARK7 in sample SRR12663677'),            sample=c('SRR12663677'))

plotTranscripts("merge.3987", m_bg_filt,
                main='merge.3987 transcripts', 
                sample=c('SRR12663677','SRR12663681'))

plotTranscripts(ballgown::geneIDs(m_bg_filt)[200], m_bg_filt,
                main= paste(as.character(ballgown::geneIDs(m_bg_filt)[200]),    'transcripts'), 
                sample=c('SRR12663677','SRR12663681'))


sampleNames(m_bg_filt)

plotMeans("merge.24191",
          m_bg_filt,groupvar="group",legend=FALSE)

plotMeans(ballgown::geneIDs(m_bg_filt)[200],
          m_bg_filt,groupvar="group",legend=FALSE)

edgeR基因表达差异分析

参考：
一个比较详细的例子：http://www.iwhgao.com/edger%E7%AE%80%E5%8D%95%E4%BD%BF%E7%94%A8%E6%89%8B%E5%86%8C/
bioconductor主页：http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
RNAseq123:http://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html

官方文档总结

注意⚠️：

• edgeR设计用于实际读取计数。我们不建议将预测的转录本丰度输入到edgeR而不是实际计数中。

读取read数

可以使用read.table()，read.delim(),如果文件比较多，readDGE(files, columns=)来一次性读取多个文件。注意，readDGE()需要指定两列，一列用于计数，一列用于基因标识符

x <- readDGE(files, columns=c(1,3))

DGEList对象、构建分组

DGEList是一个可以包含多种内容和统计的列表。DGEList至少需要的元素：counts、samples（包含group分组信息和lib.size文库大小），counts用来存放表达矩阵，samples用来标记样本信息和库的大小，group声明组别 。

# 构建一个含有组别标记的DGEList
y<- DGEList(counts=x)
group <-  c(1,1,2,2)
y<- DGEList(counts=x,group=group)

DGEList中分组是必要的，也可以添加其他信息进去，例如lane道，基因长度，基因注释等信息
一个注释过的DGEList

过滤，删除低表达基因

低表达基因不仅用不到，而且会干扰结果，所以要去除在任何样本中都没有足够多的序列片段的基因应该从下游分析中过滤掉，因为：

1. 低表达没有生物学意义
2. 去除低表达数据可以对数据中均值-方差关系有更精确的估计
3. 减少了观察差异表达下游分析中的运算量

edgeR包中的filterByExpr函数提供了自动过滤基因的方法，可保留尽可能多的有足够表达计数的基因。此函数默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有10个或更多序列片段计数的基因。
过滤标准是，以最小对组内样本数为标准，（此例最小组内样本为3），如果有基因在所有样本中表达数（count）小于10的个数超过最小组内样本数，就剔除该基因。换算为cpm即cut.off.cpm=10/<sum_col_counts>

CPM标度转换

常用的标度转换有CPM（counts per million）、log-CPM、FPKM、RPKM.
CPM是将counts转变为CPM指数counts per million，消除测序深度影响。
log-CPM是将CPMlog2化。cpm函数会在进行log2转换前给CPM值加上一个弥补值。默认的弥补值是2/L其中2是“预先计数”，而L是样本总序列数（以百万计）的平均值，所以log-CPM值是根据CPM值通过log2(CPM + 2/L)计算得到的。
对于一个基因，CPM值为1相当于在序列数量约2千万的样品中，有20个计数，或者在序列数量约7.6千万有76个计数。

cpm<- edgeR::cpm(x)
lcpm <- edgeR::cpm(x, log=TRUE)

手动过滤

经验设置为cpm=1位为cutoff点。但是，这并不是最精准。因为随着测序深度增加，例如20million（2 千万），cpm=1 意味着 counts=20。阈值可能会有点高。测序深度低的话，例如2million（2百万），cpm=1 意味着counts=1。阈值可能会太低。此时可以使用自动过滤，或者根据cut.off.cpm=10/<sum_col_counts>来计算。

cpm(y)
keep_cpm <- rowSums(cpm(y)>1) >= 3      #此例设置最小组内样本数为3
x<- x[keep_cpm,,keep.lib.sizes=FALSE]

自动过滤

使用edgeR::filterByExpr()自动过滤

keep.exprs<- edgeR::filterByExpr(x,group=group)
x<- x[keep.exprs,,keep.lib.sizes=FALSE]

归一化

归一化不是绝对必要的，但是推荐进行归一化。
有重复的样本中，应该不具备生物学意义的外部因素会影响单个样品的表达，例如中第一批制备的样品会总体上表达高于第二批制备的样品，假设所有样品表达值的范围和分布都应当相似，需要进行归一化来确保整个实验中每个样本的表达分布都相似。可以使用edgeR::calcNormFctors()函数来用TMM方法进行归一化,通过计算归一化系数来决定文库大小的缩放系数当进行归一化时，如果对象是DGElist，那么这些归一化系数会被自动存在x$samples$norm.factors

测序深度

不同文库大小代表不同测序深度。这是基本建模过程的一部分，并自动进入倍数变化或p值计算。它始终存在，不需要任何用户干预。

有效库大小

在某些情况下，如果一小部分高表达基因消耗了一个特定样本文库大小的很大一部分，这将导致该样本剩余基因的采样不足。除非对这种效应进行调整，否则在该样本中其余的基因可能会错误地表现为下调。（因为高表达基因比较多，消耗了大量的测序资源，导致非高表达基因会相对偏低。）
使用函数edgeR::calcNormFactors()，默认使用TMM方法进行归一化，归一化后，会给样品分配缩放系数。
将原始库大小与缩放因子的乘积称为有效库大小。有效的库大小会在所有下游分析中替换原始库大小。

y <- calcNormFactors(y)

注意⚠️：
归一化并不会直接在counts数值上修改，而是归一化系数会被自动存在x$samples$norm.factors。

GC含量

基因长度

因为GC含量和基因长度在每个基因间不会改变，因此这些都是相对的，因此对差异分析影响很小。

MDS图形展示 样本无监督聚类

这种图表使用无监督聚类方法展示出了样品间的相似性和不相似性，能让我们在进行正式的检验之前对于能检测到多少差异表达基因有个大致概念。

# 图形展示
library(RColorBrewer)
lcpm <- cpm(y, log=TRUE)
col.group <- group
levels(col.group) <-  brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
plotMDS(lcpm,labels=group,col=col.group)

#或者直接使用 
plotMDS(y)

负二项式模型

计算生物变异系数

edgeR中对于离散度的估计主要有两种模型，其中比较经典的是qCML。但qCML方法仅适用于单因素的数据集。例如，比较癌和正常样本间的基因表达差异，而不考虑其他因素（如年龄、性别等）的影响。qCML后续会使用与fisher精确检验比较类似的 exact test 进行差异表达分析。
如果使用qCML估计离散度，则不需要预先设定实验矩阵！(因为就只看一个因素，不必大费周折)
而对于更复杂的实验设计，则推荐使用GLMs模型，此时就需要提供实验矩阵design。
qCML方法：仅适用于具有单因素设计的数据集，但是有着更可靠的性能，尤其在小样本下有很好的表现，比较切合NGS测序。
通过estimateDisp()来估计公共离散度和 tagwise离散度（一个命令来运行）

y < -estimateDisp（y）

或者先计算common离散度，再计算tagwise离散度

y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

计算差异基因

et<- exactTest(y)
topTags(et)

广义线性模型（Glm）

对于更复杂的实验设计（有多个因素），可以通过广义线性模型来及性能你和。

计算离散度

通过下面来估计common离散度、trended离散度、tagwise离散度。

y<- estimateDisp(y,design)

或者分开依次进行

y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTrendedDisp(y, design)  #估计trended离散度
y <- estimateGLMTagwiseDisp(y, design)

关于design的构建(?)：

#group是区别分组的factor对应不同样本，lane是不同样本的的lane道
design <- model.matrix(~0+group+lane) 
colnames(design) <- gsub("group", "", colnames(design))

计算DE基因

edgrR涉及到差异表达分析的函数有很多： exactTest、glmFit、glmLRT、glmQLFit、glmQLFTest。
qCML估计离散度需要搭配 exact test 进行差异表达分析，对应 exactTest 函数。
而其他四个glm*都是与GLM模型搭配使用的函数。其中，glmFit 和 glmLRT 函数是配对使用的，用于 likelihood ratio test (似然比检验)，而 glmQLFit和 glmQLFTest则配对用于 quasi-likelihood F test (拟极大似然F检验)。
有两个方法glmQLFit()和glmQLFTest(),在两个检验方法中，首选QLFit，因为它反映了估计每个基因的离散度时的不确定性。当重复次数较少时，它提供了更强大和可靠的错误率控制

# 例子
group <- factor(c(1,1,2,2,3,3))
design <- model.matrix(~group)
fit <- glmQLFit(y, design)

如果没有重复样本？

首先确定一个良好的BCV，人类数据一般设定为0.4
一个例子：

bcv <- 0.2
counts <- matrix( rnbinom(40,size=1/bcv^2,mu=10), 20,2)
y <- DGEList(counts=counts, group=1:2)
et <- exactTest(y, dispersion=bcv^2)

如果是人类数据：

y_bcv <- y
bcv <- 0.4
et <- exactTest(y_bcv,dispersion = bcv ^2)

注：
如果差异基因过少，可以尝试下调BCV

输出结果

result = topTags(et, n = nrow(et$table))$table

查看统计

summary(de <- decideTestsDGE(et))

我想建立并管理一个高质量的生信&统计相关的微信讨论群，如果你想参与讨论，可以添加我微信： veryqun 。我会拉你进群，当然有问题也可以微信咨询我。

因为我的笔记分布CSDN、简书、知乎专栏等比较零散，管理起来比较麻烦，思考再三申请了一个 微信公众号，会更加方便地发布更多有关生信息、统计方面内容，如果你觉得有需要可以尝试性和我缔结契约 ^ . ^ 。公众号如下：

对我唯一的精神鼓励可能就是下方的点赞了吧 ^ ^