-
Notifications
You must be signed in to change notification settings - Fork 0
/
deseq2
48 lines (35 loc) · 1.45 KB
/
deseq2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
setwd('G:/博士/project_qu_chiip8/')
df <- read.table("rna_8.count",header = TRUE,row.names = 1) #读入数据文件
colnames(df)<-c('142','143','145','146','149','150')
#df_fh_x=df[,2:5]
#df_fh_y=df[, c(2,3,6,7)]
condition <- factor(c("control","control","control","case","case",'case'))
metadata <- data.frame(sample_id = colnames(df),
condition = condition
)
dds <- DESeqDataSetFromMatrix(df, metadata, design= ~ condition )
head(dds) #查
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds,contrast = c("condition","case","control"))
# summary一下,看一下结果的概要信息
summary(res)
# 获取padj(p值经过多重校验校正后的值)小于0.05,表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
table(res$padj<0.05) #取P值小于0.05的结果
res <- res[order(res$padj),]
diff_df <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene <- row.names(diff_df)
write.csv(diff_df,file= "diff_df_fdr005_fc2.csv",quote = FALSE)
write.csv(diff_gene,file= "diff_gene_fdr005_fc2.csv",quote = FALSE)
#heatmap of degene
name <- read.csv("fh_x_def_padj0.01_fc2.csv")
rn=name[,1]
pheatmap::pheatmap(df[rn,],
fontsize = 8,
show_colnames =T,
show_rownames = F,
cluster_cols = T,
width = 8,
height = 6,
angle_col=45)
#----------pca