Artificial Intelligence Programming Lab(AIPLab) 討論區

105學年(下)課程2017Spring => 105(下)中西醫結合的大數據分析 => R packages => 主題作者是: admin 於 六月 24, 2019, 01:13:59 am



主題: Salmon+DeSeq2
作者: admin六月 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")