搜索
查看: 2413|回复: 7

[mRNA-seq] 转录组入门分析7_差异表达分析

[复制链接]

9

主题

15

帖子

126

积分

注册会员

Rank: 2

积分
126
发表于 2017-9-1 06:52:08 | 显示全部楼层 |阅读模式
本帖最后由 laofuzi 于 2017-9-1 06:52 编辑

[AppleScript] 纯文本查看 复制代码
library(DESeq2)
#########################
#Build the Count Matrix
#load data
dataFrame <- read.csv("/mnt/d/AKAP95/matrix/m_raw_count.csv", header = TRUE, sep =",")  
#view the head of data
head(dataFrame)
#convert the datatype of dataFrame from dataframe into matrix, and extract the column 2-5
countMatrix <- as.matrix(dataFrame[2:5])
#add the rowname for countMatrix using the first column of dataFrame
rownames(countMatrix) <- dataFrame$gene_id
#view the head of count matrix
head(countMatrix)
########################
#Build the condition factor
condition <- factor(c(rep("control",2), rep("akap95",2)), levels = c("control","akap95"))
########################
#Build the Sample Information Table. First method, sample information table, also named as colData, or design matrix
#sample names
sampleNames <- c("m_control1", "m_control2", "akap95_1", "akap95_2")
#Build the Object colData
colData <- data.frame(name= c("m_control1", "m_control2", "akap95_1", "akap95_2"), condition) 
#add the rowname for colData using the vector sampleNames
rownames(colData) <- sampleNames
#view the head of colData
head(colData)
##########
# Build the Sample Information table. Second method, sample information table, also named as colData, or design matrix
#colData1 <- data.frame(row.names=colnames(countMatrix), name= c("m_control1", "m_control2", "akap95_1", "akap95_2"), condition)
#head(colData1)
########################
#Build the Object dds (DESeqDataSet)
dds <- DESeqDataSetFromMatrix(countMatrix, colData, design= ~ condition)
# view the structure of dds
dds
########################
# filter the row which the count numbers of all samples are 0
dds <- dds[rowSums(counts(dds)) > 1,]
# view the structure of dds
dds
########################
# normalize data
dds2 <- DESeq(dds)
# view results name
resultsNames(dds2)
# acquire the results using function results(), and assign to res
res <- results(dds2)
# view the summary of results
summary(res)
########################
# DESeq2 uses the so-called Benjamini-Hochberg (BH) adjustment for multiple testing problem
# padj is BH-adjusted p values which are given in the column padj of the results object
res <- res[order(res$padj),]
resdata <- merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)
diff_gene <-subset(resdata, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
write.csv(resdata, file= "/mnt/d/AKAP95/matrix/all_gene.csv",row.names = F)
write.csv(diff_gene, file= "/mnt/d/AKAP95/matrix/diff_gene.csv",row.names = F)
########################
# Plot MA, M means log fold change, A means the count mean. 
library(geneplotter)
plotMA(res, main="Plot MA", ylim=c(-5,5))
# Dispersion plot
plotDispEsts( dds2, ylim = c(1e-6, 1e1) )
#Histogram of the p-value
hist( res$pvalue, breaks=20, col="grey" )
#Heatmap
sum(res$padj < 0.05, na.rm=TRUE)
library("pheatmap")
select <- order(rowMeans(counts(dds2, normalized=TRUE)), decreasing=TRUE)[1:1000]
nt <- normTransform(dds2) # defaults to log2(x+1)
log2.norm.counts <- assay(nt)[select,]
df <- as.data.frame(colData(dds2)[,c("name", "condition")])
pdf('/mnt/d/AKAP95/matrix/heatmap1000.pdf', width = 6, height = 7)
pheatmap(log2.norm.counts, cluster_rows=TRUE, show_rownames=FALSE, cluster_cols=TRUE, annotation_col=df)
dev.off()




上一篇:Trinity的安装与使用
下一篇:《生信菜鸟团》公众号文章合辑
回复

使用道具 举报

6

主题

26

帖子

363

积分

中级会员

Rank: 3Rank: 3

积分
363
发表于 2017-9-11 10:40:08 | 显示全部楼层
hello,按照你的结果流程重新跑了一次,找到了其中一个原因:在上一步reads计数合并表达矩阵,生成的文件中是这样的:
我没有把红色框内的信息去掉,所有导致了下面找不到差异基因。
去掉上述内容后,重新跑了一次,发现结果依然找不到差异基因。然后重新认真审查了合并表达矩阵流程,
发现了第二个错误原因——对照组和实验组搞混了,如下:

正确的合并矩阵已经是下面的:

唉~~终于搞定这一步了,粗心耶

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

1

主题

10

帖子

422

积分

中级会员

Rank: 3Rank: 3

积分
422
发表于 2017-9-5 08:38:44 | 显示全部楼层
学习一下~
回复 支持 反对

使用道具 举报

6

主题

26

帖子

363

积分

中级会员

Rank: 3Rank: 3

积分
363
发表于 2017-9-8 14:33:11 | 显示全部楼层
我照着你的步骤的代码走了一遍(从reads计数开始),到获得差异表达结果,然后统计了一下,发现上调和下调的基因数很少喔。不知当时你跑出来的结果如何呢?能参考一下吗?
以下是我分析到的结果:


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

9

主题

15

帖子

126

积分

注册会员

Rank: 2

积分
126
 楼主| 发表于 2017-9-9 03:12:27 | 显示全部楼层
蚊子 发表于 2017-9-8 14:33
我照着你的步骤的代码走了一遍(从reads计数开始),到获得差异表达结果,然后统计了一下,发现上调和下调 ...

Hi,下面是我的详细步骤及其结果。差异基因几百个吧,没有你显示的那么少。
不过,最后富集作图的时候,有些GO和KEGG没有富集出结果,估计和富集的设置条件及P值校正有关,我没有细究那步。
> dataFrame <- read.csv("/home/R/AKAP95/matrix/m_raw_count.csv", header = TRUE, sep =",")
> head(dataFrame)
                gene_id m_control1 m_control2 akap95_1 akap95_2
1  ENSMUSG00000000001.4       3338       4488     5807     5449
2 ENSMUSG00000000003.15          0          0        0        0
3 ENSMUSG00000000028.14       1683       1867     2806     2039
4 ENSMUSG00000000031.15         98        104       73       92
5 ENSMUSG00000000037.16        136        108      181      129
6 ENSMUSG00000000049.11          0          0        1        0
> countMatrix <- as.matrix(dataFrame[2:5])
> rownames(countMatrix) <- gsub("\\.\\d*", "", dataFrame$gene_id)
> head(countMatrix)
                   m_control1 m_control2 akap95_1 akap95_2
ENSMUSG00000000001       3338       4488     5807     5449
ENSMUSG00000000003          0          0        0        0
ENSMUSG00000000028       1683       1867     2806     2039
ENSMUSG00000000031         98        104       73       92
ENSMUSG00000000037        136        108      181      129
ENSMUSG00000000049          0          0        1        0
> condition <- factor(c(rep("control",2), rep("akap95",2)), levels = c("control","akap95"))
> sampleNames <- c("m_control1", "m_control2", "akap95_1", "akap95_2")
> colData <- data.frame(name= c("m_control1", "m_control2", "akap95_1", "akap95_2"), condition)
> rownames(colData) <- sampleNames
> head(colData)
                 name condition
m_control1 m_control1   control
m_control2 m_control2   control
akap95_1     akap95_1    akap95
akap95_2     akap95_2    akap95
> dds <- DESeqDataSetFromMatrix(countMatrix, colData, design= ~ condition)
> dds
class: DESeqDataSet
dim: 48820 4
metadata(1): version
assays(1): counts
rownames(48820): ENSMUSG00000000001 ENSMUSG00000000003 ... ENSMUSG00000110423
  ENSMUSG00000110424
rowData names(0):
colnames(4): m_control1 m_control2 akap95_1 akap95_2
colData names(2): name condition
> dds <- dds[rowSums(counts(dds)) > 1,]
> dds
class: DESeqDataSet
dim: 26424 4
metadata(1): version
assays(1): counts
rownames(26424): ENSMUSG00000000001 ENSMUSG00000000028 ... ENSMUSG00000110422
  ENSMUSG00000110424
rowData names(0):
colnames(4): m_control1 m_control2 akap95_1 akap95_2
colData names(2): name condition
> dds2 <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> resultsNames(dds2)
[1] "Intercept"        "conditioncontrol" "conditionakap95"
> res <- results(dds2)
> summary(res)

out of 26424 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 804, 3%
LFC < 0 (down)   : 672, 2.5%
outliers [1]     : 0, 0%
low counts [2]   : 10758, 41%
(mean count < 20)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
回复 支持 反对

使用道具 举报

9

主题

15

帖子

126

积分

注册会员

Rank: 2

积分
126
 楼主| 发表于 2017-9-11 21:10:42 | 显示全部楼层
蚊子 发表于 2017-9-11 10:40
hello,按照你的结果流程重新跑了一次,找到了其中一个原因:在上一步reads计数合并表达矩阵,生成的文件中 ...

没关系,都是多折腾几次就明白了。
回复 支持 反对

使用道具 举报

3

主题

20

帖子

156

积分

注册会员

Rank: 2

积分
156
发表于 2017-11-25 18:02:00 | 显示全部楼层
蚊子 发表于 2017-9-11 10:40
hello,按照你的结果流程重新跑了一次,找到了其中一个原因:在上一步reads计数合并表达矩阵,生成的文件中 ...

红框内容不去掉 完全没问题的,我都不去。估计就是control 和 treat 弄反了。
回复 支持 反对

使用道具 举报

1

主题

54

帖子

672

积分

高级会员

Rank: 4

积分
672
发表于 2018-1-10 09:41:19 | 显示全部楼层
没有可视化吗
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-7-22 11:10 , Processed in 0.116635 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.