library(cluster)
library(Biobase)
library(qvalue)
library(fastcluster)
options(stringsAsFactors = FALSE)
NO_REUSE = T

# try to reuse earlier-loaded data if possible
if (file.exists("all.htseq-count.unstranded.tsv.RData") && ! NO_REUSE) {
    print('RESTORING DATA FROM EARLIER ANALYSIS')
    load("all.htseq-count.unstranded.tsv.RData")
} else {
    print('Reading matrix file.')
    primary_data = read.table("/home/jforment/biovice/internal_projects/200120_mblazquez_rnaseq/04-diffexpr.Ler/all.htseq-count.unstranded.tsv", header=T, com='', row.names=1, check.names=F, sep='\t')
    primary_data = as.matrix(primary_data)
}
source("/home/jforment/install/trinityrnaseq-Trinity-v2.8.4/Analysis/DifferentialExpression/R/heatmap.3.R")
source("/home/jforment/install/trinityrnaseq-Trinity-v2.8.4/Analysis/DifferentialExpression/R/misc_rnaseq_funcs.R")
source("/home/jforment/install/trinityrnaseq-Trinity-v2.8.4/Analysis/DifferentialExpression/R/pairs3.R")
source("/home/jforment/install/trinityrnaseq-Trinity-v2.8.4/Analysis/DifferentialExpression/R/vioplot2.R")
data = primary_data
myheatcol = colorpanel(75, 'purple','black','yellow')
samples_data = read.table("/home/jforment/biovice/internal_projects/200120_mblazquez_rnaseq/00-rawdata/samples.tsv", header=F, check.names=F, fill=T)
samples_data = samples_data[samples_data[,2] != '',]
colnames(samples_data) = c('sample_name', 'replicate_name')
sample_types = as.character(unique(samples_data[,1]))
rep_names = as.character(samples_data[,2])
data = data[, colnames(data) %in% rep_names, drop=F ]
nsamples = length(sample_types)
sample_colors = rainbow(nsamples)
names(sample_colors) = sample_types
sample_type_list = list()
for (i in 1:nsamples) {
    samples_want = samples_data[samples_data[,1]==sample_types[i], 2]
    sample_type_list[[sample_types[i]]] = as.vector(samples_want)
}
sample_factoring = colnames(data)
for (i in 1:nsamples) {
    sample_type = sample_types[i]
    replicates_want = sample_type_list[[sample_type]]
    sample_factoring[ colnames(data) %in% replicates_want ] = sample_type
}
pdf("all.htseq-count.unstranded.tsv.barplot_sum_counts.pdf")
op <- par(mar = c(10,10,10,10))
barplot(colSums(data), las=2, main=paste("Sums of Frags"), ylab='', cex.names=0.7)
dev.off()
par(op)
pdf("all.htseq-count.unstranded.tsv.boxplot_log2_dist.pdf")
# set up barplot colors:
sample_cols = rainbow(nsamples)
barplot_cols = c()
for (i in 1:nsamples) {
    barplot_cols[ sample_factoring %in% sample_types[i] ] = sample_cols[i]
}
boxplot_data = data
boxplot_data[boxplot_data<10.0] = NA
boxplot_data = log2(boxplot_data+1)
num_data_points = apply(boxplot_data, 2, function(x) sum(! is.na(x)))
write.table(num_data_points, file="all.htseq-count.unstranded.tsv.feature_per_sample_count_min10.0.dat", quote=F, sep="	")
num_features_per_boxplot = 100
for(i in 1:ceiling(ncol(boxplot_data)/num_features_per_boxplot)) {
    from = (i-1)*num_features_per_boxplot+1; to = min(from+num_features_per_boxplot-1, ncol(boxplot_data));
    op <- par(mar = c(0,4,2,2), mfrow=c(2,1))
    boxplot(boxplot_data[,from:to], outline=F, main=paste('boxplot log2 >', 10.0, ', reps:', from, '-', to), xaxt='n')
    par(mar = c(7,4,2,2))
    barplot(num_data_points[from:to], las=2, main=paste('Count of features > ', 10.0, ', reps:', from, '-', to), cex.names=0.7, col=barplot_cols[from:to])
    par(op)
}
dev.off()
data = data[rowSums(data)>=10.0,]
initial_matrix = data # store before doing various data transformations
cs = colSums(data)
data = t( t(data)/cs) * 1e6;
data = log2(data+1)
sample_factoring = colnames(data)
for (i in 1:nsamples) {
    sample_type = sample_types[i]
    replicates_want = sample_type_list[[sample_type]]
    sample_factoring[ colnames(data) %in% replicates_want ] = sample_type
}
sampleAnnotations = matrix(ncol=ncol(data),nrow=nsamples)
for (i in 1:nsamples) {
  sampleAnnotations[i,] = colnames(data) %in% sample_type_list[[sample_types[i]]]
}
sampleAnnotations = apply(sampleAnnotations, 1:2, function(x) as.logical(x))
sampleAnnotations = sample_matrix_to_color_assignments(sampleAnnotations, col=sample_colors)
rownames(sampleAnnotations) = as.vector(sample_types)
colnames(sampleAnnotations) = colnames(data)
data = as.matrix(data) # convert to matrix
MA_plot = function(x, y, ...) {
    M = log( (exp(x) + exp(y)) / 2)
    A = x - y;
    res = list(x=M, y=A)
    return(res)
}
MA_color_fun = function(x,y) {
    col = sapply(y, function(y) ifelse(abs(y) >= 1, 'red', 'black')) # color 2-fold diffs
    return(col)
}
Scatter_color_fun = function(x,y) {
    col = sapply(abs(x-y), function(z) ifelse(z >= 1, 'red', 'black')) # color 2-fold diffs
    return(col)
}
for (i in 1:nsamples) {
    sample_name = sample_types[[i]]
    cat('Processing replicate QC analysis for sample: ', sample_name, "
")
    samples_want = sample_type_list[[sample_name]]
    samples_want = colnames(data) %in% samples_want
    if (sum(samples_want) > 1) {
        pdf(file=paste(sample_name, '.rep_compare.pdf', sep=''))
        d = data[,samples_want]
        initial_matrix_samples_want = initial_matrix[,samples_want]
        op <- par(mar = c(10,10,10,10))
        barplot(colSums(initial_matrix_samples_want), las=2, main=paste("Sum of Frags for replicates of:", sample_name), ylab='', cex.names=0.7)
        par(op)
        pairs3(d, pch='.', CustomColorFun=Scatter_color_fun, main=paste('Replicate Scatter:', sample_name)) # scatter plots
        pairs3(d, XY_convert_fun=MA_plot, CustomColorFun=MA_color_fun, pch='.', main=paste('Replicate MA:', sample_name)); # MA plots
        reps_cor = cor(d, method="pearson", use='pairwise.complete.obs')
        hc_samples = hclust(as.dist(1-reps_cor), method="complete")
        heatmap.3(reps_cor, dendrogram='both', Rowv=as.dendrogram(hc_samples), Colv=as.dendrogram(hc_samples), col = cm.colors(256), scale='none', symm=TRUE, key=TRUE,density.info='none', trace='none', symbreaks=F, margins=c(10,10), cexCol=1, cexRow=1, main=paste('Replicate Correlations:', sample_name) )
        dev.off()
    }
}
write.table(data, file="all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.dat", quote=F, sep='	');
if (nrow(data) < 2) { stop("

**** Sorry, at least two rows are required for this matrix.

");}
if (ncol(data) < 2) { stop("

**** Sorry, at least two columns are required for this matrix.

");}
sample_cor = cor(data, method='pearson', use='pairwise.complete.obs')
write.table(sample_cor, file="all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.sample_cor.dat", quote=F, sep='	')
sample_dist = dist(t(data), method='euclidean')
hc_samples = hclust(sample_dist, method='complete')
pdf("all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.sample_cor_matrix.pdf")
sample_cor_for_plot = sample_cor
heatmap.3(sample_cor_for_plot, dendrogram='both', Rowv=as.dendrogram(hc_samples), Colv=as.dendrogram(hc_samples), col = myheatcol, scale='none', symm=TRUE, key=TRUE,density.info='none', trace='none', symkey=FALSE, symbreaks=F, margins=c(10,10), cexCol=1, cexRow=1, cex.main=0.75, main=paste("sample correlation matrix
", "all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2") , ColSideColors=sampleAnnotations, RowSideColors=t(sampleAnnotations))
dev.off()
pdf("all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.prcomp.principal_components.pdf")
prin_comp_data = data
pca = prcomp(prin_comp_data, center = FALSE, scale. = FALSE)
pc_pct_variance = (pca$sdev^2)/sum(pca$sdev^2)
def.par <- par(no.readonly = TRUE) # save default, for resetting...
gridlayout = matrix(c(1:4),nrow=2,ncol=2, byrow=TRUE);
layout(gridlayout, widths=c(1,1));
write.table(pca$rotation, file="all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.PCA.prcomp.scores", quote=F, sep="	")
write.table(pca$x, file="all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.PCA.prcomp.loadings", quote=F, sep="	")
PCA.loadings=pca$x
PCA.scores = pca$rotation
for (i in 1:(max(3,2)-1)) {
    xrange = range(PCA.scores[,i])
    yrange = range(PCA.scores[,i+1])
    samples_want = rownames(PCA.scores) %in% sample_type_list[[sample_types[1]]]
    pc_i_pct_var = sprintf("(%.2f%%)", pc_pct_variance[i]*100)
    pc_i_1_pct_var = sprintf("(%.2f%%)", pc_pct_variance[i+1]*100)
    plot(PCA.scores[samples_want,i], PCA.scores[samples_want,i+1], xlab=paste('PC',i, pc_i_pct_var), ylab=paste('PC',i+1, pc_i_1_pct_var), xlim=xrange, ylim=yrange, col=sample_colors[1])
    for (j in 2:nsamples) {
        samples_want = rownames(PCA.scores) %in% sample_type_list[[sample_types[j]]]
        points(PCA.scores[samples_want,i], PCA.scores[samples_want,i+1], col=sample_colors[j], pch=j)
    }
    plot.new()
    legend('topleft', as.vector(sample_types), col=sample_colors, pch=1:nsamples, ncol=2)
}

par(def.par)
pcloadings_mat_vals = PCA.loadings[,1:3]
print(dim(pcloadings_mat_vals))
pcloadings_mat = matrix_to_color_assignments(pcloadings_mat_vals, col=colorpanel(256,'purple','black','yellow'), by='col')
print(dim(pcloadings_mat))
colnames(pcloadings_mat) = paste('PC', 1:ncol(pcloadings_mat))
dev.off()
gene_cor = NULL
gene_dist = dist(data, method='euclidean')
if (nrow(data) <= 1) { message('Too few genes to generate heatmap'); quit(status=0); }
hc_genes = hclust(gene_dist, method='complete')
heatmap_data = data
pdf("all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.genes_vs_samples_heatmap.pdf")
heatmap.3(heatmap_data, dendrogram='both', Rowv=as.dendrogram(hc_genes), Colv=as.dendrogram(hc_samples), col=myheatcol, scale="none", density.info="none", trace="none", key=TRUE, keysize=1.2, cexCol=1, margins=c(10,10), cex.main=0.75, main=paste("samples vs. features
", "all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2" ) , ColSideColors=sampleAnnotations)
dev.off()
pdf("all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.mean_vs_sd.pdf")
data_mean = apply(data, 1, mean)
data_sd = apply(data, 1, sd)
data_mean_vs_sd_table = cbind(data_mean, data_sd)
data_mean_vs_sd_table = as.data.frame(data_mean_vs_sd_table)
write.table(data_mean_vs_sd_table, file="all.htseq-count.unstranded.tsv.minRow10.0.CPM.log2.mean_vs_sd.dat", quote=F, sep="	")
smoothScatter(data_mean_vs_sd_table)
dev.off()