先了解完整的分析流程和工具：http://blog.sina.com.cn/s/blog_9cf2d3640102x9kx.html
這是一個系列文遍略，包括：從標準的workflow開始掖鱼，到更高級的數(shù)據(jù)操作和workflow個性化逐沙，最后是DESeq2的統(tǒng)計學原理以及常見的question解答
本文介紹在差異表達分析之前的操作步驟帆喇，主要是DESeq2導入數(shù)據(jù)和預處理，DESeq2對導入不同數(shù)據(jù)類型兼容性很好贱勃。
導入數(shù)據(jù)
Why un-normalized counts? DESeq2要求輸入的數(shù)據(jù)必須是沒有標準化的raw count纹因，第i行第j列代表著在樣本j分配到featrue(i)
的reads數(shù)或fragemnts數(shù)（feature可以是基因壕曼、轉錄本、ChIP-seq等NGS測的對應區(qū)段bin寸齐、定量質譜中的某肽段序列）至于怎么獲取count data參考：http://master.bioconductor.org/packages/release/workflows/html/rnaseqGene.html
因為DESeq2模型內(nèi)部本身會對文庫大小做校正欲诺，所以不要用已經(jīng)對文庫做過scaled的數(shù)據(jù)
關于DESeqDataSet對象
這個對象實際上是繼承于RangedSummarizedExperiment對象的（來自SummarizedExperiment package）也就是說從assay里獲取的每一行(counts)都可以和genomic ranges關系上(比如gene exon/transcript)
這個關聯(lián)對于下游進一步個性化分析非常有用，比如可以通過DESeq2獲取了差異表達基因渺鹦，對比ChIP-seq數(shù)據(jù)尋找在基因組坐標上距離該差異基因最近的ChIP-seq peaks
modeling中重要的一步就是design formula,添加適當?shù)淖兞浚―ESeq2基于負二項分布線性模型）默認是把control放在前面扰法，感興趣變量放在后面：design~control+variable
四種構建DESeqDataSet對象的方法
Transcript abundance files and tximport / tximeta
pipeline里推薦的一步是使用快速轉錄本定量工具獲取counts，比如
- Salmon (Patro et al. 2017)
- Sailfish (Patro, Mount, and Kingsford 2014)
- kallisto (Bray et al. 2016)
- RSEM (Li and Dewey 2011)
再用定量數(shù)據(jù)導入工具tximport導入毅厚，直接可以為DESeq2所用
前三者是alignment-free或只是“輕度比對”塞颁，具有可以校正在不同樣本中基因長度不同（不同isoform)速度快、省內(nèi)存吸耿、省bam文件存儲祠锣、靈敏度高的優(yōu)點，參考https://www.bioinfo-scrounger.com/archives/411/
和 http://www.reibang.com/p/6d4cba26bb60
Salmon軟件使用說明：https://combine-lab.github.io/salmon/getting_started/
示例數(shù)據(jù)演示( gene-level DESeqDataSet object from Salmon quant.sf files)珍语，分析自己的數(shù)據(jù)時dir更改為/path/to/quant/files锤岸，這里tximportData包只是為了提供演示數(shù)據(jù)而已，做自己的數(shù)據(jù)分析時不需要
library("tximport")
library("readr")
library("tximportData")
dir <- system.file("extdata", package="tximportData")
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
samples$condition <- factor(rep(c("A","B"),each=3))
# pop center assay sample experiment run condition
# 1 TSI UNIGE NA20503.1.M_111124_5 ERS185497 ERX163094 ERR188297 A
# 2 TSI UNIGE NA20504.1.M_111124_7 ERS185242 ERX162972 ERR188088 A
# 3 TSI UNIGE NA20505.1.M_111124_6 ERS185048 ERX163009 ERR188329 A
# 4 TSI UNIGE NA20507.1.M_111124_7 ERS185412 ERX163158 ERR188288 B
# 5 TSI UNIGE NA20508.1.M_111124_2 ERS185362 ERX163159 ERR188021 B
# 6 TSI UNIGE NA20514.1.M_111124_4 ERS185217 ERX163062 ERR188356 B
rownames(samples) <- samples$run
samples[,c("pop","center","run","condition")]
## pop center run condition
## ERR188297 TSI UNIGE ERR188297 A
## ERR188088 TSI UNIGE ERR188088 A
## ERR188329 TSI UNIGE ERR188329 A
## ERR188288 TSI UNIGE ERR188288 B
## ERR188021 TSI UNIGE ERR188021 B
## ERR188356 TSI UNIGE ERR188356 B
#讀取每個run數(shù)據(jù)
files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
#把轉錄本和gene關聯(lián)起來的table
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))
#用tximport function導入轉錄本定量數(shù)據(jù)
#關于tximport函數(shù)的詳細說明和怎么構建tx2gene把轉錄本和基因對應上板乙，參考：http://bioconductor.org/packages/release/bioc/html/tximport.html
txi <- tximport(files, type="salmon", tx2gene=tx2gene)
# reading in files with read_tsv
# 1 2 3 4 5 6
# summarizing abundance
# summarizing counts
# summarizing length
library("DESeq2")#依賴包非常豐富是偷，其中就包括了GenomicRanges拳氢、SummarizedExperiment等
ddsTxi <- DESeqDataSetFromTximport(txi,
colData = samples,
design = ~ condition)
#using counts and average transcript lengths from tximport
#tximeta這個包還可以自動添加annotation metadata(會自動下載gtf文件)，不需要tx2gene蛋铆，詳見https://bioconductor.org/packages/tximeta
coldata <- samples
coldata$files <- files
coldata$names <- coldata$run
library("tximeta")
se <- tximeta(coldata)
ddsTxi <- DESeqDataSet(se, design = ~ condition)
#ddsTxi object here can then be used as dds in the following analysis steps
Count matrix input
關于DESeqDataSetFromMatrix函數(shù), 需要提供的輸入對象包括表達矩陣(count-matrix)馋评、樣本信息samples (the columns of the count matrix) ，樣本信息制作成dataframe刺啦，然后就是design formula.
示例數(shù)據(jù)來自pasilla包http://bioconductor.org/packages/pasilla
library("pasilla")#示例數(shù)據(jù)來自pasilla包http://bioconductor.org/packages/pasilla
pasCts <- system.file("extdata",
"pasilla_gene_counts.tsv",
package="pasilla", mustWork=TRUE)
pasAnno <- system.file("extdata",
"pasilla_sample_annotation.csv",
package="pasilla", mustWork=TRUE)
cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))
coldata <- read.csv(pasAnno, row.names=1)
coldata <- coldata[,c("condition","type")]
#注意檢查colData和count matrix的rowname是否完全一致
head(cts,2)
## untreated1 untreated2 untreated3 untreated4 treated1 treated2
## FBgn0000003 0 0 0 0 0 0
## FBgn0000008 92 161 76 70 140 88
## treated3
## FBgn0000003 1
## FBgn0000008 70
coldata
## condition type
## treated1fb treated single-read
## treated2fb treated paired-end
## treated3fb treated paired-end
## untreated1fb untreated single-read
## untreated2fb untreated single-read
## untreated3fb untreated paired-end
## untreated4fb untreated paired-end
rownames(coldata) <- sub("fb", "", rownames(coldata))
all(rownames(coldata) %in% colnames(cts))#修改rownames
## [1] TRUE
all(rownames(coldata) == colnames(cts))#但order不一致
## [1] FALSE
cts <- cts[, rownames(coldata)]#重排一下
all(rownames(coldata) == colnames(cts))
## [1] TRUE
library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ condition)
dds
# class: DESeqDataSet
# dim: 14599 7
# metadata(1): version
# assays(1): counts
# rownames(14599): FBgn0000003 FBgn0000008 ... FBgn0261574 FBgn0261575
# rowData names(0):
# colnames(7): treated1 treated2 ... untreated3 untreated4
# colData names(2): condition type
#如果有額外的meatdata要添加留特，可以直接添加到DESeqDataSet，(這里只是為了演示）
featureData <- data.frame(gene=rownames(cts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
mcols(dds)
# DataFrame with 14599 rows and 1 column
# gene
# <factor>
# FBgn0000003 FBgn0000003
# FBgn0000008 FBgn0000008
# FBgn0000014 FBgn0000014
# FBgn0000015 FBgn0000015
# FBgn0000017 FBgn0000017
# ... ...
# FBgn0261571 FBgn0261571
# FBgn0261572 FBgn0261572
# FBgn0261573 FBgn0261573
# FBgn0261574 FBgn0261574
# FBgn0261575 FBgn0261575
htseq-count input
如果是用HTSeq做的轉錄本定量玛瘸，可以用DESeqDataSetFromHTSeqCount導入數(shù)據(jù)蜕青，同樣這里用的演示數(shù)據(jù)：
directory <- system.file("extdata", package="pasilla",
mustWork=TRUE)
#用list.files選擇要讀入的數(shù)據(jù)，并規(guī)定只讀取包含"treated"字符串的文件
sampleFiles <- grep("treated",list.files(directory),value=TRUE)
sampleCondition <- sub("(.*treated).*","\\1",sampleFiles)#condition用于制作colData
sampleTable <- data.frame(sampleName = sampleFiles,
fileName = sampleFiles,
condition = sampleCondition)#需要把每一個文件和condition對應起來
library("DESeq2")
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
directory = directory,
design= ~ condition)
ddsHTSeq
# class: DESeqDataSet
# dim: 70463 7
# metadata(1): version
# assays(1): counts
# rownames(70463): FBgn0000003:001 FBgn0000008:001 ... FBgn0261575:001
# FBgn0261575:002
# rowData names(0):
# colnames(7): treated1fb.txt treated2fb.txt ... untreated3fb.txt
# untreated4fb.txt
# colData names(1): condition
SummarizedExperiment input
示例數(shù)據(jù)來自airway糊渊，這是一個非常經(jīng)典的RNA-seq測試數(shù)據(jù)右核，導入后實際上是一個RangedSummarizedExperiment 對象，簡稱為"se"渺绒，如果不熟悉請參考：http://bioconductor.org/packages/release/data/experiment/html/airway.html
library("airway")
data("airway")
se <- airway
library("DESeq2")
ddsSE <- DESeqDataSet(se, design = ~ cell + dex)
ddsSE
# class: DESeqDataSet
# dim: 64102 8
# metadata(2): '' version
# assays(1): counts
# rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
# rowData names(0):
# colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
# colData names(9): SampleName cell ... Sample BioSample
數(shù)據(jù)預處理
盡管在使用DESeq2函數(shù)前過濾低count的gene并不是必須的贺喝，但預過濾數(shù)據(jù)的好處是，去除那些只有很少reads的行以后宗兼，可以減少dds的存儲躏鱼，增加程序運行速度。這里演示一下簡單的過濾殷绍，只保留那些至少有10條reads的基因染苛。至于更嚴格的過濾方法（ independent filtering ），以后我會在推文里寫到篡帕。
keep <- rowSums(counts(dds)) >= 10#對行求和函數(shù)rowSums()
dds <- dds[keep,]
關于因子水平
和limma里面的contrast 類似的是殖侵，DESeq2同樣最好說明清楚哪個因子level對應的control，避免后面解釋數(shù)據(jù)遇到麻煩（你不知道到底logFC是以誰做參照的镰烧，我們希望是以control作為參照拢军，這樣logFC>1就表示實驗組表達大于對照組；但是R不知道怔鳖，R默認按字母順序排列因子順序茉唉，所以我們要對這個因子順序set一下：
#factor levels，寫在前面的level作為參照
dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
#或者用relevel函數(shù)
dds$condition <- relevel(dds$condition, ref = "untreated")
#如果對應的level沒有樣本结执，可以這樣把對應level刪除
dds$condition <- droplevels(dds$condition)
合并技術重復的數(shù)據(jù)
DESeq2提供了一個collapseReplicates函數(shù)來把技術重復的樣本數(shù)據(jù)合并到表達矩陣的一個列中度陆；技術重復就是對同一樣本測多次（multiple sequencing runs of the same library）。注意：不能對生物學重復數(shù)據(jù)用這種方法合并
major reference
https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
