DESeq2中怎么實(shí)現(xiàn)兩組間差異分析操作

這篇文章將為大家詳細(xì)講解有關(guān)DESeq2中怎么實(shí)現(xiàn)兩組間差異分析操作，文章內(nèi)容質(zhì)量較高，因此小編分享給大家做個(gè)參考，希望大家閱讀完這篇文章后對(duì)相關(guān)知識(shí)有一定的了解。

創(chuàng)新互聯(lián)公司是一家集成都網(wǎng)站建設(shè)、網(wǎng)站設(shè)計(jì)、網(wǎng)站頁(yè)面設(shè)計(jì)、網(wǎng)站優(yōu)化SEO優(yōu)化為一體的專(zhuān)業(yè)網(wǎng)站設(shè)計(jì)公司，已為成都等多地近百家企業(yè)提供網(wǎng)站建設(shè)服務(wù)。追求良好的瀏覽體驗(yàn),以探求精品塑造與理念升華,設(shè)計(jì)最適合用戶(hù)的網(wǎng)站頁(yè)面。 合作只是第一步，服務(wù)才是根本,我們始終堅(jiān)持講誠(chéng)信，負(fù)責(zé)任的原則，為您進(jìn)行細(xì)心、貼心、認(rèn)真的服務(wù),與眾多客戶(hù)在蓬勃發(fā)展的市場(chǎng)環(huán)境中,互促共生。

1. 讀取數(shù)據(jù)

讀取基因的表達(dá)量表格和樣本的分組信息兩個(gè)文件，其中表達(dá)量的文件示例如下

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

每一行為一個(gè)基因，每一列代表一個(gè)樣本。

分組信息的文件示例如下

sample  condition
ctrl-1    control
ctrl-2    control
ctrl-3    control
case-1  case
case-2  case
case-3  case

第一列為樣本名，第二列為樣本的分組信息。

讀取文件的代碼如下

# 讀取表達(dá)量的表格
count <- read.table(
  "gene.counts.tsv",
  header=T,
  sep="\t",
  row.names=1,
  comment.char="",
  check.names=F)
# 預(yù)處理，過(guò)濾低豐度的數(shù)據(jù)
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)
# 構(gòu)建DESeq2中的對(duì)象
dds <- DESeqDataSetFromMatrix(
   countData = countData,
   colData = colData,
   design = ~ condition)
# 指定哪一組作為control
dds$condition <- relevel(dds$condition, ref = "control")

在讀取數(shù)據(jù)的過(guò)程中，有兩點(diǎn)需要注意，第一個(gè)就是根據(jù)表達(dá)量對(duì)基因進(jìn)行過(guò)濾，通常是過(guò)濾低表達(dá)量的基因，這一步是可選的，閾值可以自己定義；另外一個(gè)就是指定哪一組作為control組，在計(jì)算log2FD時(shí) ，需要明確control組，默認(rèn)會(huì)字符串順序?qū)Ψ纸M的名字進(jìn)行排序，排在前面的作為control組，這種默認(rèn)行為選出的control可能與我們的實(shí)驗(yàn)設(shè)計(jì)不同，所以必須明確指定control組。

2. 歸一化

計(jì)算每個(gè)樣本的歸一化系數(shù)，代碼如下

dds <- estimateSizeFactors(dds)

將原始的表達(dá)量除以每個(gè)樣本的歸一化系數(shù)，就得到了歸一化之后的表達(dá)量。

3. 估計(jì)基因的離散程度

DESeq2假定基因的表達(dá)量符合負(fù)二項(xiàng)分布，有兩個(gè)關(guān)鍵參數(shù)，總體均值和離散程度α值， 如下圖所示

這個(gè)α值衡量的是均值和方差之間的關(guān)系，表達(dá)式如下

通過(guò)下列代碼估算每個(gè)基因的α值

dds <- estimateDispersions(dds)

4. 差異分析

代碼如下

dds <- nbinomWaldTest(dds)
res <- results(dds)

為了簡(jiǎn)化調(diào)用，將第二部到第四部封裝到了DESeq這個(gè)函數(shù)中，代碼如下

dds <- DESeq(dds)
res <- results(dds)
write.table(
  res,
  "DESeq2.diff.tsv",
  sep="\t",
  quote=F,
  col.names = NA)

在DESeq2差異分析的過(guò)程中，已經(jīng)考慮到了樣本之間已有的差異，所以可以發(fā)現(xiàn)，最終結(jié)果里的log2FD值和我們拿歸一化之后的表達(dá)量計(jì)算出來(lái)的不同,  示意如下

> head(results(dds)[, 1:2])
log2 fold change (MLE): condition B vs A
DataFrame with 6 rows and 2 columns
        baseMean log2FoldChange
       <numeric>      <numeric>
gene1   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

提取歸一化之后的表達(dá)量，自己計(jì)算baseMean和logFoldChange, 示例數(shù)據(jù)包含了12個(gè)樣本，其中前6個(gè)樣本為control, 后6個(gè)樣本為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

對(duì)比DESeq2和自己計(jì)算的結(jié)果，可以發(fā)現(xiàn)，baseMeans是一致的，而log2Foldchange 差異很大，甚至連數(shù)值的正負(fù)都發(fā)生了變化。

log2FD 反映的是不同分組間表達(dá)量的差異，這個(gè)差異由兩部分構(gòu)成，一種是樣本間本身的差異，比如生物學(xué)重復(fù)樣本間基因的表達(dá)量就有一定程度的差異，另外一部分就是我們真正感興趣的，由于分組不同或者實(shí)驗(yàn)條件不同造成的差異。

關(guān)于DESeq2中怎么實(shí)現(xiàn)兩組間差異分析操作就分享到這里了，希望以上內(nèi)容可以對(duì)大家有一定的幫助，可以學(xué)到更多知識(shí)。如果覺(jué)得文章不錯(cuò)，可以把它分享出去讓更多的人看到。

新聞標(biāo)題：DESeq2中怎么實(shí)現(xiàn)兩組間差異分析操作
文章源于：http://muchs.cn/article46/phoceg.html

成都網(wǎng)站建設(shè)公司_創(chuàng)新互聯(lián)，為您提供網(wǎng)站改版、、品牌網(wǎng)站設(shè)計(jì)、網(wǎng)站建設(shè)、微信公眾號(hào)、面包屑導(dǎo)航

廣告

聲明：本網(wǎng)站發(fā)布的內(nèi)容（圖片、視頻和文字）以用戶(hù)投稿、用戶(hù)轉(zhuǎn)載內(nèi)容為主，如果涉及侵權(quán)請(qǐng)盡快告知，我們將會(huì)在第一時(shí)間刪除。文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng)，如需處理請(qǐng)聯(lián)系客服。電話：028-86922220；郵箱：631063699@qq.com。內(nèi)容未經(jīng)允許不得轉(zhuǎn)載，或轉(zhuǎn)載時(shí)需注明來(lái)源： 創(chuàng)新互聯(lián)