library(massdataset)
data('expression_data',package='massdataset');data('sample_info',package='massdataset');data('variable_info',package='massdataset')
object<-create_mass_dataset(expression_data,sample_info,variable_info)
out<-'tutorial-figures'
dir.create(out,recursive=TRUE,showWarnings=FALSE)
x<-as.matrix(extract_expression_data(object))
png(file.path(out,'missingness.png'),width=1600,height=760,res=150)
par(mfrow=c(1,2),mar=c(8,4.8,3,1))
barplot(colMeans(is.na(x))*100,names.arg=colnames(x),las=2,col='#5365a5',border=NA,ylab='Missing intensity cells (%)',main='Packaged example: 1,000 features × 8 samples')
hist(rowMeans(is.na(x))*100,breaks=seq(0,100,by=12.5),col='#6e9d9e',border='white',main='Feature missingness',xlab='Missing samples per feature (%)',ylab='Number of features')
dev.off()
subject<-sample_info$sample_id[sample_info$class=='Subject'];z<-x[,subject,drop=FALSE];keep<-which(rowSums(is.finite(z))>=2)
filtered<-object[keep,subject];imputed<-masscleaner::impute_mv(filtered,method='median');normalized<-masscleaner::normalize_data(imputed,method='median')
y<-as.matrix(extract_expression_data(imputed));n<-as.matrix(extract_expression_data(normalized))
png(file.path(out,'normalization.png'),width=1600,height=760,res=150)
par(mfrow=c(1,2),mar=c(8,4.8,3,1))
boxplot(log10(y+1),las=2,col='#b9c2e2',outline=FALSE,main='After filtering and median imputation',ylab='log10(intensity + 1)')
boxplot(log10(n+1),las=2,col='#92bdbe',outline=FALSE,main='After median normalization',ylab='log10(intensity + 1)')
dev.off()
ok<-apply(n,1,function(v) all(is.finite(v)) && sd(v)>0);pca<-prcomp(t(n[ok,,drop=FALSE]),center=TRUE,scale.=TRUE);pct<-100*pca$sdev^2/sum(pca$sdev^2)
png(file.path(out,'pca.png'),width=1000,height=800,res=150)
par(mar=c(5,5,4,2));plot(pca$x[,1],pca$x[,2],xlim=extendrange(pca$x[,1],f=.15),ylim=extendrange(pca$x[,2],f=.18),pch=19,col='#303c78',cex=1.4,xlab=sprintf('PC1 (%.1f%%)',pct[1]),ylab=sprintf('PC2 (%.1f%%)',pct[2]),main='Four Subject samples: descriptive PCA');text(pca$x[,1],pca$x[,2],labels=rownames(pca$x),pos=3,cex=.85);abline(h=0,v=0,lty=3,col='gray80');dev.off()
writeLines(capture.output(sessionInfo()),file.path(out,'session-info.txt'))
cat('rows retained:',length(keep),' missing before:',sum(is.na(extract_expression_data(filtered))),' after:',sum(is.na(y)),'\n')