Advanced Interdisciplinary Projects Lab(AIPLab) 討論區

Please login or register.

請輸入帳號, 密碼以及預計登入時間

新聞:

[開學]106學年第1學期的課程看版開張了 歡迎同學問問題-20170917

作者 主題: DEG比較  (閱讀 36 次)

admin

  • 管理員
  • Hero Member
  • *****
  • 文章: 1715
    • 檢視個人資料
DEG比較
« 於: 五月 07, 2018, 11:51:13 am »
程式碼: [Select]
#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")
程式碼: [Select]
#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))
程式碼: [Select]
#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")
程式碼: [Select]
#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
程式碼: [Select]
#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"))

已記錄
 

SimplePortal Classic 2.0.5