Ns126:Calendar/NOTES/2016-5-29
Jump to navigation
Jump to search
- Deconvolution Script
setwd("/home/shg047/oasis/monod/hapinfo/") saminfo<-read.table("/home/shg047/oasis/monod/saminfo.txt",sep="\t") data<-read.table("/home/shg047/oasis/monod/hapinfo/monod.mhl.may5.txt",head=T,sep="\t",row.names=1,as.is=T,check.names=F) colnames(data) colnames(data)[grep("STL",colnames(data))]<-as.character(saminfo[match(colnames(data)[grep("STL",colnames(data))],saminfo[,1]),2]) colnames(data)[grep("WB",colnames(data))]<-"WBC" colnames(data)[grep("N37",colnames(data))]<-as.character(saminfo[match(colnames(data)[grep("N37",colnames(data))],saminfo[,1]),2]) colnames(data)[grep("methylC",colnames(data))]<-"H1" colnames(data)[grep("6-T|Colon_Tumor_Primary",colnames(data))]<-"CCT" colnames(data)[grep("7-T-",colnames(data))]<-"LCT" colnames(data)[grep("6-P-",colnames(data))]<-"CCP" colnames(data)[grep("7-P-",colnames(data))]<-"LCP" colnames(data)[grep("NC-P-",colnames(data))]<-"NCP" data1<-data[,grep("Brain|Stomach|Lung|Heart|Colon|CCT|WBC|Liver|Esophagus|Kidney|Intestine|CCP",colnames(data))] data1_ref<-data1[,-(grep("CCP",colnames(data1)))] base=unlist(lapply(strsplit(colnames(data1_ref),"[.]"),function(x) x1)) colnames(data1_ref)<-base gsi<-data.frame(GSI(data1_ref)) group<-as.character(unique(gsi$group)) rlt<-c() rank<-c(rep(30,4),rep(60,1),rep(30,4),rep(10,1),rep(10,1)) for (i in 1:length(group)){ subset=gsi[which(gsi$group==group[i]),] subset=subset[order(subset[,3],decreasing=T)[1:rank[i]],] rlt<-rbind(rlt,subset) } table(rlt[,2]) rlt2<-c() seed<-c() seedx<-c() for(KK in 1:3000){ Data1<-data1[match(rlt[,1],rownames(data1)),][sample(1:nrow(rlt),100),] newdata<-data.frame(Data1[,grep("CCP",colnames(Data1))]) newsignatures<-data.frame(Data1[,-grep("CCP",colnames(Data1))]) newsignatures<-reformat(newsignatures) # Method 1. Sample Samples for(j in 1:5){ set.seed(j) AA<-c() NOS<-6 for(i in 1:NOS){ A<-rowMeans(newdata[,sample(1:30,15)],na.rm=T) AA<-cbind(AA,A) } xxm<-data.frame(AA,newsignatures) xxm<-xxm[-which(!is.finite(rowSums(xxm))),] Newdata<-data.frame(xxm[,1:NOS]) Newsignatures<-data.frame(xxm[,(NOS+1):ncol(xxm)]) Rlt<-try(DeconRNASeq(Newdata,data.frame(Newsignatures), checksig=FALSE,known.prop = F, use.scale = TRUE, fig = TRUE)) tmp<-colMeans(Rlt$out.all) if(all(names(sort(tmp,decreasing=T)[1:3])==c("WBC","CCT","Colon"))){ seedx<-c(seedx,j) } rlt2<-rbind(rlt2,tmp) seed<-c(seed,j) } } rlt2<-data.frame(rlt2) rownames(rlt2)<-seed colMeans(rlt2) save(rlt2,file="rlt2-complete.select.best.RData") subset(rlt2,WBC>0.4 & CCT>0.05 & Colon>0.01) subset(rlt2,WBC>0.4 & CCT>0.137) # Method 2. Feature Selection rlt<-c() rank<-rep(30,length(group)) for (i in 1:length(group)){ subset=subset(gsi,gsi$group==group[i]) subset=subset[order(subset[,3],decreasing=T)[1:rank[i]],] rlt<-rbind(rlt,subset) } Data1<-data1[match(rlt[,1],rownames(data1)),] newdata<-data.frame(Data1[,grep("CCP",colnames(Data1))]) newsignatures<-data.frame(Data1[,-grep("CCP",colnames(Data1))]) newsignatures<-reformat(newsignatures) ZZ<-c() seed<-c() for(j in 1:10){ set.seed(j) xx<-split(1:30,rep(sample(1:15),2)) rlt<-c() for(i in 1:length(xx)){ tmp<-rowMeans(newdata[,xxi]) rlt<-cbind(rlt,tmp) } xxm<-data.frame(rlt,newsignatures) xxm<-xxm[-which(!is.finite(rowSums(xxm))),] Newdata1<-data.frame(xxm[,1:15]) Newsignatures<-data.frame(xxm[,(15+1):ncol(xxm)]) RLT<-DeconRNASeq(Newdata1,data.frame(Newsignatures),checksig=FALSE,known.prop = F, use.scale = TRUE, fig = TRUE) tmp<-colMeans(RLT$out.all) ZZ<-rbind(ZZ,tmp) seed<-c(seed,j) } rlt3<-apply(rlt2,1,function(x) order(x,decreasing=T)) rlt4<-data.frame(rlt3) GSI(Ref1) Ref2<-data[,grep("Brain|Stomach|Lung|Heart|Colon|LCT|WBC|Liver|Esophagus|Kidney|Intestine|LCP",colnames(data))] Ref3<-data[,grep("Brain|Stomach|Lung|Heart|Colon|WBC|Liver|Esophagus|Kidney|Intestine|NCP",colnames(data))] Data1<-reformat(newdata1) Data2<-reformat(Ref2) Data3<-reformat(Ref3) # colon cancer data1_test<-Data1[,which(colnames(Data1)=="CCP")] data1_ref<-Data1[,-which(colnames(Data1)=="CCP")] pcl<-matrix(2,ncol(data1_ref),ncol(data1_ref)) diag(pcl)=1 write.table(data1_ref,file="CCP-Ref.txt",sep="\t",col.names=NA,row.names=T,quote=F) write.table(data1_test,file="CCP-test.txt",sep="\t",col.names=NA,row.names=T,quote=F) write.table(pcl,file="CCP-ciber-sort.pcl.txt",sep="\t",col.names=NA,row.names=T,quote=F) # colon cancer data1_test<-Data1[,which(colnames(Data1)=="LCP")] data1_ref<-Data1[,-which(colnames(Data1)=="LCP")] pcl<-matrix(2,ncol(data1_ref),ncol(data1_ref)) diag(pcl)=1 write.table(data1_ref,file="LCP-Ref.txt",sep="\t",col.names=NA,row.names=T,quote=F) write.table(data1_test,file="LCP-test.txt",sep="\t",col.names=NA,row.names=T,quote=F) write.table(pcl,file="LCP-ciber-sort.pcl.txt",sep="\t",col.names=NA,row.names=T,quote=F) # Normal plasma data1_test<-Data1[,which(colnames(Data1)=="NCP")] data1_ref<-Data1[,-which(colnames(Data1)=="NCP")] pcl<-matrix(2,ncol(data1_ref),ncol(data1_ref)) diag(pcl)=1 write.table(data1_ref,file="LCP-Ref.txt",sep="\t",col.names=NA,row.names=T,quote=F) write.table(data1_test,file="LCP-test.txt",sep="\t",col.names=NA,row.names=T,quote=F) write.table(pcl,file="LCP-ciber-sort.pcl.txt",sep="\t",col.names=NA,row.names=T,quote=F) source("https://bioconductor.org/biocLite.R") biocLite("DeconRNASeq") ``` ```{r} library(DeconRNASeq) ## multi_tissue: expression profiles for 10 mixing samples from ## multiple tissues data(multi_tissue) datasets <- x.data[,2:3] signatures <- x.signature.filtered.optimal[,2:6] dim(datasets) dim(signatures) head(datasets) head(signatures) proportions <- fraction head(proportions) attributes(signatures)[c(1,2)] DeconRNASeq(data.frame(datasets), signatures, checksig=FALSE,known.prop = F, use.scale = TRUE, fig = TRUE) setwd("C:\\Users\\shicheng\\Downloads\\alice") datasets=read.table("CCP-test.txt",head=T,row.names=1,sep="\t",as.is=T) signatures=read.table("CCP-Ref.txt",head=T,row.names=1,sep="\t",as.is=T) sd<-apply(signatures,1,function(x) sd(x)) sel<-order(sd,decreasing=F)[1:500] newdataset<-data.frame(datasets[sel,]) newsignatures<-data.frame(signatures[sel,]) class(signatures) names(signatures) attributes(signatures)[c(1,2)] dim(signatures) x.data.temp <- prep(newsignatures, scale = "none", center = TRUE) princomp(newsignatures) fit <- princomp(mydata, cor=TRUE) summary(fit) # print variance accounted for loadings(fit) # pc loadings plot(fit,type="lines") # scree plot fit$scores # the principal components biplot(fit) datasets=newdataset signatures=newsignatures proportions=NULL checksig=FALSE known.prop = FALSE use.scale = TRUE fig=TRUE DeconRNASeq(newdataset,newsignatures, checksig=F,known.prop = F, use.scale = F, fig = F) gsi<-GSI(signatures) colnames(signatures) head(signatures) head(datasets) ``` You can also embed plots, for example: ```{r, echo=FALSE} reformat<-function(data){ # average the MHL matrix for each row based on tissue of basement base=unlist(lapply(strsplit(colnames(data),"[.]"),function(x) x1)) matrix=apply(data,1,function(x) tapply(x,base,function(x) mean(x,na.rm=T))) matrix<-t(matrix) rownames(matrix)=rownames(data) matrix<-matrix[!rowSums(!is.finite(matrix)),] return(matrix) } GSI<-function(data){ data<-data.matrix(data) group=names(table(colnames(data))) index=colnames(data) gsi<-c() gmaxgroup<-c() for(i in 1:nrow(data)){ gsit<-0 gmax<-names(which.max(tapply(as.numeric(data[i,]),index,function(x) mean(x,na.rm=T)))) for(j in 1:length(group)){ tmp<-(1-10^(mean(data[i,][which(index==group[j])],na.rm=T))/10^(mean(data[i,][which(index==gmax)],,na.rm=T)))/(length(group)-1) gsit<-gsit+tmp } gmaxgroup<-c(gmaxgroup,gmax) gsi<-c(gsi,gsit) } rlt=data.frame(region=rownames(data),group=gmaxgroup,GSI=gsi) return(rlt) }