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

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

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

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

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)
# 預(yù)處理，過濾低豐度的數(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中的對象
dds <- DESeqDataSetFromMatrix(
   countData = countData,
   colData = colData,
   design = ~ condition)
# 指定哪一組作為control
dds$condition <- relevel(dds$condition, ref = "control")

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

2. 歸一化

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

dds <- estimateSizeFactors(dds)

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

3. 估計基因的離散程度

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

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

通過下列代碼估算每個基因的α值

dds <- estimateDispersions(dds)

4. 差異分析

代碼如下

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

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

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

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

> 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

提取歸一化之后的表達量，自己計算baseMean和logFoldChange, 示例數(shù)據(jù)包含了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和自己計算的結(jié)果，可以發(fā)現(xiàn)，baseMeans是一致的，而log2Foldchange 差異很大，甚至連數(shù)值的正負都發(fā)生了變化。

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

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

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

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

廣告

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