> 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)
+ }
[1] "Reading matrix file."
> 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()
null device 
          1 
> 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()
pdf 
  2 
> 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()
+     }
+ }
Processing replicate QC analysis for sample:  5X 
Processing replicate QC analysis for sample:  Ler 
Processing replicate QC analysis for sample:  Mp 
Processing replicate QC analysis for sample:  PRO 
Processing replicate QC analysis for sample:  Pa 
Processing replicate QC analysis for sample:  RGA 
Processing replicate QC analysis for sample:  Sm 
> 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 
  2 
> 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))
[1] 22917     3
> pcloadings_mat = matrix_to_color_assignments(pcloadings_mat_vals, col=colorpanel(256,'purple','black','yellow'), by='col')
> print(dim(pcloadings_mat))
[1] 22917     3
> colnames(pcloadings_mat) = paste('PC', 1:ncol(pcloadings_mat))
> dev.off()
pdf 
  2 
> 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 
  2 
> 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()
pdf 
  2 
>