setwd('xx_result')

library("DESeq2")

directory <-'xx/06_featurecounts'

directory

sampleFiles <- grep("wy",list.files(directory),value=TRUE)

sampleFiles

sampleCondition <- c("Het","Het","KO","KO","KO")

sampleCondition

sampleTable <- data.frame(sampleName= sampleFiles,fileName = sampleFiles,condition = sampleCondition)

sampleTable

dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,directory = directory,design= ~ condition)

dds

dds <- dds [ rowSums(counts(dds)) > 1, ]

#PCA#

rld<-rlog(dds)

plotPCA(rld)

dds<-DESeq(dds)

res <- results(dds)

head(res)

summary(res)

resOrdered <- res[order(res$padj),]

resOrdered=as.data.frame(resOrdered)

head(resOrdered)

resOrdered=na.omit(resOrdered)

DEmiRNA=resOrdered[abs(resOrdered$log2FoldChange)>log2(1.5) & resOrdered$padj <0.01 ,]

head(resOrdered)

write.csv(resOrdered,"hisat2_samtools_htseq_DESeq2_output.csv")

#MA圖#

library("geneplotter")

plotMA(res,main="DESeq2",ylim=c(-2,2))

#heatmap#

select<-order(rowMeans(counts(dds,normalized=TRUE)),decreasing = TRUE)[1:666]

nt<-normTransform(dds)

log2.norm.counts<-assay(nt)[select,]

df<-as.data.frame(colData(dds))

library(pheatmap)

pheatmap(log2.norm.counts,cluster_rows = TRUE,show_rownames = FALSE,cluster_cols = TRUE,annotation_col = df)