.Rhistory

library(rtracklayer)
combn(c(1,2,3),2)
as.list(combn(c(1,2,3),2))
combn(c(1,2,3),2)[,1]
library(dplyr)
?left_join
args <- c("--count_matrix_dir", "./tests", "--gff_file_dir", "./tests", "--metadata", "./tests/metadata.csv")
suppressPackageStartupMessages(library("DEXSeq"))
#load files
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
flattenedFile <- list.files(args[grep("--gff_file_dir", args)+1], pattern = "gff$", full.names = TRUE)
sampleTable <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
df <- read.table(countFiles[1], sep = "\t")
View(df)
head(df)
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
View(df)
df <- full_join(df, tmp.df, by = "V1")
for (x in countFiles[2:length(countFiles)]){
tmp.df <- read.table(x, sep = "\t")
tmp.df <- tmp.df[!tmp.df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
df <- full_join(df, tmp.df, by = "V1")
}
View(df)
basename(countFiles[1])
sub("_htseq_count.csv","",basename(countFiles[1]))
df <- read.table(countFiles[1], sep = "\t")
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(df) <- c("V1", sub("_htseq_count.csv","",basename(countFiles[1])))
for (x in countFiles[2:length(countFiles)]){
tmp.df <- read.table(x, sep = "\t")
tmp.df <- tmp.df[!tmp.df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(tmp.df) <- c("V1", sub("_htseq_count.csv","",basename(x)))
df <- full_join(df, tmp.df, by = "V1")
}
View(df)
View(df)
df[is.na(df)] <- 0
View(df)
rownames(df) <- df$V1
df$V1 <- NULL
View(df)
DEXSeqDataSet()
?DEXSeqDataSet
#construct an DEXSeqDataSet object
DEXSeqDataSet( df, sampleTable, design= ~ sample + exon + condition:exon)
#construct an DEXSeqDataSet object
DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon)
countData <- matrix( rpois(10000, 100), nrow=1000 )
sampleData <- data.frame(
condition=rep( c("untreated", "treated"), each=5 ) )
design <- formula( ~ sample + exon + condition:exon )
groupID <- rep(
paste0("gene", 1:10),
each= 100 )
featureID <- rep(
paste0("exon", 1:100),
times= 10 )
featureID
groupID
library(GenomicRanges)
library(GenomicFeatures)
library(GenomicAlignments)
makeTxDbFromGFF(file = flattenedFile, format = "auto")
flattenedFile <- list.files(args[grep("--gff_file_dir", args)+1], pattern = "gtf$", full.names = TRUE)
makeTxDbFromGFF(file = flattenedFile, format = "auto")
args
flattenedFile <- list.files(args[grep("--gff_file_dir", args)+1], pattern = "gtf$", full.names = TRUE)
flattenedFile
makeTxDbFromGFF(file = "./tests/test.gtf", format = "auto")
makeTxDbFromGFF
txdb <- makeTxDbFromGFF(file = "./tests/test.gtf", format = "auto")
as.data.frame(txdb)
txdb@.xData
txdb@.xData$.self
columns(txdb)
keytypes(txdb)
select(txdb, keys = keys, columns="TXNAME", keytype="GENEID")
select(txdb, keys = rownames(df), columns="GENEID", keytype="select")
rownames(df)
select(txdb, keys = rownames(df), columns="GENEID", keytype="TXNAME")
sub("....$", "", rownames(df))
select(txdb, keys = sub("....$", "", rownames(df)), columns="GENEID", keytype="TXNAME")
#construct an DEXSeqDataSet object
DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = rownames(df), groupID = sub("....$", "", rownames(df)))
strsplit(":",rownames(df))
?strsplit
strsplit(rownames(df),":")
as.data.frame(strsplit(rownames(df),":"))
t(as.data.frame(strsplit(rownames(df),":")))
#construct an DEXSeqDataSet object
DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[,2], groupID = ids[,1])
ids <- t(as.data.frame(strsplit(rownames(df),":")))
#construct an DEXSeqDataSet object
DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[,2], groupID = ids[,1])
?DEXSeqDataSet
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[,2], groupID = ids[,1])
as.data.frame(dxd)
counts(dxd)
df
sampleTable
df
counts(dxd)
dxd
colnames(dxd) <- colnames(df)
colData(dxd)
# normalisation
dxd = estimateSizeFactors(dxd)
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd)
# test for each exon in each gene
dxd = testForDEU( dxd)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, fitExpToVar="condition")
countData <- matrix( rpois(10000, 100), nrow=1000 )
sampleData <- data.frame(
condition=rep( c("untreated", "treated"), each=5 ) )
design <- formula( ~ sample + exon + condition:exon )
groupID <- rep(
paste0("gene", 1:10),
each= 100 )
featureID <- rep(
paste0("exon", 1:100),
times= 10 )
DEXSeqDataSet( countData, sampleData, design,
featureID, groupID )
args <- c("--count_matrix_dir", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/output/hisat2_samtools_htseq/", "--gtf_file_dir", "/project/data/rnaseq/GenomeIndex/Homo_Sapiens/Homo_sapiens.GRCh38.95.gtf", "--metadata", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/output/metadata.csv")
#load files
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
flattenedFile <- list.files(args[grep("--gtf_file_dir", args)+1], pattern = "gtf$", full.names = TRUE)
sampleTable <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
sampleTable
args <- c("--count_matrix_dir", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/output/hisat2_samtools_htseq/", "--gtf_file_dir", "/project/data/rnaseq/GenomeIndex/Homo_Sapiens/Homo_sapiens.GRCh38.95.gtf", "--metadata", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/metadata.csv")
sampleTable <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
txdb <- makeTxDbFromGFF(file = "./tests/test.gtf", format = "auto")
txdb
#load files
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
flattenedFile <- list.files(args[grep("--gtf_file_dir", args)+1], pattern = "gtf$", full.names = TRUE)
flattenedFile
args[grep("--gtf_file_dir", args)+1]
flattenedFile <- args[grep("--gtf_file_dir", args)+1]
flattenedFile
txdb <- makeTxDbFromGFF(file = "./tests/test.gtf", format = "auto")
args <- c("--count_matrix_dir", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/output/hisat2_samtools_htseq/", "--gtf_file_dir", "/project/data/rnaseq/GenomeIndex/Homo_Sapiens/Homo_sapiens.GRCh38.95.gtf", "--metadata", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/metadata.csv")
txdb
flattenedFile <- args[grep("--gtf_file_dir", args)+1]
flattenedFile
txdb <- makeTxDbFromGFF(file = flattenedFile, format = "auto")
if( "--condition" %in% args ){
condition.idx <- grep("--condition", args)
condition <- args[ condition.idx + 1 ]
colnames(sampleTable) <- sub(condition, "condition",colnames(sampleTable))
}
df <- read.table(countFiles[1], sep = "\t")
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(df) <- c("V1", sub("_htseq_count.csv","",basename(countFiles[1])))
df
df <- read.table(countFiles[1], sep = "\t")
View(df)
df <- read.table(countFiles[1], sep = " ")
df <- read.table(countFiles[1], sep = "\t")
df <- read.table(countFiles[1], sep = "  ")
df <- read.table(countFiles[1], sep = "  ", skip = 5)
df <- read.table(countFiles[1], sep = "\t", skip = 5)
View(df)
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(df) <- c("V1", sub("_htseq_count.csv","",basename(countFiles[1])))
for (x in countFiles[2:length(countFiles)]){
tmp.df <- read.table(x, sep = "\t", skip = 5)
tmp.df <- tmp.df[!tmp.df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(tmp.df) <- c("V1", sub("_htseq_count.csv","",basename(x)))
df <- full_join(df, tmp.df, by = "V1")
}
df[is.na(df)] <- 0
rownames(df) <- df$V1
df$V1 <- NULL
View(df)
ids <- t(as.data.frame(strsplit(rownames(df),":")))
ids <- as.data.frame(strsplit(rownames(df),":"))
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[1,], groupID = ids[2,])
ids <- t(as.data.frame(strsplit(rownames(df),":")))
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[,1], groupID = ids[,2])
View(ids)
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[,1], groupID = ids[,2])
# normalisation
dxd = estimateSizeFactors(dxd)
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd)
threads <- 4
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# Set up workers
BPPARAM = MultiCoreParam(workers = as.numeric(threads))
?MultiCoreParam
??MultiCoreParam
suppressPackageStartupMessages(library(BiocParallel))
# Set up workers
BPPARAM = MultiCoreParam(workers = as.numeric(threads))
# Set up workers
BPPARAM = BiocParallel::MulticoreParam(workers = as.numeric(threads))
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# test for each exon in each gene
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM)
View(df)
sampleTable
countData <- matrix( rpois(10000, 100), nrow=1000 )
sampleData <- data.frame(
condition=rep( c("untreated", "treated"), each=5 ) )
design <- formula( ~ sample + exon + condition:exon )
groupID <- rep(
paste0("gene", 1:10),
each= 100 )
featureID <- rep(
paste0("exon", 1:100),
times= 10 )
DEXSeqDataSet( countData, sampleData, design,
featureID, groupID )
countData <- matrix( rpois(10000, 100), nrow=1000 )
sampleData <- data.frame(
condition=rep( c("untreated", "treated"), each=5 ) )
design <- formula( ~ sample + exon + condition:exon )
groupID <- rep(
paste0("gene", 1:10),
each= 100 )
featureID <- rep(
paste0("exon", 1:100),
times= 10 )
dxd <- DEXSeqDataSet( countData, sampleData, design,
featureID, groupID )
# normalisation
dxd = estimateSizeFactors( dxd )
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# test for each exon in each gene
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM)
dxd
counts(dxd)
View(counts(dxd))
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ condition + exon + condition:exon, featureID = ids[,1], groupID = ids[,2])
args
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
flattenedFile <- args[grep("--gtf_file_dir", args)+1]
sampleTable <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
if( "--condition" %in% args ){
condition.idx <- grep("--condition", args)
condition <- args[ condition.idx + 1 ]
colnames(sampleTable) <- sub(condition, "condition",colnames(sampleTable))
}
df <- read.table(countFiles[1], sep = "\t", skip = 5)
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(df) <- c("V1", sub("_htseq_count.csv","",basename(countFiles[1])))
for (x in countFiles[2:length(countFiles)]){
tmp.df <- read.table(x, sep = "\t", skip = 5)
tmp.df <- tmp.df[!tmp.df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(tmp.df) <- c("V1", sub("_htseq_count.csv","",basename(x)))
df <- full_join(df, tmp.df, by = "V1")
}
df[is.na(df)] <- 0
rownames(df) <- df$V1
df$V1 <- NULL
ids <- t(as.data.frame(strsplit(rownames(df),":")))
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ sample + exon + condition:exon, featureID = ids[,1], groupID = ids[,2])
# normalisation
dxd = estimateSizeFactors( dxd )
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# test for each exon in each gene
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, fitExpToVar="condition")
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM, fitExpToVar="condition")
dxd_norm <- counts(dxd,normalized=T)
# create a results object
dxr = DEXSeqResults( dxd )
dxd_norm2 <- data.frame("name"=rownames(dxd_norm),as.data.frame(dxd_norm))
View(dxr)
dxr_dataframe = as.data.frame(dxr)
View(dxr_dataframe)
ids
head(ids)
#construct an DEXSeqDataSet object
dxd <- DEXSeqDataSet( df, sampleTable, design= ~ sample + exon + condition:exon, featureID = ids[,2], groupID = ids[,1])
# normalisation
dxd = estimateSizeFactors( dxd )
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# test for each exon in each gene
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM, fitExpToVar="condition")
dxd_norm <- counts(dxd,normalized=T)
# create a results object
dxr = DEXSeqResults( dxd )
dxr_dataframe = as.data.frame(dxr)
View(dxr)
View(dxr_dataframe)
head(dxr_dataframe)
head(dxd_norm2)
?estimateExonFoldChanges
args
suppressPackageStartupMessages(library("DEXSeq"))
suppressPackageStartupMessages(library("BiocParallel"))
library(GenomicRanges)
library(GenomicFeatures)
library(GenomicAlignments)
#load files
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
metadata <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
if( "--condition" %in% args ){
condition.idx <- grep("--condition", args)
condition <- args[ condition.idx + 1 ]
colnames(metadata) <- sub(condition, "condition",colnames(metadata))
}
df <- read.table(countFiles[1], sep = "\t", skip = 5)
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(df) <- c("V1", sub("_htseq_count.csv","",basename(countFiles[1])))
for (x in countFiles[2:length(countFiles)]){
tmp.df <- read.table(x, sep = "\t", skip = 5)
tmp.df <- tmp.df[!tmp.df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(tmp.df) <- c("V1", sub("_htseq_count.csv","",basename(x)))
df <- full_join(df, tmp.df, by = "V1")
}
df[is.na(df)] <- 0
rownames(df) <- df$V1
df$V1 <- NULL
comb <- combn(unique(as.character(metadata[,"condition"])), 2)
for(i in 1:ncol(comb)){
metadata.f <- metadata[metadata$condition %in% comb[,i],]
count.f <- df[,rownames(metadata.f)]
#construct an DEXSeqDataSet object
ids <- t(as.data.frame(strsplit(rownames(count.f),":")))
dxd <- DEXSeqDataSet(count.f, metadata, design= ~ sample + exon + condition:exon, featureID = ids[,2], groupID = ids[,1])
if("--threads" %in% args){
threads.idx <- grep("--threads", args)
threads <- args[ threads.idx + 1 ]
# Set up workers
BPPARAM = BiocParallel::MulticoreParam(workers = as.numeric(threads))
# normalisation
dxd = estimateSizeFactors( dxd )
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# test for each exon in each gene
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM, fitExpToVar="condition")
} else {
# normalisation
dxd = estimateSizeFactors(dxd)
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd)
# test for each exon in each gene
dxd = testForDEU( dxd)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, fitExpToVar="condition")
}
# create a results object
contrast <- gsub(".$","",paste0(paste0(unique(metadata.f$condition)),sep="-", collapse = ""))
dxr = DEXSeqResults( dxd )
dxr_dataframe = as.data.frame(dxr)
dxr_dataframe <- data.frame("name"=rownames(dxr_dataframe),dxr_dataframe)
#write.csv(dxr_dataframe, paste0(contrast,"_DEE_results.csv"), row.names = FALSE)
dxd_norm <- counts(dxd,normalized=T)
dxd_norm2 <- data.frame("name"=rownames(dxd_norm),as.data.frame(dxd_norm))
#write.csv(dxd_norm2, paste0(contrast,"_norm_count.csv"), row.names = FALSE)
}
history()
args <- c("--count_matrix_dir", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/output/hisat2_samtools_htseq/", "--gtf_file_dir", "/project/data/rnaseq/GenomeIndex/Homo_Sapiens/Homo_sapiens.GRCh38.95.gtf", "--metadata", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/metadata.csv")
#load files
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
metadata <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
countFiles
metadata
args <- c("--count_matrix_dir", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/output/hisat2_samtools_htseq/", "--gtf_file_dir", "/project/data/rnaseq/GenomeIndex/Homo_Sapiens/Homo_sapiens.GRCh38.95.gtf", "--metadata", "/project/data/rnaseq/Data/51c0bbfd-5f53-4f56-be15-f18dbc221b05/metadata.csv", "--threads", "4")
#load files
countFiles <- list.files(args[grep("--count_matrix_dir", args)+1], pattern = "htseq_count.csv$", full.names = TRUE)
metadata <- read.table(args[grep("--metadata", args)+1], row.names = 1, header = TRUE, sep = ",")
if( "--condition" %in% args ){
condition.idx <- grep("--condition", args)
condition <- args[ condition.idx + 1 ]
colnames(metadata) <- sub(condition, "condition",colnames(metadata))
}
df <- read.table(countFiles[1], sep = "\t", skip = 5)
df <- df[!df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(df) <- c("V1", sub("_htseq_count.csv","",basename(countFiles[1])))
for (x in countFiles[2:length(countFiles)]){
tmp.df <- read.table(x, sep = "\t", skip = 5)
tmp.df <- tmp.df[!tmp.df$V1 %in% c("_ambiguous", "_ambiguous_readpair_position", "_empty", "_lowaqual", "_notaligned"),]
colnames(tmp.df) <- c("V1", sub("_htseq_count.csv","",basename(x)))
df <- full_join(df, tmp.df, by = "V1")
}
View(df)
df[is.na(df)] <- 0
rownames(df) <- df$V1
df$V1 <- NULL
View(df)
comb <- combn(unique(as.character(metadata[,"condition"])), 2)
for(i in 1:ncol(comb)){
metadata.f <- metadata[metadata$condition %in% comb[,i],]
count.f <- df[,rownames(metadata.f)]
#construct an DEXSeqDataSet object
ids <- t(as.data.frame(strsplit(rownames(count.f),":")))
dxd <- DEXSeqDataSet(count.f, metadata, design= ~ sample + exon + condition:exon, featureID = ids[,2], groupID = ids[,1])
if("--threads" %in% args){
threads.idx <- grep("--threads", args)
threads <- args[ threads.idx + 1 ]
# Set up workers
BPPARAM = BiocParallel::MulticoreParam(workers = as.numeric(threads))
# normalisation
dxd = estimateSizeFactors( dxd )
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
# test for each exon in each gene
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM, fitExpToVar="condition")
} else {
# normalisation
dxd = estimateSizeFactors(dxd)
# dispertion estimation. estimate the variability of the data not explained by the biological variation.
dxd = estimateDispersions( dxd)
# test for each exon in each gene
dxd = testForDEU( dxd)
# calculate relative exon usage fold change
dxd = estimateExonFoldChanges(dxd, fitExpToVar="condition")
}
# create a results object
contrast <- gsub(".$","",paste0(paste0(unique(metadata.f$condition)),sep="-", collapse = ""))
dxr = DEXSeqResults( dxd )
dxr_dataframe = as.data.frame(dxr)
dxr_dataframe <- data.frame("name"=rownames(dxr_dataframe),dxr_dataframe)
#write.csv(dxr_dataframe, paste0(contrast,"_DEE_results.csv"), row.names = FALSE)
dxd_norm <- counts(dxd,normalized=T)
dxd_norm2 <- data.frame("name"=rownames(dxd_norm),as.data.frame(dxd_norm))
#write.csv(dxd_norm2, paste0(contrast,"_norm_count.csv"), row.names = FALSE)
}
head(dxr_dataframe)
colnames(dxr_dataframe)
colnames(dxr_dataframe)
length(c("name", "groupID", "featureID","exonBaseMean","dispersion","stat","pvalue","padj","Normal","Tumour","log2fold_Tumour_Normal","genomicData"))
dxr_dataframe <- dxr_dataframe[,1:12]
colnames(dxr_dataframe)
library("ballgown")
?expr
data(bg)
norm_count <- gexpr(bg)
View(norm_count)