RNA-seq-123

1
2
3
4

library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)

1
2
3
4
5
6
7

#使用DGEList读取counts矩阵
genelist <- DGEList(counts = salmon_counts[2:31], group = group)
#使用tximport读入数据
## 读入分组信息

##读入基因注释信息
genelist$genes

1
2
3
4
5
6
7

#直接根据所有样品中的表达量进行过滤
table(rowSums(genelist$counts==0)==30)

#使用
keep.exprs <- filterByExpr(genelist, group = group)
genelist <- genelist[keep.exprs,, keep.lib.sizes=FALSE]
dim(genelist)

1
2
3

x <- calcNormFactors(genelist, method = "TMM")

x$genes <- x$genes[rownames(x$counts),]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# 输入设计矩阵
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
logcpm <- cpm(x, log = TRUE, prior.count = 2)

contr.matrix <- makeContrasts(
  IOZ07.vs.1621 = Mature - Mature1621,
  Mature.vs.vitro = Mature - Invitro,
  vitro.vs.1621 = Proliferation - Mature,
  Mature.vs.Transition = Mature - Transition,
  Transition.vs.Mycelium = Transition - Mycelium,
  levels = colnames(design)
)
v <- voom(x, design, plot = TRUE)
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
summary(decideTests(efit))