先了解完整的分析流程和工具：http://blog.sina.com.cn/s/blog_9cf2d3640102x9kx.html
这是一个系列文，包括：从标准的workflow开始，到更高级的数据操作和workflow个性化，最后是DESeq2的统计学原理以及常见的question解答
本文介绍在差异表达分析之前的操作步骤，主要是DESeq2导入数据和预处理，DESeq2对导入不同数据类型兼容性很好。
导入数据
Why un-normalized counts? DESeq2要求输入的数据必须是没有标准化的raw count，第i行第j列代表着在样本j分配到featrue(i)
的reads数或fragemnts数（feature可以是基因、转录本、ChIP-seq等NGS测的对应区段bin、定量质谱中的某肽段序列）至于怎么获取count data参考：http://master.bioconductor.org/packages/release/workflows/html/rnaseqGene.html
因为DESeq2模型内部本身会对文库大小做校正，所以不要用已经对文库做过scaled的数据
关于DESeqDataSet对象
这个对象实际上是继承于RangedSummarizedExperiment对象的（来自SummarizedExperiment package）也就是说从assay里获取的每一行(counts)都可以和genomic ranges关系上(比如gene exon/transcript)
这个关联对于下游进一步个性化分析非常有用，比如可以通过DESeq2获取了差异表达基因，对比ChIP-seq数据寻找在基因组坐标上距离该差异基因最近的ChIP-seq peaks
modeling中重要的一步就是design formula,添加适当的变量（DESeq2基于负二项分布线性模型）默认是把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)
再用定量数据导入工具tximport导入，直接可以为DESeq2所用
前三者是alignment-free或只是“轻度比对”，具有可以校正在不同样本中基因长度不同（不同isoform)速度快、省内存、省bam文件存储、灵敏度高的优点，参考https://www.bioinfo-scrounger.com/archives/411/
和 https://www.jianshu.com/p/6d4cba26bb60
Salmon软件使用说明：https://combine-lab.github.io/salmon/getting_started/
示例数据演示( gene-level DESeqDataSet object from Salmon quant.sf files)，分析自己的数据时dir更改为/path/to/quant/files，这里tximportData包只是为了提供演示数据而已，做自己的数据分析时不需要
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数据
files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
#把转录本和gene关联起来的table
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))
#用tximport function导入转录本定量数据
#关于tximport函数的详细说明和怎么构建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函数, 需要提供的输入对象包括表达矩阵(count-matrix)、样本信息samples (the columns of the count matrix) ，样本信息制作成dataframe，然后就是design formula.
示例数据来自pasilla包http://bioconductor.org/packages/pasilla
library("pasilla")#示例数据来自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导入数据，同样这里用的演示数据：
directory <- system.file("extdata", package="pasilla",
mustWork=TRUE)
#用list.files选择要读入的数据，并规定只读取包含"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
示例数据来自airway，这是一个非常经典的RNA-seq测试数据，导入后实际上是一个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
数据预处理
尽管在使用DESeq2函数前过滤低count的gene并不是必须的，但预过滤数据的好处是，去除那些只有很少reads的行以后，可以减少dds的存储，增加程序运行速度。这里演示一下简单的过滤，只保留那些至少有10条reads的基因。至于更严格的过滤方法（ independent filtering ），以后我会在推文里写到。
keep <- rowSums(counts(dds)) >= 10#对行求和函数rowSums()
dds <- dds[keep,]
关于因子水平
和limma里面的contrast 类似的是，DESeq2同样最好说明清楚哪个因子level对应的control，避免后面解释数据遇到麻烦（你不知道到底logFC是以谁做参照的，我们希望是以control作为参照，这样logFC>1就表示实验组表达大于对照组；但是R不知道，R默认按字母顺序排列因子顺序，所以我们要对这个因子顺序set一下：
#factor levels，写在前面的level作为参照
dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
#或者用relevel函数
dds$condition <- relevel(dds$condition, ref = "untreated")
#如果对应的level没有样本，可以这样把对应level删除
dds$condition <- droplevels(dds$condition)
合并技术重复的数据
DESeq2提供了一个collapseReplicates函数来把技术重复的样本数据合并到表达矩阵的一个列中；技术重复就是对同一样本测多次（multiple sequencing runs of the same library）。注意：不能对生物学重复数据用这种方法合并
major reference
https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
