189 8069 5689
这篇文章将为大家详细讲解有关DESeq2中怎么实现两组间差异分析操作，文章内容质量较高，因此小编分享给大家做个参考，希望大家阅读完这篇文章后对相关知识有一定的了解。
创新互联公司是一家集成都网站建设、网站设计、网站页面设计、网站优化SEO优化为一体的专业网站设计公司，已为成都等多地近百家企业提供网站建设服务。追求良好的浏览体验,以探求精品塑造与理念升华,设计最适合用户的网站页面。 合作只是第一步，服务才是根本,我们始终坚持讲诚信，负责任的原则，为您进行细心、贴心、认真的服务,与众多客户在蓬勃发展的市场环境中,互促共生。
读取基因的表达量表格和样本的分组信息两个文件，其中表达量的文件示例如下
gene_id ctrl-1 ctrl-2 ctrl-3 case-1 case-2 case-3 geneA 14 0 11 4 0 12 geneB 125 401 442 175 59 200
每一行为一个基因，每一列代表一个样本。
分组信息的文件示例如下
sample condition ctrl-1 control ctrl-2 control ctrl-3 control case-1 case case-2 case case-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) # 指定哪一组作为control dds$condition <- relevel(dds$condition, ref = "control")
在读取数据的过程中，有两点需要注意，第一个就是根据表达量对基因进行过滤，通常是过滤低表达量的基因，这一步是可选的，阈值可以自己定义；另外一个就是指定哪一组作为control组，在计算log2FD时 ，需要明确control组，默认会字符串顺序对分组的名字进行排序，排在前面的作为control组，这种默认行为选出的control可能与我们的实验设计不同，所以必须明确指定control组。
计算每个样本的归一化系数，代码如下
dds <- estimateSizeFactors(dds)
将原始的表达量除以每个样本的归一化系数，就得到了归一化之后的表达量。
DESeq2假定基因的表达量符合负二项分布，有两个关键参数，总体均值和离散程度α值， 如下图所示
这个α值衡量的是均值和方差之间的关系，表达式如下
通过下列代码估算每个基因的α值
dds <- estimateDispersions(dds)
代码如下
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 A DataFrame with 6 rows and 2 columns baseMean log2FoldChangegene1 7.471250 -0.8961954 gene2 18.145279 0.4222174 gene3 2.329461 -2.3216915 gene4 165.634256 -0.1974001 gene5 38.300621 1.3573162 gene6 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 log2FoldChange
gene1 7.471250 0.5380191
gene2 18.145279 1.3404422
gene3 2.329461 0.1991979
gene4 165.634256 0.8719078
gene5 38.300621 2.5621035
gene6 7.904819 3.5365201对比DESeq2和自己计算的结果，可以发现，baseMeans是一致的，而log2Foldchange 差异很大，甚至连数值的正负都发生了变化。
log2FD 反映的是不同分组间表达量的差异，这个差异由两部分构成，一种是样本间本身的差异，比如生物学重复样本间基因的表达量就有一定程度的差异，另外一部分就是我们真正感兴趣的，由于分组不同或者实验条件不同造成的差异。
关于DESeq2中怎么实现两组间差异分析操作就分享到这里了，希望以上内容可以对大家有一定的帮助，可以学到更多知识。如果觉得文章不错，可以把它分享出去让更多的人看到。