###############################################################
## DEG pipeline with edgeR, DESeq2 and NOISeq #
## Modified from the template of Dr. Kristian Ullrich 2014 #
## @fabianhaas, AG Rensing, Marburg, Jan 2016 #
###############################################################
# The pipeline expects sample triplicates.
# If you don't have three samples of one condition, use face files with zero counts.
# RPKM calculation is optional.
# Additional program to sum up the tsv files: /mnt/nas_temp/home/haasf/Program/sumup_htseq_counts.pl folder file_prefix
###########################################
## please define your variables #
###########################################
saveat <- "/mnt/NAS_coruscant_datashare/haasf/madland_RNA-seq_Hoecker/DEGs_0.9"
file <- '/mnt/NAS_coruscant_datashare/haasf/madland_RNA-seq_Hoecker/31315.p.sort.counts'
# Table with the length and the GC content
LenGCTable <- "/mnt/NAS_coruscant_datashare/haasf/GeneAtlas_RNASeq/DEG_2nd_round/cosmoss_V3.3.release_2016-07-22.IsoformV31.fasta.length_GC_Contig_Daten_for_34840_1nd.csv"
# Column position of the gene name and the width/gc content (first column is gene name!)
myLenPos <- 1
myGCPos <- 2
#
# hoo <- c("","")
# sampleLists <- c()
# for(a in hoo) {
# for(b in hoo) {
# if(a!=b) {
# tmp <- c(a,b)
# sampleLists <- c(sampleLists, list(tmp))
# }
# }
# }
sampleLists <- list(c("cop_D", "cop_L"),
c("cop_D", "spa_D"),
c("cop_D", "spa_L"),
c("cop_D", "WT_D"),
c("cop_D", "WT_L"),
c("cop_L", "spa_D"),
c("cop_L", "spa_L"),
c("cop_L", "WT_D"),
c("cop_L", "WT_L"),
c("spa_D", "spa_L"),
c("spa_D", "WT_D"),
c("spa_D", "WT_L"),
c("spa_L", "WT_D"),
c("spa_L", "WT_L"),
c("WT_D", "WT_L")
)
# Do you want a RPKM value table (1/0)
# Be warned, it takes a while to do this. approx. 10 minutes for 35000 genes and 6 samples.
rpkm <- 0
# Are the data original HTSeq files? yes <- 1 ; no <- 0
htseq <- 0
# threshold settings for the R packages
#edgeR
pvtedgeR <- 0.001
#DESeq2
pvtDESeq2 <- 0.001
#NOISeq
pvtNOISeq <- 0.9
# Output folder names
edgeRDir <- "edgeR"
DESeq2Dir <- "DESeq2"
NOISeqDir <- "NOISeq"
#########################
## stop editing here!! ##
#########################
###########################################
## executiv function #
###########################################
workflow <- function(file, saveat, SampleNames, LenGCTable, myLenPos, myGCPos, rpkm, pvtedgeR, pvtDESeq2, pvtNOISeq, edgeRDir, DESeq2Dir, NOISeqDir, htseq)
{
# setwd(workingFolder) # set workspace
sampleNameI <- SampleNames[1]
sampleNameII <- SampleNames[2]
# Comparison name
cName <- paste(sampleNameI, sampleNameII, sep = " vs. ")
filename <- strsplit(file, "/")
filename <- filename[[1]][length(filename[[1]])] # long but it's working
sampleName <- gsub(", ","_",toString(SampleNames))
samples <- SampleNames
# create output dir
dir.create(saveat, showWarnings = FALSE)
outdir <- paste0(saveat,"/",sampleName)
dir.create(outdir)
# Start the log file
sink(paste0(outdir, "/DEGpipe_",cName, ".log"))
# print parameter into the log file
print("set parameters:")
print(paste0("Counts file:", file))
print(paste0("Results stored at:", saveat))
print(paste0("SampleNames: ",SampleNames))
print(paste0("LenGCTable: ",LenGCTable))
print(paste0("myLenPos: ",myLenPos))
print(paste0("myGCPos: ",myGCPos))
print(paste0("rpkm: ",rpkm))
print(paste0("htseq: ",htseq))
print(paste0("pvtedgeR: ",pvtedgeR))
print(paste0("pvtDESeq2: ",pvtDESeq2))
print(paste0("pvtNOISeq: ",pvtNOISeq))
print(paste0("edgeRDir: ",edgeRDir))
print(paste0("DESeq2Dir: ",DESeq2Dir))
print(paste0("NOISeqDir: ",NOISeqDir))
# session and package infos for the log file
print(sessionInfo())
print(Sys.info())
# load spike-in count tables
data0 <- read.table(file, header=T, sep="\t", row.names=1)
# just if the gene ids not ending with .1 (it is nessessery for NOISeqbio iteration step)
#row.names(data0) <- paste0(row.names(data0), ".1")
# remove HTSeq rows like "__no_feature"
data0 <- data0[ with(data0, !grepl('__', rownames(data0))) , ]
data0 <- data0[ with(data0, !grepl('ERCC', rownames(data0))) , ]
#data0 <- data0[ , with(data0, !grepl('X.1', colnames(data0)))]
#rename col names (library name to roman number)
# librariesName <- c("56754_WT_Naturally_3.bam.sort.fastq.unmapped.sam.sort.bam", "56753_WT_Naturally_2.bam.sort.fastq.unmapped.sam.sort.bam", "56752_WT_Naturally_1.bam.sort.fastq.unmapped.sam.sort.bam", "56751_tt_6_days_3.bam.sort.fastq.unmapped.sam.sort.bam", "56750_tt_6_days_2.bam.sort.fastq.unmapped.sam.sort.bam", "56749_tt_6_days_1.bam.sort.fastq.unmapped.sam.sort.bam", "56748_tt_0_days_3.bam.sort.fastq.unmapped.sam.sort.bam", "56747_tt_0_days_2.bam.sort.fastq.unmapped.sam.sort.bam", "56746_tt_0_days_1.bam.sort.fastq.unmapped.sam.sort.bam", "56745_WT_10_days_3.bam.sort.fastq.unmapped.sam.sort.bam", "56744_WT_10_days_2.bam.sort.fastq.unmapped.sam.sort.bam", "56743_WT_10_days_1.bam.sort.fastq.unmapped.sam.sort.bam", "56742_WT_6_days_3.bam.sort.fastq.unmapped.sam.sort.bam", "56741_WT_6_days_2.bam.sort.fastq.unmapped.sam.sort.bam", "56740_WT_6_days_1.bam.sort.fastq.unmapped.sam.sort.bam", "56739_WT_0_days_3.bam.sort.fastq.unmapped.sam.sort.bam", "56738_WT_0_days_2.bam.sort.fastq.unmapped.sam.sort.bam", "56737_WT_0_days_1.bam.sort.fastq.unmapped.sam.sort.bam", "56736_tt_T1_tt_6d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56735_tt_T1_tt_6d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56734_tt_T1_tt_6d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56733_tt_T1_tt_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56732_tt_T1_tt_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56731_tt_T1_tt_0d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56730_tt_T1_tt_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56729_tt_T1_tt_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56728_tt_T1_tt_0d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56727_WT_T1_NA_3.bam.sort.fastq.unmapped.sam.sort.bam", "56726_WT_T1_NA_2.bam.sort.fastq.unmapped.sam.sort.bam", "56725_WT_T1_NA_1.bam.sort.fastq.unmapped.sam.sort.bam", "56724_WT_T1_WT_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56723_WT_T1_WT_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56722_WT_T1_WT_0d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56721_WT_T1_WT_10d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56720_WT_T1_WT_10d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56719_WT_T1_WT_10d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56718_WT_T1_WT_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56717_WT_T1_WT_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56716_WT_T1_WT_0d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56715_WT_T1_WT_6d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56714_WT_T1_WT_6d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56713_WT_T1_WT_6d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56712_WT_T1_WT_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56711_WT_T1_WT_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56710_WT_T1_WT_0d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56709_WT_T1_WT_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56708_WT_T1_WT_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56707_WT_T1_WT_0d_1.bam.sort.fastq.unmapped.sam.sort.bam", "56706_WT_T1_tt_0d_3.bam.sort.fastq.unmapped.sam.sort.bam", "56705_WT_T1_tt_0d_2.bam.sort.fastq.unmapped.sam.sort.bam", "56704_WT_T1_tt_0d_1.bam.sort.fastq.unmapped.sam.sort.bam")
# #librariesName <- c("80484_Celery_3_day_imbibed_embryo_20C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80483_Celery_3_day_imbibed_embryo_20C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80482_Celery_3_day_imbibed_embryo_20C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80481_Celery_1_day_imbibed_embryo_20C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80480_Celery_1_day_imbibed_embryo_20C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80479_Celery_1_day_imbibed_embryo_20C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80478_Celery_5_day_imbibed_whole_fruit_29C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80477_Celery_5_day_imbibed_whole_fruit_29C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80460_Celery_3_day_imbibed_whole_fruit_20C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80459_Celery_3_day_imbibed_whole_fruit_20C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80458_Celery_1_day_imbibed_whole_fruit_20C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80457_Celery_1_day_imbibed_whole_fruit_20C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80456_Celery_1_day_imbibed_whole_fruit_20C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80455_Celery_Dry_whole_fruit_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80454_Celery_Dry_whole_fruit_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80453_Celery_Dry_whole_fruit_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80452_Celery_5_day_imbibed_whole_fruit_29C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80451_Celery_5_day_imbibed_whole_fruit_100uM_ABA_20C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80450_Celery_5_day_imbibed_whole_fruit_100uM_ABA_20C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80449_Celery_5_day_imbibed_whole_fruit_100uM_ABA_20C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80448_Celery_5_day_imbibed_whole_fruit_20C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80447_Celery_5_day_imbibed_whole_fruit_20C_rep_2_2.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80446_Celery_5_day_imbibed_whole_fruit_20C_rep_1_1.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam", "80445_Celery_3_day_imbibed_whole_fruit_20C_rep_3_3.bam.sort.fastq.unmapped.concordant_uniqsort.bam_reads.Celery_genome_GG_trinity.500bp_new.sam.sort.bam.featureCounts.p.txt")
# romanNumber <- c("T17_Dry_Naturally_Freezedryseed", "T17_Dry_Naturally_Freezedryseed.1", "T17_Dry_Naturally_Freezedryseed.2", "T16_Dry_6days_Freezeafter6dageing", "T16_Dry_6days_Freezeafter6dageing.1", "T16_Dry_6days_Freezeafter6dageing.2", "T15_Dry_0days_Freezedryseed", "T15_Dry_0days_Freezedryseed.1", "T15_Dry_0days_Freezedryseed.2", "T14_Dry_10days_Freezeafter10dageing", "T14_Dry_10days_Freezeafter10dageing.1", "T14_Dry_10days_Freezeafter10dageing.2", "T13_Dry_6days_Freezeafter6dageing", "T13_Dry_6days_Freezeafter6dageing.1", "T13_Dry_6days_Freezeafter6dageing.2", "T12_Dry_0days_Freezedryseed", "T12_Dry_0days_Freezedryseed.1", "T12_Dry_0days_Freezedryseed.2", "T11_Imbibe_6days_12h", "T11_Imbibe_6days_12h.1", "T11_Imbibe_6days_12h.2", "T10_Imbibe_6days_5h", "T10_Imbibe_6days_5h.1", "T10_Imbibe_6days_5h.2", "T9_Imbibe_0dyas_5h", "T9_Imbibe_0dyas_5h.1", "T9_Imbibe_0dyas_5h.2", "T8_Imbibe_Naturally_47h", "T8_Imbibe_Naturally_47h.1", "T8_Imbibe_Naturally_47h.2", "T7_Imbibe_Naturally_24h", "T7_Imbibe_Naturally_24h.1", "T7_Imbibe_Naturally_24h.2", "T6_Imbibe_10days_72h", "T6_Imbibe_10days_72h.1", "T6_Imbibe_10days_72h.2", "T5_Imbibe_10days_24h", "T5_Imbibe_10days_24h.1", "T5_Imbibe_10days_24h.2", "T4_Imbibe_6days_47h", "T4_Imbibe_6days_47h.1", "T4_Imbibe_6days_47h.2", "T3_Imbibe_6days_24h", "T3_Imbibe_6days_24h.1", "T3_Imbibe_6days_24h.2", "T2_Imbibe_0days_24h", "T2_Imbibe_0days_24h.1", "T2_Imbibe_0days_24h.2", "T1_Imbibe_0days_5h", "T1_Imbibe_0days_5h.1", "T1_Imbibe_0days_5h.2")
# #romanNumber <- c("emb_3d_20.2", "emb_3d_20.1", "emb_3d_20", "emb_1d_20.2", "emb_1d_20.1", "emb_1d_20", "wf_5d_29.2", "wf_5d_29.1", "wf_3d_20.1", "wf_3d_20", "wf_1d_20.2", "wf_1d_20.1", "wf_1d_20", "wf_dry.2", "wf_dry.1", "wf_dry", "wf_5d_29", "wf_5d_ABA_20.2", "wf_5d_ABA_20.1", "wf_5d_ABA_20", "wf_5d_20.2", "wf_5d_20.1", "wf_5d_20", "wf_3d_20.2")
#
#
# for(i in 1:length(librariesName)) {
# colnames(data0)[colnames(data0)==grep(librariesName[i], names(data0), value = TRUE)] <- romanNumber[i]
# }
# creats a list of the column names/samples you would like to compare
SampleNames <- makeSampleList(sampleNameI, length(grep(paste0("^",sampleNameI,"[^A-Z]"), colnames(data0))) + 1,
sampleNameII, length(grep(paste0("^",sampleNameII,"[^A-Z]"), colnames(data0))) + 1)
print(SampleNames)
## mal anschauen als ersatz für das darüber ##
## newenw <- table.new[, grepl(paste(LIST,collapse="|"), colnames(table.new))]
# exctract the samples from the data0 if needed
if(length(SampleNames) > 0) {
RTableInFiles <- data0[,SampleNames]
data <- data0[,SampleNames]
}
print(head(RTableInFiles))
# sort by column name
# data <- data[ , order(names(data))] # dangeros, the order of query and referenz might be not the same
# see how many rows and columns are in the tables
print(paste("Raw data row number: ",nrow(data)))
print(paste("Raw data column number: ",ncol(data)))
# 2) A typical differential expression analysis workflow
# 2.1) Filtering and exploratory data analysis
# filter: One single gene/spike-in must have at least 5 counts in minimum 2 libraries
filter <- apply(data, 1, function(x) length(x[x>5])>=2)
dataFiltered <- data[filter,]
# see how many rows and columns are left
print(paste("Filtered data row number: ",nrow(dataFiltered)))
print(paste("Filtered data row number: ",ncol(dataFiltered)))
# filter gene names
genes <- rownames(dataFiltered)[grep("^Pp3", rownames(dataFiltered))]
# filter spike-in names
spikes <- rownames(dataFiltered)[grep("^ERCC", rownames(dataFiltered))]
# get the column name list as.factors
xtmp <- colnames(dataFiltered)
x <- split2list(xtmp)
sampleNameList <- gsub(".[1234567890]+$", "", SampleNames)
columnPos <- makeIndexList(SampleNames)
# run ...
### edgeR ###
print(paste0("edgeR is running with ", cName))
packageDescription("edgeR")
edgeRFunction(outdir, edgeRDir, cName, RTableInFiles, pvtedgeR, sampleNameI, sampleNameII, sampleNameList)
#### DESeq2 ####
print(paste0("DESeq2 is running with ", cName))
packageDescription("DESeq2")
DESeq2Function(outdir, DESeq2Dir, cName, RTableInFiles, pvtDESeq2, sampleNameI, sampleNameII, sampleNameList)
### NOISeq ###
print(paste0("NOISeq is running with ", cName))
packageDescription("NOISeq")
NOISeqFunction(outdir, NOISeqDir, cName, RTableInFiles, pvtNOISeq, sampleNameI, sampleNameII, myLenPos, myGCPos, LenGCTable, columnPos, sampleNameList, SampleNames)
### Overlap of all three methods ###
# intersection with DEG direction check => 1, without => 0
intersection_4.file <- intersection(outdir,1,4)
write.table(intersection_4.file, file = paste0(outdir, "/", cName,"_intersection_four_methods_",nrow(intersection_4.file),".tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
# edgeR, DESeq2, NOISeq
intersection_3.file <- intersection(outdir,1,3)
write.table(intersection_3.file, file = paste0(outdir, "/", cName,"_intersection_three_methods_",nrow(intersection_3.file),".tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
sink() # End of log file
}
####################
## main functions #
####################
plotGraphics <- function(data, librariesName, cName, outdir, romanNumber)
{
library("DESeq2")
library("ggplot2")
library("RColorBrewer")
library("pheatmap")
library("genefilter")
# create sample annotation lists
# sample.name.list <- gsub(".[12]", "", librariesName[colnames(data)])
sample.name.list <- gsub(".[12]", "", romanNumber)
sample.name.uniq <- unique(unlist(sample.name.list, use.names = FALSE))
col.data.condition = data.frame(row.names = librariesName, condition = sample.name.list)
# run DESeq2 functions
dds.matrix <- DESeqDataSetFromMatrix(countData = data, colData = col.data.condition, design = ~condition) # data must be integer! (Ganze Zahlen)
colData(dds.matrix)$condition <- factor(colData(dds.matrix)$condition, levels = sample.name.uniq)
dds.matrix.deseq <- DESeq(dds.matrix)
dds.matrix.res <- results(dds.matrix.deseq)
dds.matrix.res <- dds.matrix.res[order(dds.matrix.res$padj), ]
dds.matrix.res.fdr <- dds.matrix.res[!is.na(dds.matrix.res$padj), ]
dds.matrix.res.fdr <- dds.matrix.res.fdr[dds.matrix.res.fdr$padj < 0.0001, ]
dds.matrix.rld <- rlogTransformation(dds.matrix, blind = TRUE)
# data for PCA plot
pcaData <- plotPCA(dds.matrix.rld, intgroup = "condition", returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
# color for samples
sample.colors <- c("gray60","#D73027","#4575B4", "purple", "orange","#113FF7","green", "black")
# plot the PCA
# https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html
png(filename=paste0(outdir,"PCA_all_samples_DEseq2_normalized.png"))
ggplot(
pcaData,
aes(
PC1,
PC2,
color=condition,
shape=condition)
) +
scale_shape_manual(values = c(4, 1, 2, 3, 0, 5, 6, 7)) +
geom_point(size=2.5) +
xlab(
paste0("PC1: ",percentVar[1],"% variance")
) +
ylab(
paste0("PC2: ",percentVar[2],"% variance")
) +
theme() + #, face="bold"
scale_colour_manual(
values=sample.colors # dodgerblue3
) +
ggtitle("PCA of all samples (DESeq2 normalized)") +
theme(
plot.title = element_text(lineheight=.8),
panel.background = element_rect(fill = "gray95")
)
dev.off()
# how many genes will be involved in the heatmap calculation
numbGenes <- 100
# pheatmap data
topVarGenes <- head(order(rowVars(assay(dds.matrix.rld)),decreasing=TRUE),numbGenes)
mat <- assay(dds.matrix.rld)[ topVarGenes, ]
mat <- mat - rowMeans(mat)
# pheatmap data annotation
anno <-as.data.frame(colData(dds.matrix.rld)[,c("condition")])
rownames(anno) <- colnames(data) # check out the row names of annotation
colnames(anno) <- "condition"
# pheatmap data annotation color
condition <- sample.colors
names(condition) <- sample.name.uniq
anno_colors <- list(condition = condition)
#heatmap(mat)
#pheatmap(mat)
# ftp://cran.r-project.org/pub/R/web/packages/pheatmap/pheatmap.pdf
png(filename=paste0(outdir,"HC_heatmap_all_samples_DESeq2_normalized.png"))
pheatmap(mat ,
annotation = anno,
annotation_colors = anno_colors,
fontsize = 10,
fontsize_row = 6,
fontsize_col = 10,
main = paste0("Heatmap of all samples (DESeq2 normalized) ", numbGenes," genes\n"),
annotation_names_col = F,
show_rownames = F,
cluster_rows = F,
border_color = F,
treeheight_col = 20
)
dev.off()
}
edgeRFunction <- function(odir, edgeRDir, cName, RTableInFiles, pvtedgeR, sampleNameI, sampleNameII, sampleNameList)
{
library(edgeR)
outdir <- paste0(odir, "/",edgeRDir,"/") # set outdir
dir.create(outdir, showWarnings = FALSE) # create output folder
dir.create(paste0(outdir,"plots"), showWarnings = FALSE) # create subfolder "plots"
comparison <- cName
counts <- RTableInFiles
pvt <- pvtedgeR
group <- factor(sampleNameList)
# edgeR 'normal'
y.normal <- DGEList(counts = counts, group = group)
y.normal[,1]
png(filename=paste0(outdir,"plots/",paste0(cName, " edgeR - plotMDS")))
plotMDS(y.normal, method = "bcv", main = paste0(cName, "\nedgeR - plotMDS"))
dev.off()
y.normal <- calcNormFactors(y.normal)
y.normal <- estimateCommonDisp(y.normal)
y.normal <- estimateTagwiseDisp(y.normal)
png(filename=paste0(outdir,"plots/",paste0(cName, " edgeR - plotBCV")))
plotBCV(y.normal, main = paste0(cName, "\nedgeR - plotBCV"))
dev.off()
et.normal <- exactTest(y.normal, pair = c(sampleNameI, sampleNameII))
topTags(et.normal)
# adjust.method ==> character string specifying p-value adjustment method. Possible values are "none", "BH", "fdr" (equivalent to "BH"), "BY" and "holm". See p.adjust for details.
summary(de.normal <- decideTestsDGE(et.normal, p = pvt, adjust = "BH"))
et.normal.fdr <- cbind(et.normal$table, p.adjust(et.normal$table[, 3], method = "BH"))
et.normal.fdr <- et.normal.fdr[et.normal.fdr[, 4] < pvt, ]
write.table(et.normal.fdr, file = paste0(outdir, cName, "_", nrow(et.normal.fdr), "_edgeR_fdr.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
detags.normal <- rownames(de.normal)[as.logical(de.normal)]
png(filename=paste0(outdir,"plots/",paste0(cName, " edgeR - plotSmear")))
plotSmear(et.normal, de.normal, tags = detags.normal, main = paste0(cName, "\nedgeR - plotSmear"))
abline(h = c(-2, 2), col = "blue")
dev.off()
# edgeR 'GLM'
y.glm <- DGEList(counts = counts, group = group)
design <- model.matrix(~0 + group, data = y.glm$samples)
design
y.glm <- estimateGLMCommonDisp(y.glm, design)
y.glm <- estimateGLMTrendedDisp(y.glm, design)
## Loading required package: splines
y.glm <- estimateGLMTagwiseDisp(y.glm, design)
png(filename=paste0(outdir,"plots/",paste0(cName, " edgeR - plotBCV - GLM")))
plotBCV(y.glm, main = paste0(cName, "\nedgeR - plotBCV - GLM"))
dev.off()
y.fit <- glmFit(y.glm, design)
y.lrt <- glmLRT(y.fit)
FDR <- p.adjust(y.lrt$table$PValue, method = "BH")
sum(FDR < pvt)
topTags(y.lrt)
summary(de.glm <- decideTestsDGE(y.lrt, p = pvt, adjust = "BH"))
}
DESeq2Function <- function(odir, DESeq2Dir, cName, RTableInFiles, pvtDESeq2, sampleNameI, sampleNameII, sampleNameList)
{
library(DESeq2)
library(RColorBrewer)
library(gplots)
library(vsn)
outdir <- paste0(odir, "/",DESeq2Dir,"/") # set outdir
dir.create(outdir, showWarnings = FALSE)
dir.create(paste0(outdir,"plots"), showWarnings = FALSE) # create subfolder "plots"
pvt <- pvtDESeq2
# DESeqDataSet(se, design, ignoreRank = FALSE)
# DESeqDataSetFromMatrix(countData, colData, design, tidy = FALSE, ignoreRank = FALSE, ...)
# DESeqDataSetFromHTSeqCount(sampleTable, directory = ".", design, ignoreRank = FALSE, ...)
# DESeqDataSetFromTximport(txi, colData, design, ...)
coldata = data.frame(row.names = colnames(RTableInFiles), condition = sampleNameList)
ddsHTSeq <- DESeqDataSetFromMatrix(countData = RTableInFiles, colData = coldata, design = ~ condition)
colData(ddsHTSeq)$condition <- factor(colData(ddsHTSeq)$condition, levels = c(sampleNameI, sampleNameII))
ddsHTSeq <- DESeq(ddsHTSeq)
ddsHTSeq.res <- results(ddsHTSeq)
ddsHTSeq.res <- ddsHTSeq.res[order(ddsHTSeq.res$padj), ]
head(ddsHTSeq.res)
# png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - plotMA")))
# plotMA(ddsHTSeq, main = paste0(cName, "\nDESeq2 - plotMA"), ylim = c(-2, 2))
# dev.off()
ddsHTSeq.res.fdr <- ddsHTSeq.res[!is.na(ddsHTSeq.res$padj), ]
ddsHTSeq.res.fdr <- ddsHTSeq.res.fdr[ddsHTSeq.res.fdr$padj < pvt, ]
write.table(as.matrix(ddsHTSeq.res.fdr), file = paste0(outdir, cName, "_", nrow(as.matrix(ddsHTSeq.res.fdr)), "_DESeq2_fdr.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
ddsHTSeq.rld <- rlogTransformation(ddsHTSeq, blind = TRUE)
ddsHTSeq.vsd <- varianceStabilizingTransformation(ddsHTSeq, blind = TRUE)
par(mfrow = c(1, 3))
notAllZero <- (rowSums(counts(ddsHTSeq)) > 0)
png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - SdPlot.1")))
meanSdPlot(log2(counts(ddsHTSeq, normalized = TRUE)[notAllZero, ] + 1)) # , ylim = c(0, 2.5)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - SdPlot.2"))) # , ylim = c(0, 2.5)
meanSdPlot(assay(ddsHTSeq.rld[notAllZero, ]))
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - SdPlot.3"))) # , ylim = c(0, 2.5)
meanSdPlot(assay(ddsHTSeq.vsd[notAllZero, ]))
dev.off()
hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100)
distRL <- dist(t(assay(ddsHTSeq.rld)))
mat <- as.matrix(distRL)
rownames(mat) <- colnames(mat) <- with(colData(ddsHTSeq), paste0(condition))
# png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - heatmap")))
# heatmap(mat, trace = "none", col = rev(hmcol), margin = c(13, 13))
# dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - PCA")))
print(plotPCA(ddsHTSeq.rld, intgroup = "condition"))
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " DESeq2 - DispEsts")))
plotDispEsts(ddsHTSeq, main = paste0(cName, "\nDESeq2 - plotDispEsts"))
dev.off()
}
NOISeqFunction <- function(odir, NOISeqDir, cName, RTableInFiles, pvtNOISeq, sampleNameI, sampleNameII, myLenPos, myGCPos, LenGCTable, columnPos, sampleNameList, SampleNames)
{
library(NOISeq)
outdir <- paste0(odir, "/",NOISeqDir,"/") # set outdir
dir.create(outdir, showWarnings = FALSE)
dir.create(paste0(outdir,"plots"), showWarnings = FALSE) # create subfolder "plots"
pvt <- pvtNOISeq
counts <- RTableInFiles[, columnPos]
head(counts)
# must be changed if there are no triplicates
# !! Please change this setting if needed !!
myfactors <- data.frame(Treatment = sampleNameList, TreatmentRun = SampleNames)
# sort counts by row names
mycounts <- counts[order(rownames(counts)), ]
# filter out low counts
mycounts.filter <- filtered.data(mycounts, factor = myfactors$Treatment, norm = FALSE, depth = NULL, method = 1, cv.cutoff = 100, cpm = 1)
mycounts.rowname <- rownames(mycounts.filter)
# Extra file with length and GC content
myLengthTable <- read.table(LenGCTable, sep = "\t", header = TRUE, stringsAsFactor = FALSE, row.names=1)
myLengthTable <- myLengthTable[order(rownames(myLengthTable)),] # sort by rownames
myLengthTable.filter <- myLengthTable[mycounts.rowname,]
mylength <- myLengthTable.filter[,myLenPos]
names(mylength) <- rownames(myLengthTable.filter)
mygc <- myLengthTable.filter[,myGCPos]
names(mygc) <- rownames(myLengthTable.filter)
mydata <- readData(data = mycounts.filter, length = mylength, gc = mygc, factors = myfactors)
mydata
mydata.GCbias <- dat(mydata, type = "GCbias")
#-# png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - GCbias.1")))
#-# explo.plot(mydata.GCbias)
#-# dev.off()
mydata.GCbias.group <- dat(mydata, factor = "Treatment", type = "GCbias")
#-# png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - GCbias.2")))
#-# explo.plot(mydata.GCbias.group)
#-# dev.off()
show(mydata.GCbias.group)
mydata.lengthbias <- dat(mydata, type = "lengthbias")
#-# png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - Lenbias.1")))
#-# explo.plot(mydata.lengthbias)
#-# dev.off()
mydata.lengthbias.group <- dat(mydata, factor = "Treatment", type = "lengthbias")
#-# png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - Lenbias.2")))
#-# explo.plot(mydata.lengthbias.group)
#-# dev.off()
show(mydata.lengthbias.group)
# looks different
mydata.cd <- dat(mydata, type = "cd")
#-# png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - cd")))
#-# explo.plot(mydata.cd)
#-# dev.off()
mylength
myRPKM <- NOISeq::rpkm(assayData(mydata)$exprs)
myUQUA <- NOISeq::uqua(assayData(mydata)$exprs, long = mylength, lc = 0.5, k = 0)
myTMM <- NOISeq::tmm(assayData(mydata)$exprs, long = 1000, lc = 0)
mynoiseq.rpkm <- noiseq(mydata, k = 0.5, norm = "rpkm", factor = "Treatment", pnr = 0.2, nss = 5, v = 0.02, lc = 1, replicates = "biological")
mynoiseq.nnorm <- noiseq(mydata, k = 0.5, norm = "n", factor = "Treatment", pnr = 0.2, nss = 5, v = 0.02, lc = 1, replicates = "biological")
mynoiseq.tmm <- noiseq(mydata, k = 0.5, norm = "tmm", factor = "Treatment", pnr = 0.2, nss = 5, v = 0.02, lc = 1, replicates = "biological")
# RPKM
mynoiseq.rpkm.prob <- mynoiseq.rpkm@results[[1]][mynoiseq.rpkm@results[[1]]$prob > pvt, ]
mynoiseq.rpkm.prob <- mynoiseq.rpkm.prob[!is.na(mynoiseq.rpkm.prob$prob), ]
write.table(mynoiseq.rpkm.prob, file = paste0(outdir, cName, "_", nrow(mynoiseq.rpkm.prob), "_NOISeq_rpkm_prob.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
# none
mynoiseq.nnorm.prob <- mynoiseq.nnorm@results[[1]][mynoiseq.nnorm@results[[1]]$prob > pvt, ]
mynoiseq.nnorm.prob <- mynoiseq.nnorm.prob[!is.na(mynoiseq.nnorm.prob$prob), ]
write.table(mynoiseq.nnorm.prob, file = paste0(outdir, cName,"_", nrow(mynoiseq.nnorm.prob), "_NOISeq_nnorm_prob.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
# TMM
mynoiseq.tmm.prob <- mynoiseq.tmm@results[[1]][mynoiseq.tmm@results[[1]]$prob > pvt, ]
mynoiseq.tmm.prob <- mynoiseq.tmm.prob[!is.na(mynoiseq.tmm.prob$prob), ]
write.table(mynoiseq.tmm.prob, file = paste0(outdir, cName,"_", nrow(mynoiseq.tmm.prob), "_NOISeq_tmm_prob.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
############ NOISeq bio ################
pvtnsb <- 0.995
a0per_value <- 0.99
iterations <- 5
x5 <- c(sample(10000:99999, iterations, replace=F), 12345)
x5
collector.rpkm <- NULL
collector.tmm <- NULL
# for (iter in x5) {
# # RPKM
# mynoiseqbio.rpkm <- noiseqbio(mydata, k = 0.5, norm = "rpkm", factor = "Treatment", lc = 1, r = 20, adj = 1.5, a0per = a0per_value, random.seed = iter, filter = 1)#, plot = TRUE)
# mynoiseqbio.rpkm.prob <- mynoiseqbio.rpkm@results[[1]][mynoiseqbio.rpkm@results[[1]]$prob > pvtnsb, ]
# mynoiseqbio.rpkm.prob <- mynoiseqbio.rpkm.prob[!is.na(mynoiseqbio.rpkm.prob$prob), ]
# collector.rpkm <- rbind(collector.rpkm, mynoiseqbio.rpkm.prob)
# #collector <- rbind(collector, mynoiseqbio.rpkm.prob[rownames(mynoiseqbio.rpkm.prob) %in% rownames(collector),])
# # TMM
# mynoiseqbio.tmm <- noiseqbio(mydata, k = 0.5, norm = "tmm", factor = "Treatment", lc = 1, r = 20, adj = 1.5, a0per = a0per_value, random.seed = iter, filter = 1)#, plot = TRUE)
# mynoiseqbio.tmm.prob <- mynoiseqbio.tmm@results[[1]][mynoiseqbio.tmm@results[[1]]$prob > pvtnsb, ]
# mynoiseqbio.tmm.prob <- mynoiseqbio.tmm.prob[!is.na(mynoiseqbio.tmm.prob$prob), ]
# collector.tmm <- rbind(collector.tmm, mynoiseqbio.tmm.prob)
# }
#
# # RPKM
# print(sort(row.names(collector.rpkm)))
# library(stringr)
# mynoiseqbio.rpkm.prob.overlap <- collector.rpkm[str_split_fixed(row.names(collector.rpkm), "\\.",2)[,2]>=iterations+10,]
# row.names(mynoiseqbio.rpkm.prob.overlap) <- substr(row.names(mynoiseqbio.rpkm.prob.overlap),1,str_length(row.names(mynoiseqbio.rpkm.prob.overlap))-1)
#
# write.table(mynoiseqbio.rpkm.prob.overlap, file = paste0(outdir, cName,"_", nrow(mynoiseqbio.rpkm.prob.overlap),"_NOISeqbio_rpkm_prob_seed_iteration_overlap.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
#
# #TMM
# print(sort(row.names(collector.tmm)))
# library(stringr)
# mynoiseqbio.tmm.prob.overlap <- collector.tmm[str_split_fixed(row.names(collector.tmm), "\\.",2)[,2]>=iterations+10,]
# row.names(mynoiseqbio.tmm.prob.overlap) <- substr(row.names(mynoiseqbio.tmm.prob.overlap),1,str_length(row.names(mynoiseqbio.tmm.prob.overlap))-1)
#
# write.table(mynoiseqbio.tmm.prob.overlap, file = paste0(outdir, cName,"_", nrow(mynoiseqbio.tmm.prob.overlap),"_NOISeqbio_tmm_prob_seed_iteration_overlap.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
# # none
# mynoiseqbio.nnorm.prob <- mynoiseqbio.nnorm@results[[1]][mynoiseqbio.nnorm@results[[1]]$prob > pvtnsb, ]
# mynoiseqbio.nnorm.prob <- mynoiseqbio.nnorm.prob[!is.na(mynoiseqbio.nnorm.prob$prob), ]
# write.table(mynoiseqbio.nnorm.prob, file = paste0(outdir, cName, "_NOISeqbio_nnorm_prob.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
# # TMM
# mynoiseqbio.tmm.prob <- mynoiseqbio.tmm@results[[1]][mynoiseqbio.tmm@results[[1]]$prob > pvtnsb, ]
# mynoiseqbio.tmm.prob <- mynoiseqbio.tmm.prob[!is.na(mynoiseqbio.tmm.prob$prob), ]
# write.table(mynoiseqbio.tmm.prob, file = paste0(outdir, cName, "_NOISeqbio_tmm_prob.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - noiseqbiorpkm")))
mynoiseqbio.rpkm <- noiseqbio(mydata, k = 0.5, norm = "rpkm", factor = "Treatment", lc = 1, r = 20, adj = 1.5, a0per = a0per_value, random.seed = 12345, filter = 1, plot = TRUE)
dev.off()
mynoiseqbio.rpkm.prob <- mynoiseqbio.rpkm@results[[1]][mynoiseqbio.rpkm@results[[1]]$prob > pvtnsb, ]
mynoiseqbio.rpkm.prob <- mynoiseqbio.rpkm.prob[!is.na(mynoiseqbio.rpkm.prob$prob), ]
write.table(mynoiseqbio.rpkm.prob, file = paste0(outdir, cName,"_", nrow(mynoiseqbio.rpkm.prob),"_NOISeqbio_rpkm_prob_seed_12345.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - noiseqbiontmm")))
mynoiseqbio.tmm <- noiseqbio(mydata, k = 0.5, norm = "tmm", factor = "Treatment", lc = 1, r = 20, adj = 1.5, a0per = a0per_value, random.seed = 12345, filter = 1, plot = TRUE)
dev.off()
mynoiseqbio.tmm.prob <- mynoiseqbio.tmm@results[[1]][mynoiseqbio.tmm@results[[1]]$prob > pvtnsb, ]
mynoiseqbio.tmm.prob <- mynoiseqbio.tmm.prob[!is.na(mynoiseqbio.tmm.prob$prob), ]
write.table(mynoiseqbio.tmm.prob, file = paste0(outdir, cName,"_", nrow(mynoiseqbio.tmm.prob),"_NOISeqbio_tmm_prob_seed_12345.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - noiseqbionnorm")))
mynoiseqbio.nnorm <- noiseqbio(mydata, k = 0.5, norm = "n", factor = "Treatment", lc = 1, r = 20, adj = 1.5, a0per = a0per_value, random.seed = 12345, filter = 1, plot = TRUE)
dev.off()
head(mynoiseqbio.rpkm@results[[1]])
head(mynoiseqbio.nnorm@results[[1]])
mynoiseqbio.rpkm.deg <- degenes(mynoiseqbio.rpkm, q = pvtnsb, M = NULL)
mynoiseqbio.rpkm.deg.up <- degenes(mynoiseqbio.rpkm, q = pvtnsb, M = "up")
mynoiseqbio.rpkm.deg.down <- degenes(mynoiseqbio.rpkm, q = pvtnsb, M = "down")
mynoiseqbio.nnorm.deg <- degenes(mynoiseqbio.nnorm, q = pvtnsb, M = NULL)
mynoiseqbio.nnorm.deg.up <- degenes(mynoiseqbio.nnorm, q = pvtnsb, M = "up")
mynoiseqbio.nnorm.deg.down <- degenes(mynoiseqbio.nnorm, q = pvtnsb, M = "down")
par(mfrow = c(3, 2))
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.rpkm_0.9")))
DE.plot(mynoiseq.rpkm, q = 0.9, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.rpkm_0.8")))
DE.plot(mynoiseq.rpkm, q = 0.8, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.nnorn_0.9")))
DE.plot(mynoiseq.nnorm, q = 0.9, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.nnorn_0.8")))
DE.plot(mynoiseq.nnorm, q = 0.8, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.tmm_0.9")))
DE.plot(mynoiseq.tmm, q = 0.9, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.tmm_0.8")))
DE.plot(mynoiseq.tmm, q = 0.8, graphic = "expr", log.scale = TRUE)
dev.off()
par(mfrow = c(2, 2))
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeqbio - DE.rpkm_",pvt)))
DE.plot(mynoiseqbio.rpkm, q = pvt, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeqbio - DE.rpkm_",pvtnsb)))
DE.plot(mynoiseqbio.rpkm, q = pvtnsb, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeqbio - DE.nnorm_",pvt)))
DE.plot(mynoiseqbio.nnorm, q = pvt, graphic = "expr", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeqbio - DE.nnorm_",pvtnsb)))
DE.plot(mynoiseqbio.nnorm, q = pvtnsb, graphic = "expr", log.scale = TRUE)
dev.off()
par(mfrow = c(1, 2))
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.rpkm_0.9_MD")))
DE.plot(mynoiseq.rpkm, q = 0.9, graphic = "MD", log.scale = TRUE)
dev.off()
png(filename=paste0(outdir,"plots/",paste0(cName, " NOISeq - DE.nnorm_0.9_MD")))
DE.plot(mynoiseq.nnorm, q = 0.9, graphic = "MD", log.scale = TRUE)
dev.off()
}
###########################################
## additional function for more usability #
###########################################
calcRPKM <- function(c,n,l)
{
RPKM <- (10^9* as.numeric(c)) / (as.numeric(n) * as.numeric(l))
return(RPKM)
}
RPKMMatrix <- function(InTable,directory,LenGCTable,myLenPos)
{
myLengthTable <- read.table(paste0(directory,LenGCTable), sep = "\t", header = TRUE, stringsAsFactor = FALSE, row.names = 1)
newtable <- matrix(data = c(1), nrow = nrow(InTable), ncol(InTable), byrow= TRUE)
colnames(newtable) <- colnames(InTable)
rownames(newtable) <- rownames(InTable)
sumCol <- apply(InTable,2,sum) # list of the all col sums
for(s in colnames(InTable)) {
for(r in rownames(InTable)) {
w <- InTable[r, s]
len <- myLengthTable[r, myLenPos]
sum <- as.numeric(sumCol[s])
newtable[r, s] <- calcRPKM(w, sum, len)
}
}
return(newtable)
}
# rmHTSeqLines removes the last five rows
rmHTSeqLines <- function(counts)
{
counts <- head(counts, -5)
return(counts)
}
###############
## functions ##
###############
# split list to a factor list
split2list <- function(tmp)
{
xx <- c()
for(i in 1:length(tmp)) {
if(grepl("\\.", tmp[i])) {
foo <- strsplit(tmp[i], '\\.')[[1]]
fooo <- strsplit(foo, " ")[[1]]
xx <- c(xx, fooo)
} else {
xx <- c(xx, tmp[i])
}
}
return(as.factor(xx))
}
# makeGroups
# https://github.com/drisso/RUVSeq/blob/master/R/helper.R
makeGroups <- function(xs)
{
xs <- factor(xs)
groups <- matrix(-1, nrow = length(levels(xs)), ncol = max(table(xs)))
for (i in 1:length(levels(xs))) {
idxs <- which(xs == levels(xs)[i])
groups[i,1:length(idxs)] <- idxs
}
return(groups)
}
# The function makeSampleList creates a list of the columns(samples) you would like to compare
makeSampleList <- function(A, Al, B, Bl)
{
l <- c()
if(exists("Al", mode = "numeric")) {
l <- makeSampleList.helper(A, Al, l)
l <- makeSampleList.helper(B, Bl, l)
}
return(l)
}
makeSampleList.helper <- function(v,k,l)
{
for(h in 1:k) {
if(h == 1) {
l <- c(l, v)
} else {
l <- c(l, paste0(v,".",h-1))
}
}
return(l)
}
# function creates a list of numbers out of another list. The List is the index of the given list.
makeIndexList <- function(l)
{
n <- 1
jj <- c()
for(v in l) {
jj <- c(jj, n)
n <- n + 1
}
return(jj)
}
intersection <- function(l,check,t)
{
l.list <- list.files(path = l, pattern = "*(rpkm|edgeR|DESeq2).*tsv", recursive = TRUE)
if(grepl("bio", l.list[4],fixed = TRUE)) {
l.list.noiseq <- read.table(paste0(l,"/",l.list[3]), header=T, sep="\t", row.names=1)
l.list.noiseqbio <- read.table(paste0(l,"/",l.list[4]), header=T, sep="\t", row.names=1)
} else {
l.list.noiseq <- read.table(paste0(l,"/",l.list[4]), header=T, sep="\t", row.names=1)
l.list.noiseqbio <- read.table(paste0(l,"/",l.list[3]), header=T, sep="\t", row.names=1)
}
l.list.edgeR <- read.table(paste0(l,"/",l.list[2]), header=T, sep="\t", row.names=1)
l.list.deseq2 <- read.table(paste0(l,"/",l.list[1]), header=T, sep="\t", row.names=1)
print("DEGs edgeR")
print(nrow(l.list.edgeR))
print("DEGs deseq2")
print(nrow(l.list.deseq2))
print("DEGs noiseq")
print(nrow(l.list.noiseq))
print("DEGs noiseqbio")
print(nrow(l.list.noiseqbio))
interect.list <- ""
interect.list.direct <- ""
if(t==4) {
# List of overlaps without direction check
interect.list <- intersect(rownames(l.list.edgeR), intersect(rownames(l.list.noiseq), intersect(rownames(l.list.noiseqbio), rownames(l.list.deseq2))))
print("DEGs intersection without direction check (four methods)")
print(length(interect.list))
print(interect.list)
print("DEGs intersection with direction check")
# List of overlaps with direction check
if(length(interect.list)!=0) {
interect.list.direct <- interDirect(l.list.edgeR, l.list.deseq2, l.list.noiseq, l.list.noiseqbio, interect.list)
print(length(interect.list.direct))
print(interect.list.direct)
} else {
print("0")
}
}
if(t==3) {
# List of overlaps without direction check
interect.list <- intersect(rownames(l.list.edgeR), intersect(rownames(l.list.noiseq), rownames(l.list.deseq2)))
print("DEGs intersection without direction check (three methods)")
print(length(interect.list))
print(interect.list)
print("DEGs intersection with direction check")
# List of overlaps with direction check
if(length(interect.list)!=0) {
# fake data frame with 0 rows => nrow == 0
df_fake <- data.frame(Date=as.Date(character()),
File=character(),
User=character(),
stringsAsFactors=FALSE)
interect.list.direct <- interDirect(l.list.edgeR, l.list.deseq2, l.list.noiseq, df_fake, interect.list)
print(length(interect.list.direct))
print(interect.list.direct)
} else {
print("0")
}
}
if(check == 1){
return(l.list.noiseq[interect.list.direct,])
} else {
return(l.list.noiseq[interect.list,])
}
}
interDirect <- function(e,d,n,nb,i)
{
if(nrow(nb) == 0) {
nbx <- apply(n[i,], 2, function(x) {ifelse(x > 0, "-", "+")}) # use noiseq twice to simulate noiseqbio
nby <- paste0(nbx[,"M"],rownames(nbx))
} else {
nbx <- apply(nb[i,], 2, function(x) {ifelse(x > 0, "-", "+")})
nby <- paste0(nbx[,"log2FC"],rownames(nbx))
}
ex <- apply(e[i,], 2, function(x) {ifelse(x < 0, "-", "+")})
dx <- apply(d[i,], 2, function(x) {ifelse(x < 0, "-", "+")})
nx <- apply(n[i,], 2, function(x) {ifelse(x > 0, "-", "+")}) # twisted, because edgeR and DESeq2 are calculating vice versa
ey <- paste0(ex[,"logFC"],rownames(ex))
dy <- paste0(dx[,"log2FoldChange"],rownames(dx))
ny <- paste0(nx[,"M"],rownames(nx))
interect.list.direct <- intersect(ey, intersect(dy, intersect(ny, nby)))
return(substr(interect.list.direct,2,100))
}
###########################################
## executiv part of the script #
###########################################
for(SampleNames in sampleLists) {
workflow(file, saveat, SampleNames, LenGCTable, myLenPos, myGCPos, rpkm, pvtedgeR, pvtDESeq2, pvtNOISeq, edgeRDir, DESeq2Dir, NOISeqDir, htseq)
}
print("Thank you for using the DEG pipeline :) \n")
# fin
### compare DEGs analysis ###
# used venny and Excel