Artificial Intelligence Programming Lab(AIPLab) 討論區

Please login or register.

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

新聞:

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

作者 主題: Salmon+DeSeq2  (閱讀 1321 次)

admin

  • 管理員
  • Hero Member
  • *****
  • 文章: 1752
    • 檢視個人資料
Salmon+DeSeq2
« 於: 六月 24, 2019, 01:13:59 am »
install.packages("BiocManager")
install.packages("devtools")
devtools::install_github("stephenturner/annotables")
BiocManager::install("SummarizedExperiment")
BiocManager::install("DESeq2")a
BiocManager::install("pathview")
BiocManager::install("gage")
BiocManager::install("genefilter")
BiocManager::install("AnnotationDbi")
BiocManager::install("org.Hs.eg.db")


library("tximport")
library("readr")
library("tximportData")
dir <- system.file("extdata", package="tximportData")
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))
curpath <- "I:/06 Salmon/SRP068976-HCC"
setwd(curpath)
files = c("SRR3129887.sf","SRR3129888.sf","SRR3129889.sf","SRR3129890.sf","SRR3129891.sf","SRR3129892.sf", "SRR3129893.sf", "SRR3129894.sf", "SRR3129895.sf",
            "SRR3129913.sf","SRR3129921.sf", "SRR3129922.sf", "SRR3129923.sf")
cond = c("A","A","A","A","A","A","A","A", "B","B","B","B","B")
txi <- tximport(files, type="salmon", tx2gene=tx2gene)
samples <- read.table("samples.txt", header=TRUE)
samples$condition <-cond


library("DESeq2")
dds <- DESeqDataSetFromTximport(txi,
                                   colData = samples,
                                   design = ~ condition)
dds <- DESeq(dds)
res <- results(dds)
summary(res)

DESeq2::plotDispEsts(dds, main="Dispersion Estimates")
plotMA(res, main="Differentially Expressed Genes ", ylim=c(-2,2))

library("annotables")
library("AnnotationDbi")
library("org.Hs.eg.db")
eids=gsub("\\..*","",row.names(res))
res$symbol = mapIds(org.Hs.eg.db,
                     keys=eids,
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")
res$entrez = mapIds(org.Hs.eg.db,
                     keys=eids,
                     column="ENTREZID",
                     keytype="ENSEMBL",
                     multiVals="first")
res$name =   mapIds(org.Hs.eg.db,
                     keys=eids,
                     column="GENENAME",
                     keytype="ENSEMBL",
                     multiVals="first")
write.csv(as.data.frame(res),file="SRP068976-HCC-de.csv")
已記錄
 

SimplePortal Classic 2.0.5