首先声明，文章中的很多代码和方法参考了微信公众号：生信宝典。
在此表示感谢，书写此文章的目的是为了进行笔记以便后续研究备用，也希望能和大家一起分享学习，里面如有错误敬请指出
批次效应是测量结果中的一部分，它们因为实验条件的不同而具有不同的表现形式，并且与我们研究的变量没有关系。一般批次效应可能在下述情形中出现：
1.一个实验的不同部分在不同时间完成；
2.一个实验的不同部分由不同的人完成；
3.试剂用量不同、芯片不同、实验仪器不同；
4.将自己测的数据与从网上下载的数据混合使用；
其中，不同的时间，不同的人，以及试剂、仪器等等构成了批次效应。
R包sva可以预测并矫正批次效应，该包通过添加两个代理变量(surrogate variables)帮助DESeq在后续分析中矫正批次效应，并不会对数据本身进行修改。
在本研究数据中，batch列是在不同时间完成的生物学重复，而condition列则相当与处理条件。
利用上一篇建立的dds数据集，进行sva分析，代码来源于参考文献(RNA-seq workflow: gene-level exploratory analysis and differential expression)，其中比较关键的是mod的建立，里面选择condition进行建模，根据参考文献的意思就是假设我们不知道有不同批次的处理，而仅以condition设计的4个处理进行建模，并设置mod0为对照，具体的算法说明可以查阅文献(svaseq: removing batch effects and other unwanted noise from sequencing data)：
dds <- estimateSizeFactors(dds)
dat <- counts(dds, normalized = TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ condition, colData(dds))
mod0 <- model.matrix(~ 1, colData(dds))
svseq <- svaseq(dat, mod, mod0, n.sv = 2)
但实际上我们是知道存在不同的处理批次batch，通过可视化batch参数来查看差别：
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
stripchart(svseq$sv[, i] ~ dds$batch, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
}
看上去在两个代理参数中r1的批次差别并不明显，r2仅仅是略有差别。但是每张图上的四个点代表不同的处理，实际的分布情况究竟如何呢。也可以通过直接查看sv1和sv2的数值来判断：
> svseq$sv
[,1] [,2]
[1,] -0.19582386 0.2958170
[2,] -0.24163478 0.5719400
[3,] -0.34376511 -0.3164447
[4,] 0.53655461 -0.1749394
[5,] 0.33288941 -0.2062186
[6,] -0.50402035 -0.4996413
[7,] 0.05640151 0.3973242
[8,] 0.35939858 -0.0678371
每一行代表一个样品，列的数值是SV1和SV2的结果，我们知道行1和2对应R0_1和R0_2，行3和4对应R16_1和R16_2，以此类推。那么这两两之间分别SV1和SV2对应的数值差异越小越好。为了更方便的展示差异，将SV1和SV2参数传入dds，并进行可视化：
#将SV1和SV2参数传入dds，并重新命名，便于后续比较：
ddssva <- dds
ddssva$SV1 <- svseq$sv[,1]
ddssva$SV2 <- svseq$sv[,2]
design(ddssva) <- ~ SV1 + SV2 + condition
#使用ggplot进行可视化：
plot_data <- as.data.frame(colData(ddssva))
plot_data$Sample <- rownames(plot_data)
library(ggplot2)
ggplot(plot_data, aes(x=SV1, y=SV2, color=condition, shape=batch)) + geom_point()+ geom_text_repel(aes(label=Sample))
这下就可以看出明显差别了，按照正常情况，应该是R0，R16，R24和R32的两个批次应该分别聚集在一起，但是从图中可以看出R16和R24的平行性并不好。本研究的四个处理是细菌在同一培养条件下，不同生长时期获取的样品，进行了两次单独的生物学重复，而R16和R24正是细菌生长的对数期，这可能是造成两批重复在这两个样品中存在较大差异的原因。也提醒我们在后续进行RNA-seq实验设计是，生物学重复实验还是尽量需要在同一批培养中实施。
当将SV1和SV2修正参数传入ddssva后，后续我们会分别分析没有修正参数dds和带有修正参数的ddssva两个数据集在最终获取表达差异基因方面的差别。
后面我会分别对这两组数据集进行分析并比较他们之间的差别。
