欢迎关注”生信修炼手册”!
DESeq2 接受raw count的定量表格,然后根据样本分组进行差异分析,具体步骤如下
1. 读取数据读取基因的表达量表格和样本的分组信息两个文件,其中表达量的文件示例如下
gene_id ctrl-1 ctrl-2 ctrl-3 case-1 case-2 case-3geneA 14 0 11 4 0 12geneB 125 401 442 175 59 200每一行为一个基因,每一列代表一个样本。
分组信息的文件示例如下
sample conditionctrl-1 controlctrl-2 controlctrl-3 controlcase-1 casecase-2 casecase-3 case第一列为样本名,第二列为样本的分组信息。
读取文件的代码如下
# 读取表达量的表格count <- read.table( "gene.counts.tsv", header=T, sep="t", row.names=1, comment.char="", check.names=F)# 预处理,过滤低丰度的数据countData <- count[apply(count, 1, sum) > 0 , ]# 读取样本分组信息colData <- read.table( "sample.group.tsv", header=T, sep="t", row.names=1, comment.char="", check.names=F)# 构建DESeq2中的对象dds <- DESeqDataSetFromMatrix( countData = countData, colData = colData, design = ~ condition)# 指定哪一组作为controldds$condition <- relevel(dds$condition, ref = "control")在读取数据的过程中,有两点需要注意,第一个就是根据表达量对基因进行过滤,通常是过滤低表达量的基因,这一步是可选的,阈值可以自己定义;另外一个就是指定哪一组作为control组,在计算log2FD时 ,需要明确control组,默认会字符串顺序对分组的名字进行排序,排在前面的作为control组,这种默认行为选出的control可能与我们的实验设计不同,所以必须明确指定control组。
2. 归一化计算每个样本的归一化系数,代码如下
dds <- estimateSizeFactors(dds)
将原始的表达量除以每个样本的归一化系数,就得到了归一化之后的表达量。
DESeq2假定基因的表达量符合负二项分布,有两个关键参数,总体均值和离散程度α值, 如下图所示
这个α值衡量的是均值和方差之间的关系,表达式如下
通过下列代码估算每个基因的α值
dds <- estimateDispersions(dds) 4. 差异分析代码如下
dds <- nbinomWaldTest(dds)res <- results(dds)为了简化调用,将第二部到第四部封装到了DESeq这个函数中,代码如下
dds <- DESeq(dds)res <- results(dds)write.table( res, "DESeq2.diff.tsv", sep="t", quote=F, col.names = NA)在DESeq2差异分析的过程中,已经考虑到了样本之间已有的差异,所以可以发现,最终结果里的log2FD值和我们拿归一化之后的表达量计算出来的不同, 示意如下
> head(results(dds)[, 1:2])log2 fold change (MLE): condition B vs ADataFrame with 6 rows and 2 columns baseMean log2FoldChange <numeric> <numeric>gene1 7.471250 -0.8961954gene2 18.145279 0.4222174gene3 2.329461 -2.3216915gene4 165.634256 -0.1974001gene5 38.300621 1.3573162gene6 7.904819 1.8384322提取归一化之后的表达量,自己计算baseMean和logFoldChange, 示例数据包含了12个样本,其中前6个样本为control, 后6个样本为case , 代码如下
> nor_count <- counts(dds, nor = T)> res <- data.frame( baseMean = apply(nor_count, 1, mean), log2FoldChange = apply(nor_count, 1, function(t){ mean(t[7:12]) / mean(t[1:6])}) )> head(res) baseMean log2FoldChangegene1 7.471250 0.5380191gene2 18.145279 1.3404422gene3 2.329461 0.1991979gene4 165.634256 0.8719078gene5 38.300621 2.5621035gene6 7.904819 3.5365201对比DESeq2和自己计算的结果,可以发现,baseMeans是一致的,而log2Foldchange 差异很大,甚至连数值的正负都发生了变化。
log2FD 反映的是不同分组间表达量的差异,这个差异由两部分构成,一种是样本间本身的差异,比如生物学重复样本间基因的表达量就有一定程度的差异,另外一部分就是我们真正感兴趣的,由于分组不同或者实验条件不同造成的差异。
用归一化之后的数值直接计算出的log2FD包含了以上两种差异,而我们真正感兴趣的只有分组不同造成的差异,DESeq2在差异分析的过程中已经考虑到了样本本身的差异,其最终提供的log2FD只包含了分组间的差异,所以会与手动计算的不同。
·end·
—如果喜欢,快分享给你的朋友们吧—
扫描关注微信号,更多精彩内容等着你!
这个α值衡量的是均值和方差之间的关系,表达式如下
通过下列代码估算每个基因的α值
dds <- estimateDispersions(dds) 4. 差异分析代码如下
dds <- nbinomWaldTest(dds)res <- results(dds)为了简化调用,将第二部到第四部封装到了DESeq这个函数中,代码如下
dds <- DESeq(dds)res <- results(dds)write.table( res, "DESeq2.diff.tsv", sep="t", quote=F, col.names = NA)在DESeq2差异分析的过程中,已经考虑到了样本之间已有的差异,所以可以发现,最终结果里的log2FD值和我们拿归一化之后的表达量计算出来的不同, 示意如下
> head(results(dds)[, 1:2])log2 fold change (MLE): condition B vs ADataFrame with 6 rows and 2 columns baseMean log2FoldChange <numeric> <numeric>gene1 7.471250 -0.8961954gene2 18.145279 0.4222174gene3 2.329461 -2.3216915gene4 165.634256 -0.1974001gene5 38.300621 1.3573162gene6 7.904819 1.8384322提取归一化之后的表达量,自己计算baseMean和logFoldChange, 示例数据包含了12个样本,其中前6个样本为control, 后6个样本为case , 代码如下
> nor_count <- counts(dds, nor = T)> res <- data.frame( baseMean = apply(nor_count, 1, mean), log2FoldChange = apply(nor_count, 1, function(t){ mean(t[7:12]) / mean(t[1:6])}) )> head(res) baseMean log2FoldChangegene1 7.471250 0.5380191gene2 18.145279 1.3404422gene3 2.329461 0.1991979gene4 165.634256 0.8719078gene5 38.300621 2.5621035gene6 7.904819 3.5365201对比DESeq2和自己计算的结果,可以发现,baseMeans是一致的,而log2Foldchange 差异很大,甚至连数值的正负都发生了变化。
log2FD 反映的是不同分组间表达量的差异,这个差异由两部分构成,一种是样本间本身的差异,比如生物学重复样本间基因的表达量就有一定程度的差异,另外一部分就是我们真正感兴趣的,由于分组不同或者实验条件不同造成的差异。
用归一化之后的数值直接计算出的log2FD包含了以上两种差异,而我们真正感兴趣的只有分组不同造成的差异,DESeq2在差异分析的过程中已经考虑到了样本本身的差异,其最终提供的log2FD只包含了分组间的差异,所以会与手动计算的不同。
·end·
—如果喜欢,快分享给你的朋友们吧—
扫描关注微信号,更多精彩内容等着你!