lncRNA组装流程的软件介绍软件推荐之DEseq2

咱们《生信技能树》的B站有一个lncRNA数据分析实战，缺乏配套笔记，所以我们安排了100个lncRNA组装案例文献分享，以及这个流程会用到的100个软件的实战笔记教程！

下面是100个lncRNA组装流程的软件的笔记教程

做转录组RNA-seq的一个重要目的就是找到差异基因，而DEseq2就是一个用于差异分析的R包

官网使用说明：http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

一、软件原理

1.RNA-Seq中的统计检验问题

代码语言：javascript复制

• 估算λg和Φg，的过程叫做estimate dispersion；

• Estimate dispersion不同软件的处理过程及策略不同；

• 通过Estimate dispersion确定λg和Φg，就可以进⾏统计检验。

2. DESeq2 Estimate dispersion的策略

第1步，通过极⼤似然估计粗略估计出各基因的dispersion参数；

若包含非常多的sample或者repeat数目，这⼀步基本上就能得到最终结果。

第2步，对极⼤似然估计的结果进⾏拟合，得到趋势线；

第3步，对于⼀些远离趋势线的点，向趋势线附近调整。

二、软件安装

代码语言：javascript复制

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

三、输入数据

代码语言：javascript复制

1.原始的整数型矩阵,就是每个样本中比对到每个基因的reads数

2.DESeqDataSet
DESeqDataSet是DESeq2包中储存read counts以及统计分析过程中的数据的一个“对象”，在代码中常表示为“dds”或deseq2.obj。

一个DESeqDataSet对象必须关联相应的design公式。design公式指明了要对哪些变量进行统计分析。该公式（上文中的design = ~batch   condition）以短波浪字符开头，中间用加号连接变量。design公式可以修改，但是相应的差异表达分析就需要重新做，因为design公式是用来估计统计模型的分散度以及log2 fold change的。

四、软件运行命令

以featureCouts结果为例：

代码语言：javascript复制

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可以用以下两个方法：

代码语言：javascript复制

#用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：

代码语言：javascript复制

dds$condition <- droplevels(dds$condition)

na repeat sample

0 人点赞