Created
October 15, 2019 14:52
-
-
Save cplaisier/1790cd5b3178d82ac0092b9d75383578 to your computer and use it in GitHub Desktop.
Example code form Systems Biology of Disease course
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #################################################### | |
| ## Systems Biology of Disease: Intro to LUSC Data ## | |
| ## ______ ______ __ __ ## | |
| ## /\ __ \ /\ ___\ /\ \/\ \ ## | |
| ## \ \ __ \ \ \___ \ \ \ \_\ \ ## | |
| ## \ \_\ \_\ \/\_____\ \ \_____\ ## | |
| ## \/_/\/_/ \/_____/ \/_____/ ## | |
| ## @Developed by: Plaisier Lab ## | |
| ## (https://plaisierlab.engineering.asu.edu/) ## | |
| ## Arizona State University ## | |
| ## 242 ISTB1, 550 E Orange St ## | |
| ## Tempe, AZ 85281 ## | |
| ## @Author: Chris Plaisier ## | |
| ## @License: GNU GPLv3 ## | |
| #################################################### | |
| # Set working directory | |
| setwd('C:/Users/cplaisie/Dropbox (ASU)/ASU/Courses/BME_598_494_SysBioDisease/8_22_2019/Fall_2019/Week6/LUSC_ConsensusClustering') | |
| ############################################## | |
| ### Load Discovery and Testing Cohort Data ### | |
| ############################################## | |
| # Data from TCGA Lung squamous cell carcinoma: https://tcga-data.nci.nih.gov/tcga/tcgaCancerDetails.jsp?diseaseType=LUSC&diseaseName=Lung%20squamous%20cell%20carcinoma | |
| gexp = t(na.omit(t(read.csv('data/lusc_mRNA.csv',header=T,row.names=1)))) | |
| d1 = gexp[,-1] # Get rid of Entrez ID in second column | |
| # Size of LUSC RNA-seq dataset | |
| dim(d1) | |
| # Snippet of data in the RNA-seq dataset | |
| head(d1[,1:20]) | |
| # Boxplots of array expression values to determine normlization | |
| pdf('boxplot_lusc_expression.pdf',width=20,height=5) | |
| boxplot(d1) | |
| dev.off() | |
| # Hold out 175 samples for testing cohort | |
| holdOut = read.table('data/holdOuts.txt')[,1] | |
| ho1 = d1[,holdOut] | |
| d1 = d1[,which(!(colnames(d1)%in%holdOut))] | |
| dim(d1) | |
| dim(ho1) | |
| ####################################################################### | |
| ### Select top 2,000 features using Median Absolute Deviation (MAD) ### | |
| ####################################################################### | |
| mads = apply(d1,1,mad) | |
| # Plot histogram of MAD before feature selection | |
| pdf('histogram_median_absolute_deviations_all_genes.pdf') | |
| hist(mads,xlab='Median Absolute Deviation (MAD)') | |
| dev.off() | |
| d1 = d1[order(mads,decreasing=T)[1:2000],] | |
| ho1 = ho1[order(mads,decreasing=T)[1:2000],] | |
| dim(d1) | |
| dim(ho1) | |
| # Plot histograms of MAD after feature selection | |
| pdf('histogram_median_absolute_deviations_selected_genes.pdf') | |
| hist(mads[order(mads,decreasing=T)[1:2000]],xlab='Median Absolute Deviation (MAD)',xlim=c(0,7)) | |
| dev.off() | |
| ##################################### | |
| ### Independent Validation Cohort ### | |
| ##################################### | |
| library(Biobase) | |
| library(GEOquery) | |
| library(limma) | |
| # Load expression data for validation set from GEO | |
| gset <- getGEO("GSE4573", GSEMatrix =TRUE)[[1]] | |
| v1 = exprs(gset) | |
| tmp = gset@featureData@data[,'Gene Symbol'] | |
| names(tmp) = gset@featureData@data[,'ID'] | |
| rownames(v1) = tmp[rownames(v1)] | |
| dim(v1) | |
| # Identify matching features (genes) that will be used to stratify in TCGA | |
| stratGenes = intersect(rownames(d1),unique(rownames(v1))) | |
| length(stratGenes) | |
| # Subset GSE4573 to only those features (genes) selected to stratify the TCGA | |
| v1 = v1[stratGenes,] | |
| dim(v1) | |
| pdf('boxplot_validation_cohort.pdf') | |
| boxplot(log10(v1)) | |
| dev.off() | |
| ### Working with clinical information ### | |
| ## Read in TCGA cinical information | |
| clin = read.csv('data/lusc_clinical.csv',header=T,row.names=1) | |
| clin = clin[gsub('\\.','-',colnames(d1)),] # Include only those with transcriptome profiles | |
| dim(clin) | |
| colnames(clin) | |
| # Summary of all columns for TCGA | |
| summary(clin) | |
| # Sex bias in TCGA | |
| table(clin[,'gender']) | |
| ## Read in GSE4573 cinical information | |
| clin_v = read.csv('data/gse4573_clinical.csv',header=T,row.names=1) | |
| clin_v = clin_v[colnames(v1),] # Include only those with transcriptome profiles | |
| dim(clin_v) | |
| colnames(clin_v) | |
| # Summary of all columns for GSE4573 | |
| summary(clin_v) | |
| # Sex bias in GSE4573 | |
| table(clin_v[,'SEX']) | |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ########################################################## | |
| ## Systems Biology of Disease: LUSC Consenus Clustering ## | |
| ## ______ ______ __ __ ## | |
| ## /\ __ \ /\ ___\ /\ \/\ \ ## | |
| ## \ \ __ \ \ \___ \ \ \ \_\ \ ## | |
| ## \ \_\ \_\ \/\_____\ \ \_____\ ## | |
| ## \/_/\/_/ \/_____/ \/_____/ ## | |
| ## @Developed by: Plaisier Lab ## | |
| ## (https://plaisierlab.engineering.asu.edu/) ## | |
| ## Arizona State University ## | |
| ## 242 ISTB1, 550 E Orange St ## | |
| ## Tempe, AZ 85281 ## | |
| ## @Author: Chris Plaisier ## | |
| ## @License: GNU GPLv3 ## | |
| ########################################################## | |
| ############################## | |
| ### Load data for analysis ### | |
| ############################## | |
| # Set working directory | |
| setwd('C:/Users/cplaisie/Dropbox (ASU)/ASU/Courses/BME_598_494_SysBioDisease/8_22_2019/Fall_2019/Week6/LUSC_ConsensusClustering') | |
| # Data from TCGA Lung squamous cell carcinoma: https://tcga-data.nci.nih.gov/tcga/tcgaCancerDetails.jsp?diseaseType=LUSC&diseaseName=Lung%20squamous%20cell%20carcinoma | |
| gexp = t(na.omit(t(read.csv('data/lusc_mRNA.csv',header=T,row.names=1)))) | |
| d1 = gexp[,-1] | |
| # Hold out 175 samples for testing | |
| holdOut = read.table('data/holdOuts.txt')[,1] | |
| ho1 = d1[,holdOut] | |
| d1 = d1[,which(!(colnames(d1)%in%holdOut))] | |
| # Feature selection using Median Absolute Deviation | |
| mads = apply(d1,1,mad) | |
| d1 = d1[order(mads,decreasing=T)[1:2000],] | |
| ho1 = ho1[order(mads,decreasing=T)[1:2000],] | |
| ############################### | |
| ### Unsupervised clustering ### | |
| ############################### | |
| # Median center each genes expression for use in consensus clustering | |
| d1 = sweep(d1,1,apply(d1,1,median,na.rm=T)) | |
| ho1 = sweep(ho1,1,apply(ho1,1,median,na.rm=T)) | |
| # Do consensus clustering | |
| library(ConsensusClusterPlus) | |
| pdf('consensusClustering.pdf') | |
| title='' | |
| results = ConsensusClusterPlus(d1,maxK=9,reps=500,pItem=0.8,pFeature=1, title=title,clusterAlg="hc",distance="pearson",seed=1262118388.71279) | |
| dev.off() | |
| # Calculate pairwise significance of clusters for 4 clusters | |
| library(sigclust) | |
| c1 = cor(as.matrix(d1)) | |
| clusts = results[[4]]$consensusClass[colnames(d1)] | |
| m1.4 = matrix(ncol=4,nrow=4,dimnames=list(1:4,1:4)) | |
| for(i in 1:4) { | |
| for(j in i:4) { | |
| if(!i==j) { | |
| tmp1 = names(which(clusts==i)) | |
| tmp2 = names(which(clusts==j)) | |
| lab1 = c(rep(1,length(tmp1)),rep(2,length(tmp2))) | |
| tmp3 = sigclust(c1[c(tmp1,tmp2),c(tmp1,tmp2)],nsim=100,labflag=1,label=lab1) | |
| m1.4[i,j] = tmp3@pvalnorm | |
| } | |
| } | |
| } | |
| # Calculate pairwise significance of clusters for 5 clusters | |
| clusts = results[[5]]$consensusClass[colnames(d1)] | |
| m1.5 = matrix(ncol=5,nrow=5,dimnames=list(1:5,1:5)) | |
| for(i in 1:5) { | |
| for(j in i:5) { | |
| if(!i==j) { | |
| tmp1 = names(which(clusts==i)) | |
| tmp2 = names(which(clusts==j)) | |
| lab1 = c(rep(1,length(tmp1)),rep(2,length(tmp2))) | |
| tmp3 = sigclust(c1[c(tmp1,tmp2),c(tmp1,tmp2)],nsim=100,labflag=1,label=lab1) | |
| m1.5[i,j] = tmp3@pvalnorm | |
| } | |
| } | |
| } | |
| # Calculate pairwise significance of clusters for 6 clusters | |
| clusts = results[[6]]$consensusClass[colnames(d1)] | |
| m1.6 = matrix(ncol=6,nrow=6,dimnames=list(1:6,1:6)) | |
| for(i in 1:6) { | |
| for(j in i:6) { | |
| if(!i==j) { | |
| tmp1 = names(which(clusts==i)) | |
| tmp2 = names(which(clusts==j)) | |
| if(length(tmp1)!=1 && length(tmp2)!=1) { | |
| lab1 = c(rep(1,length(tmp1)),rep(2,length(tmp2))) | |
| tmp3 = sigclust(c1[c(tmp1,tmp2),c(tmp1,tmp2)],nsim=100,labflag=1,label=lab1) | |
| m1.6[i,j] = tmp3@pvalnorm | |
| } else { | |
| m1.6[i,j] = 1 | |
| } | |
| } | |
| } | |
| } | |
| # Plot p-values of significance for cluster dissimilarity | |
| library(gplots) | |
| pdf('discovery_sigclust.pdf') | |
| par(mfrow=c(2,2)) | |
| image(x=1:4,y=1:4,m1.4,zlim=c(0,1),col=colorpanel(256,'blue','black','yellow'),main='4 clusters',xlab='k',ylab='k') | |
| image(x=1:5,y=1:5,m1.5,zlim=c(0,1),col=colorpanel(256,'blue','black','yellow'),main='5 clusters',xlab='k',ylab='k') | |
| image(x=1:6,y=1:6,m1.6,zlim=c(0,1),col=colorpanel(256,'blue','black','yellow'),main='6 clusters',xlab='k',ylab='k') | |
| image(matrix(data=seq(from=0,to=1,length.out=10),nrow=10,ncol=1,),col=colorpanel(256,'blue','black','yellow'),main='Legend',axes=F,xlab='P-Value') | |
| axis(1) | |
| dev.off() | |
| # Choose 5 clusters and create centroids for clustering | |
| clusts = results[[5]]$consensusClass[colnames(d1)] | |
| centroids = matrix(nrow=nrow(d1),ncol=5,dimnames=list(rownames(d1),c(1:5))) | |
| for(clust in 1:5) { | |
| tmp1 = names(which(clusts==clust)) | |
| centroids[,clust] = apply(d1[,tmp1],1,median) | |
| } | |
| ############################################################################## | |
| ### Using pamr to reduce number of tested genes (OPTIONAL Not Covered) ### | |
| ### Refer to: http://statweb.stanford.edu/~tibs/PAM/Rdist/doc/readme.html ### | |
| ### And: http://www.pnas.org/content/99/10/6567.long (PMID = 12011421) ### | |
| ############################################################################## | |
| # Load prediction analysis of microarrays (PAM) package | |
| library(pamr) | |
| pdf('pamr_lusc.pdf') | |
| # Build data structure | |
| lusc.data = list(x=as.matrix(d1), y=clusts, geneid=as.character(1:nrow(d1)), genenames=rownames(d1)) | |
| # Train classifier based on data and clusters | |
| lusc.train = pamr.train(lusc.data) | |
| new.scales = pamr.adaptthresh(lusc.train) | |
| # Conduct cross-valiation | |
| lusc.cv = pamr.cv(lusc.train, lusc.data) | |
| # Plot classification error versus threshold value (threshold = 2.0) | |
| pamr.plotcv(lusc.cv) | |
| # Compute the confusion matrix for a particular model (threshold=2.0) | |
| pamr.confusion(lusc.cv, threshold=2.0) | |
| # Plot the cross-validated class probabilities by class | |
| pamr.plotcvprob(lusc.cv, lusc.data, threshold=2.0) | |
| # Plot the class centroids | |
| pamr.plotcen(lusc.train, lusc.data, threshold=2.0) | |
| ## Make a gene plot of the most significant genes | |
| pamr.geneplot(lusc.train, lusc.data, threshold=9) | |
| # List the significant genes | |
| pamr.listgenes(lusc.train, lusc.data, threshold=2.0) | |
| dev.off() | |
| ################################################ | |
| ### Test stratificaton using hold-out cohort ### | |
| ################################################ | |
| c_ho1 = cor(cbind(centroids,ho1),method='spearman')[-c(1:5),1:5] | |
| clusts_ho1 = sapply(1:nrow(c_ho1), function(x) { which(c_ho1[x,]==max(c_ho1[x,])) }) | |
| names(clusts_ho1) = rownames(c_ho1) | |
| # Calculate pairwise significance of clusters for 5 clusters | |
| c1 = cor(as.matrix(ho1)) | |
| m1.ho5 = matrix(ncol=5,nrow=5,dimnames=list(1:5,1:5)) | |
| for(i in 1:5) { | |
| for(j in i:5) { | |
| if(!i==j) { | |
| tmp1 = names(which(clusts_ho1==i)) | |
| tmp2 = names(which(clusts_ho1==j)) | |
| lab1 = c(rep(1,length(tmp1)),rep(2,length(tmp2))) | |
| tmp3 = sigclust(c1[c(tmp1,tmp2),c(tmp1,tmp2)],nsim=100,labflag=1,label=lab1) | |
| m1.ho5[i,j] = tmp3@pvalnorm | |
| } | |
| } | |
| } | |
| # Plot significance of separability on hold-out test cohort | |
| pdf('testing_sigclust.pdf') | |
| par(mfrow=c(2,2)) | |
| image(x=1:5,y=1:5,m1.ho5,zlim=c(0,1),col=colorpanel(256,'blue','black','yellow'),main='5 clusters',xlab='k',ylab='k') | |
| image(matrix(data=seq(from=0,to=1,length.out=10),nrow=10,ncol=1,),col=colorpanel(256,'blue','black','yellow'),main='Legend',axes=F,xlab='P-Value') | |
| axis(1) | |
| dev.off() | |
| # Write out the subtype definitions for TCGA discovery and testing cohorts | |
| write.csv(cbind(names(clusts),clusts),row.names=F,'TCGA_clusters.csv') | |
| write.csv(cbind(names(clusts_ho1),clusts_ho1),row.names=F,'hold_out_clusters.csv') | |
| ############################## | |
| ### Independent Validation ### | |
| ############################## | |
| library(Biobase) | |
| library(GEOquery) | |
| library(limma) | |
| # Load expression data for validation set from GEO | |
| gset <- getGEO("GSE4573", GSEMatrix =TRUE)[[1]] | |
| v1 = exprs(gset) | |
| tmp = gset@featureData@data[,'Gene Symbol'] | |
| names(tmp) = gset@featureData@data[,'ID'] | |
| rownames(v1) = tmp[rownames(v1)] | |
| # Identify matching features (genes) that will be used to stratify in TCGA | |
| stratGenes = intersect(rownames(d1),unique(rownames(v1))) | |
| # Mean center each genes expression | |
| v1 = v1[stratGenes,] | |
| v1 = sweep(v1,1,apply(v1,1,median,na.rm=T)) | |
| # Stratify patients into 5 clusters | |
| valCentroids = centroids[stratGenes,] | |
| c_v1 = cor(cbind(valCentroids,v1),method='spearman')[-c(1:5),1:5] | |
| clusts_v = sapply(1:nrow(c_v1), function(x) { which(c_v1[x,]==max(c_v1[x,])) }) | |
| names(clusts_v) = rownames(c_v1) | |
| # Calculate pairwise significance of clusters for 5 clusters | |
| c1 = cor(as.matrix(v1)) | |
| m1.v5= matrix(ncol=5,nrow=5,dimnames=list(1:5,1:5)) | |
| for(i in 1:5) { | |
| for(j in i:5) { | |
| if(!i==j) { | |
| tmp1 = names(which(clusts_v==i)) | |
| tmp2 = names(which(clusts_v==j)) | |
| lab1 = c(rep(1,length(tmp1)),rep(2,length(tmp2))) | |
| tmp3 = sigclust(c1[c(tmp1,tmp2),c(tmp1,tmp2)],nsim=100,labflag=1,label=lab1) | |
| m1.v5[i,j] = tmp3@pvalnorm | |
| } | |
| } | |
| } | |
| # Plot significance of separability in independent validation cohort | |
| pdf('validation_sigclust.pdf') | |
| par(mfrow=c(2,2)) | |
| image(x=1:5,y=1:5,m1.v5,zlim=c(0,1),col=colorpanel(256,'blue','black','yellow'),main='5 clusters',xlab='k',ylab='k') | |
| image(matrix(data=seq(from=0,to=1,length.out=10),nrow=10,ncol=1,),col=colorpanel(256,'blue','black','yellow'),main='Legend',axes=F,xlab='P-Value') | |
| axis(1) | |
| dev.off() | |
| # Survival association in validation cohort | |
| library(survival) | |
| clin_v = read.csv('data/gse4573_clinical.csv',header=T,row.names=1) | |
| sf1_v = survfit(Surv(as.numeric(as.character(unlist(clin_v[names(clusts_v),'Survival_time']))),clin_v[names(clusts_v),'vital_status']=='1') ~ clusts_v) | |
| summary(sf1_v) | |
| sd1_v = survdiff(Surv(as.numeric(as.character(unlist(clin_v[names(clusts_v),'Survival_time']))),clin_v[names(clusts_v),'vital_status']=='1') ~ clusts_v) | |
| sd1_v | |
| # Plot Kaplan-Meier curve for subtypes | |
| pdf('kaplanMeier_LUSC.pdf') | |
| plot(sf1_v,col=c('black','red','green','blue','orange'),main='Lung Squamous Cell Carcinoma: Survival',ylab='Surviving Fraction',xlab='Time, Months',lwd=2) | |
| legend('topright',legend=c('Subtype 1','Subtype 2','Subtype 3','Subtype 4','Subtype 5'),col=c('black','red','green','blue','orange'),lty=c(1,1,1),title='Legend', inset=0.02,lwd=2) | |
| dev.off() | |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ######################################################### | |
| ## Systems Biology of Disease: Incorporating mutations ## | |
| ## ______ ______ __ __ ## | |
| ## /\ __ \ /\ ___\ /\ \/\ \ ## | |
| ## \ \ __ \ \ \___ \ \ \ \_\ \ ## | |
| ## \ \_\ \_\ \/\_____\ \ \_____\ ## | |
| ## \/_/\/_/ \/_____/ \/_____/ ## | |
| ## @Developed by: Plaisier Lab ## | |
| ## (https://plaisierlab.engineering.asu.edu/) ## | |
| ## Arizona State University ## | |
| ## 242 ISTB1, 550 E Orange St ## | |
| ## Tempe, AZ 85281 ## | |
| ## @Author: Chris Plaisier ## | |
| ## @License: GNU GPLv3 ## | |
| ######################################################### | |
| # Set working directory | |
| setwd('C:/Users/cplaisie/Dropbox (ASU)/ASU/Courses/BME_598_494_SysBioDisease/8_22_2019/Fall_2019/Week6/LUSC_ConsensusClustering') | |
| ### Load data for analysis ### | |
| # Data from TCGA Lung squamous cell carcinoma: https://tcga-data.nci.nih.gov/tcga/tcgaCancerDetails.jsp?diseaseType=LUSC&diseaseName=Lung%20squamous%20cell%20carcinoma | |
| mut = t(na.omit(t(read.csv('data/lusc_mutations.csv',header=T,row.names=1)))) | |
| # Read in cluster IDs from Tuesday | |
| clusts = read.csv('TCGA_clusters.csv',row.names=1,header=T) | |
| # Subset data to those with both mutations and expression determined subtypes | |
| n1 = intersect(colnames(mut),rownames(clusts)) | |
| tmp = t(rbind(mut[,n1],clusters=clusts[n1,1])) | |
| # Which mutations are enriched in each subtype? | |
| pv1 = matrix(ncol=ncol(tmp)-1,nrow=5,dimnames=list(1:5,colnames(tmp[,-ncol(tmp)]))) | |
| for(mutation in colnames(pv1)) { | |
| for(i in rownames(pv1)) { | |
| q = sum(as.numeric(tmp[which(tmp[,'clusters']==i),mutation])) | |
| m = sum(as.numeric(tmp[,mutation])) | |
| n = nrow(tmp)-m | |
| k = length(which(tmp[,'clusters']==i)) | |
| pv1[i,mutation] = phyper(q,m,n,k,lower.tail=F) | |
| } | |
| } | |
| pv1.adjusted = matrix(data=p.adjust(pv1,method='BH'),ncol=ncol(tmp)-1,nrow=5,dimnames=dimnames(pv1)) | |
| # List mutations enriched in each sub-type | |
| enrichedMutations = sapply(1:nrow(pv1), function(x) { names(which(pv1[x,]<=0.05)) }) | |
| enrichedMutations | |
| # List mutations enriched in each sub-type with Benjamini-Hochberg correction | |
| enrichedMutations.adjusted = sapply(1:nrow(pv1.adjusted), function(x) { names(which(pv1.adjusted[x,]<=0.05)) }) | |
| enrichedMutations.adjusted | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment