lncRNA组装流程的软件介绍软件推荐之DEseq2
咱们《生信技能树》的B站有一个lncRNA数据分析实战，缺乏配套笔记，所以我们安排了100个lncRNA组装案例文献分享，以及这个流程会用到的100个软件的实战笔记教程！
做转录组RNA-seq的一个重要目的就是找到差异基因，而DEseq2就是一个用于差异分析的R包
官网使用说明：http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
一、软件原理
1.RNA-Seq中的统计检验问题
· 估算λg和Φg，的过程叫做estimate dispersion；
· Estimate dispersion不同软件的处理过程及策略不同；
· 通过Estimate dispersion确定λg和Φg，就可以进⾏统计检验。
2. DESeq2 Estimate dispersion的策略
第1步，通过极⼤似然估计粗略估计出各基因的dispersion参数；
若包含非常多的sample或者repeat数目，这⼀步基本上就能得到最终结果。
第2步，对极⼤似然估计的结果进⾏拟合，得到趋势线；
第3步，对于⼀些远离趋势线的点，向趋势线附近调整。
二、软件安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
三、输入数据
1.原始的整数型矩阵,就是每个样本中比对到每个基因的reads数
2.DESeqDataSet
DESeqDataSet是DESeq2包中储存read counts以及统计分析过程中的数据的一个“对象”，在代码中常表示为“dds”或deseq2.obj。
一个DESeqDataSet对象必须关联相应的design公式。design公式指明了要对哪些变量进行统计分析。该公式（上文中的design = ~batch + condition）以短波浪字符开头，中间用加号连接变量。design公式可以修改，但是相应的差异表达分析就需要重新做，因为design公式是用来估计统计模型的分散度以及log2 fold change的。
四、软件运行命令
以featureCouts结果为例：
rm(list=ls())
# make count table
setwd("~/Desktop/202006-RNASeq_course/")
raw_df <- read.table(file = "./03.code_and_data/out_table/293T-RNASeq-Ctrl_vs_KD.STAR.hg38.featureCounts.tsv",header = T,skip = 1,sep = "\t")
count_df <- raw_df[,c(7:10)]
rownames(count_df) <- as.character(raw_df[,1])
colnames(count_df) <- c("Ctrl_rep1","Ctrl_rep2","METTL3_KD_rep1","METTL3_KD_rep2")
# write.table(count_df, file = "./03.code_and_data/out_table/293T-RNASeq-Ctrl_vs_KD.STAR.hg38.featureCounts.FixColName.tsv",col.names = T,row.names = T,quote = F,sep = "\t")
###############################################################################
# Part I DESeq2
###############################################################################
rm(list=ls())
library(DESeq2)
# -------------------------------------------------------->>>>>>>>>>
# make obj
# -------------------------------------------------------->>>>>>>>>>
# count table
count_df <- read.table(file = "./03.code_and_data/out_table/293T-RNASeq-Ctrl_vs_KD.STAR.hg38.featureCounts.FixColName.tsv",header = T,sep = "\t")
# filter
colnames(count_df)
count_df.filter <- count_df[rowSums(count_df) > 20 & apply(count_df,1,function(x){ all(x > 0) }),]
# condition table
sample_df <- data.frame(
condition = c(rep("ctrl",2), rep("KD",2)),
cell_line = "293T"
)
rownames(sample_df) <- colnames(count_df.filter)
deseq2.obj <- DESeqDataSetFromMatrix(countData = count_df.filter, colData = sample_df, design = ~condition) # design = ~design = ~condition,差异的变动主要存在于那些地方：condition
deseq2.obj
# -------------------------------------------------------->>>>>>>>>>
# directly get test result
# -------------------------------------------------------->>>>>>>>>>
deseq2.obj <- DESeqDataSetFromMatrix(countData = count_df.filter, colData = sample_df, design = ~condition)
deseq2.obj
# test
deseq2.obj <- DESeq(deseq2.obj)
# get result
deseq2.obj.res <- results(deseq2.obj)
deseq2.obj.res.df <- as.data.frame(deseq2.obj.res)
# -------------------------------------------------------->>>>>>>>>>
# step by step get test result
# -------------------------------------------------------->>>>>>>>>>
deseq2.obj <- DESeqDataSetFromMatrix(countData = count_df, colData = sample_df, design = ~condition)
deseq2.obj
# normalization
deseq2.obj <- estimateSizeFactors(deseq2.obj)
sizeFactors(deseq2.obj)
# dispersion
deseq2.obj <- estimateDispersions(deseq2.obj)
dispersions(deseq2.obj)
# plot dispersion
plotDispEsts(deseq2.obj, ymin = 1e-4)
# test
deseq2.obj <- nbinomWaldTest(deseq2.obj)
deseq2.obj.res <- results(deseq2.obj)
# -------------------------------------------------------->>>>>>>>>>
# MA plot
# -------------------------------------------------------->>>>>>>>>>
DESeq2::plotMA(deseq2.obj.res, alpha=0.001)
五、关于因子的说明
默认情况下，R会安装字母顺序选择一个level作为参考level。就是说，如果没有给DEseq2函数指定要与哪个level进行比较（即指定哪个level是对照组），默认会安装字母顺序选择排在最前面的level作为对照组。
如果要指定哪个level是对照组，可以直接在函数#results#的参数contrast中指明对照level（下文中会具体介绍），或者直接在factor level中指明。
注意，更改design公式中变量的factor level只能在运行DESeq2分析前，分析中或分析后都不能再进行更改。更改dds矩阵中的factor level可以用以下两个方法：
#用factor
dds$condition <- factor(dds$condition, levels = c("untreated", "treated"))
#或者用relevel，直接指定参考level
dds$condition <- relevel(dds$condition, ref = "untreated")
如果需要从dds中提取一个子集，比如说从dds中去掉部分样品，这时有可能会同时去掉所有样品的一个或者多个level。这种情况下，可以用droplevles函数去掉没有对应sample的level：
dds$condition <- droplevels(dds$condition)
文末友情推荐
与十万人一起学生信，你值得拥有下面的学习班：