|
#Merge HT-Seq count files
basedir <- "D:/SRP073050-Monocytes-counts"
setwd(basedir)
sample = c("S080","S081","S082","S083","S084","S085")
for (i in 1:length(sample) ) {
filename <- paste(sample[i], '.counts', sep = "", collapse = "")
exp_table <-read.table(filename, header = F, stringsAsFactors = FALSE)
if (i==1){#read only once
df<-read.table(filename, header = F, stringsAsFactors = FALSE)
names(df) <- c("Gene", sample[i])
}
else{
df2<-read.table(filename, header = F, stringsAsFactors = FALSE)
df[sample[i]] <-df2[2]
}
}
write.csv(df, file = "SRP073050.csv")
#Run Deseg2
library("DESeq2")
files=paste0(sample, rep(".counts", length(sample)))
cond = c("A","A","A","B","A","A")
sTable = data.frame(sampleName = files, fileName = files, condition = cond)
dds <-DESeqDataSetFromHTSeqCount(sampleTable=sTable, directory = "D:/SRP073050-Monocytes-counts", design = ~condition)
dds <- DESeq(dds)
res <- results(dds)
summary(res)
write.csv(as.data.frame(res),file="SRP073050-de.csv")
DESeq2::plotDispEsts(dds, main="Dispersion Estimates")
plotMA(res, main="Differentially Expressed Genes ", ylim=c(-2,2))
#Transform 這個步驟會比較久.....
library("annotables")
library("AnnotationDbi")
library("org.Hs.eg.db")
res$symbol = mapIds(org.Hs.eg.db,
keys=rownames(res),
column="SYMBOL",
keytype="SYMBOL",
multiVals="first")
res$entrez = mapIds(org.Hs.eg.db,
keys=row.names(res),
column="ENTREZID",
keytype="SYMBOL",
multiVals="first")
res$name = mapIds(org.Hs.eg.db,
keys=row.names(res),
column="GENENAME",
keytype="SYMBOL",
multiVals="first")
#GAG
library("dplyr")
library(pathview)
library(gage)
library(gageData)
data(kegg.sets.hs)
data(sigmet.idx.hs)
foldchanges = res$log2FoldChange
names(foldchanges) = res$entrez
kegg.sets.hs = kegg.sets.hs[sigmet.idx.hs]
# Get the results
keggres = gage(foldchanges, gsets=kegg.sets.hs, same.dir=TRUE)
# Look at both up (greater), down (less), and statatistics.
lapply(keggres, head)
# Get the pathways
keggrespathways = data.frame(id=rownames(keggres$greater), keggres$greater) %>%
tbl_df() %>%
filter(row_number()<=5) %>%
.$id %>%
as.character()
keggrespathways
#Plot KEGG pathways
detach("package:dplyr", unload=TRUE) #dplyr會和sapply衝突
keggresids = substr(keggrespathways, start=1, stop=8)
# Define plotting function for applying later
plot_pathway = function(pid) pathview(gene.data=foldchanges, pathway.id=pid, species="hsa", new.signature=FALSE)
# plot multiple pathways (plots saved to disk and returns a throwaway list object)
tmp = sapply(keggresids, function(pid) pathview(gene.data=foldchanges, pathway.id=pid, species="hsa"))
|
列出文章 - admin
