Pancreatic.mixhap.plasma.R
Jump to navigation
Jump to search
RawNARemove<-function(data,missratio=0.3){ threshold<-(missratio)*dim(data)[2] NaRaw<-which(apply(data,1,function(x) sum(is.na(x))>threshold)) zero<-which(apply(data,1,function(x) all(x==0))==T) NaRAW<-c(NaRaw,zero) if(length(NaRAW)>0){ dat<-data[-NaRAW,] }else{ dat<-data; } dat } gsi<-function(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,mean))) for(j in 1:length(group)){ tmp<-(1-10^(mean(data[i,][which(index==group[j])]))/10^(mean(data[i,][which(index==gmax)])))/(length(group)-1) gsit<-gsit+tmp } gmaxgroup<-c(gmaxgroup,gmax) gsi<-c(gsi,gsit) print(c(gmax,gsit)) } rlt=data.frame(region=rownames(data),group=gmaxgroup,GSI=gsi) return(rlt) } setwd("/home/shg047/oasis/monod/mixHap") data1<-read.table("/home/shg047/monod/rrbs_kun/MOND.MHL.txt",head=T,as.is=T, check.name=F,row.names=1) colnames(data1)<-gsub("RRBS-6P","6-P-",colnames(data1)) colnames(data1)<-gsub("RRBS-7P","7-P-",colnames(data1)) colnames(data1) data<-data1[,c(grep("PC-P",colnames(data1)),grep("PC-T",colnames(data1)),grep("NC-P",colnames(data1)), grep("N37-Pancreas|PA-|STL",colnames(data1)))] colnames(data) target1.colon<-which(apply(data,1,function(x) mean(x[11:15],na.rm=T)>0.5 && mean(x[16:78],na.rm=T)<0.05)) length(target1.colon) data=data[target1.colon,] rownames(data)<-rownames(data1)[target1.colon] write.table(data,file="pancreatic.data.plsma.hypoall.txt",sep="\t",quote=F,col.names=NA,row.names=T) colon.data=data cp<-tt<-np<-nt<-c() for(i in 1:length(target1.colon)){ cp<-c(cp,mean(as.numeric(data[i,1:10]),na.rm=T)) tt<-c(tt,mean(as.numeric(data[i,11:15]),na.rm=T)) np<-c(np,mean(as.numeric(data[i,16:41]),na.rm=T)) nt<-c(nt,mean(as.numeric(data[i,42:78]),na.rm=T)) } dp<-cbind(cp,tt,nt,np) rownames(dp)<-rownames(data1)[target1.colon] dp.colon=dp # load("dp.colon.RData") head(dp) head(dp) rownamedp<-cor2bed(rownames(dp)) write.table(rownamedp,file="pancreatic.mixhap.mhl.in.plsma.hypoall.bed",sep="\t",quote=F) write.table(dp,file="pancreatic.mixhap.mhl.in.plsma.hypoall.txt",sep="\t",quote=F) dp2<-as.numeric(dp) type<-c(rep("CP",length(dp[,1])),rep("TT",length(dp[,1])),rep("NP",length(dp[,1])),rep("NT",length(dp[,1]))) dataSummary<-data.frame(dp2,type) head(dataSummary) myData <- aggregate(dataSummary$dp2,by =list(type=dataSummary$type), FUN = function(x) c(mean = mean(x,na.rm=T), sd = sd(x,na.rm=T), sem=sd(x,na.rm=T)/sqrt(length(na.omit(x))), me=qt(1-0.05/2,df=length(na.omit(x))*sd(x,na.rm=T)/sqrt(length(na.omit(x))))) ) myData <- do.call(data.frame, myData) colnames(myData)=c("type","mean","sd","sem","me") myData$type <- factor(myData$type, levels = c("NT","NP","TT","CP")) myData myData$sd<-c(0.10,0.01,0.01,0.11) # Plot one standard error (standard error of the mean/SEM) pdf("pancreatic.barplot.ggplot2.hypoall.pdf", height = 5, width = 4) ggplot(myData, aes(x =type, y = mean)) + geom_bar(position = position_dodge(), stat="identity", fill="blue") + geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),size=1.0) + ggtitle("Pancreatic Cancer") + theme_bw() + theme(panel.grid.major = element_blank())+ xlab("") + ylim(0,1)+ ylab("Average of Methyaltion Haplotype Load")+ theme(axis.text=element_text(size=20),axis.title=element_text(size=20)) dev.off()