comparison report_clonality/RScript.r @ 44:1d8728f3ff37 draft

Uploaded

43:2325074a8461

      ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
      xlab("D genes") + 
      ylab("V Genes") +
      theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), panel.grid.major = element_line(colour = "gainsboro"))
    
    png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name)))
    print(img)
    dev.off()
    
    ggsave(paste("HeatmapVD_", unique(dat[3])[1,1] , ".pdf", sep=""), img, height=13, width=8)
    

  maxVD = data.frame(data.table(VandDCount)[, list(max=max(l)), by=c("Sample")])
  VandDCount = merge(VandDCount, maxVD, by.x="Sample", by.y="Sample", all.x=T)
  VandDCount$relLength = VandDCount$l / VandDCount$max
  check = is.nan(VandDCount$relLength)
  if(any(check)){
	VandDCount[check,"relLength"] = 0
  }
  
  cartegianProductVD = expand.grid(Top.V.Gene = Vchain$v.name, Top.D.Gene = Dchain$v.name)
  
  completeVD = merge(VandDCount, cartegianProductVD, by.x=c("Top.V.Gene", "Top.D.Gene"), by.y=c("Top.V.Gene", "Top.D.Gene"), all=TRUE)
 
  completeVD = merge(completeVD, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE)
 
  completeVD = merge(completeVD, Dchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE)
  
  fltr = is.nan(completeVD$relLength)
  if(all(fltr)){
	  completeVD[fltr,"relLength"] = 0
  }

    ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
    xlab("J genes") + 
    ylab("V Genes") +
    theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), panel.grid.major = element_line(colour = "gainsboro"))
  
  png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name)))
  print(img)
  dev.off()
  
  ggsave(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".pdf", sep=""), img, height=11, width=4)
  
  write.table(x=acast(dat, Top.V.Gene~Top.J.Gene, value.var="Length"), file=paste("HeatmapVJ_", unique(dat[3])[1,1], ".txt", sep=""), sep="\t",quote=F,row.names=T,col.names=NA)
}

VandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.J.Gene", "Sample")])

VandJCount$l = log(VandJCount$Length)
maxVJ = data.frame(data.table(VandJCount)[, list(max=max(l)), by=c("Sample")])

check = is.nan(VandJCount$relLength)
if(any(check)){
	VandJCount[check,"relLength"] = 0
}

cartegianProductVJ = expand.grid(Top.V.Gene = Vchain$v.name, Top.J.Gene = Jchain$v.name)

completeVJ = merge(VandJCount, cartegianProductVJ, all.y=TRUE)
completeVJ = merge(completeVJ, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE)
completeVJ = merge(completeVJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE)

fltr = is.nan(completeVJ$relLength)
if(any(fltr)){
	completeVJ[fltr,"relLength"] = 1
}

      ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
      xlab("J genes") + 
      ylab("D Genes") +
      theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), panel.grid.major = element_line(colour = "gainsboro"))
    
    png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name)))
    print(img)
    dev.off()
    
    ggsave(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".pdf", sep=""), img, width=4, height=7)
    

    DandJCount[check,"relLength"] = 0
  }
  
  cartegianProductDJ = expand.grid(Top.D.Gene = Dchain$v.name, Top.J.Gene = Jchain$v.name)
  
  completeDJ = merge(DandJCount, cartegianProductDJ, all.y=TRUE)
  completeDJ = merge(completeDJ, revDchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE)
  completeDJ = merge(completeDJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE)
  
  fltr = is.nan(completeDJ$relLength)
  if(any(fltr)){
	  completeDJ[fltr, "relLength"] = 1
  }

	clonaltype.in.replicates = inputdata
	clonaltype.in.replicates = clonaltype.in.replicates[clonaltype.in.replicates$Functionality %in% c("productive (see comment)","productive"),]
	clonaltype.in.replicates = clonaltype.in.replicates[!(is.na(clonaltype.in.replicates$ID) | is.na(clonaltype.in.replicates$Top.V.Gene) | is.na(clonaltype.in.replicates$Top.J.Gene)),]
	clonaltype = unlist(strsplit(clonaltype, ","))

	clonaltype.in.replicates$clonaltype = do.call(paste, c(clonaltype.in.replicates[paste_columns], sep = ":"))

	clonaltype.in.replicates = clonaltype.in.replicates[!duplicated(clonaltype.in.replicates$clonaltype),]

	clonaltype = clonaltype[-which(clonaltype == "Sample")]

44:1d8728f3ff37

      ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
      xlab("D genes") + 
      ylab("V Genes") +
      theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), panel.grid.major = element_line(colour = "gainsboro"))
    
    png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=200+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name)))
    print(img)
    dev.off()
    
    ggsave(paste("HeatmapVD_", unique(dat[3])[1,1] , ".pdf", sep=""), img, height=13, width=8)
    

  maxVD = data.frame(data.table(VandDCount)[, list(max=max(l)), by=c("Sample")])
  VandDCount = merge(VandDCount, maxVD, by.x="Sample", by.y="Sample", all.x=T)
  VandDCount$relLength = VandDCount$l / VandDCount$max
  check = is.nan(VandDCount$relLength)
  if(any(check)){
	  VandDCount[check,"relLength"] = 0
  }

  completeVD = merge(VandDCount, revVchain, by.x="Top.V.Gene", by.y="v.name", all=TRUE)
  completeVD = merge(completeVD, Dchain, by.x="Top.D.Gene", by.y="v.name", all=TRUE)
  
  fltr = is.nan(completeVD$relLength)
  if(all(fltr)){
	  completeVD[fltr,"relLength"] = 0
  }

    ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
    xlab("J genes") + 
    ylab("V Genes") +
    theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), panel.grid.major = element_line(colour = "gainsboro"))
  
  png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=200+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name)))
  print(img)
  dev.off()
  
  ggsave(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".pdf", sep=""), img, height=11, width=4)
  
  write.table(x=acast(dat, Top.V.Gene~Top.J.Gene, value.var="Length"), file=paste("HeatmapVJ_", unique(dat[3])[1,1], ".txt", sep=""), sep="\t",quote=F,row.names=T,col.names=NA)
}



VandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.J.Gene", "Sample")])

VandJCount$l = log(VandJCount$Length)
maxVJ = data.frame(data.table(VandJCount)[, list(max=max(l)), by=c("Sample")])

check = is.nan(VandJCount$relLength)
if(any(check)){
	VandJCount[check,"relLength"] = 0
}

completeVJ = merge(VandJCount, revVchain, by.x="Top.V.Gene", by.y="v.name", all=TRUE)
completeVJ = merge(completeVJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all=TRUE)

fltr = is.nan(completeVJ$relLength)
if(any(fltr)){
	completeVJ[fltr,"relLength"] = 1
}

      ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
      xlab("J genes") + 
      ylab("D Genes") +
      theme(panel.background = element_rect(fill = "white", colour="black"),text = element_text(size=15, colour="black"), panel.grid.major = element_line(colour = "gainsboro"))
    
    png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=200+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name)))
    print(img)
    dev.off()
    
    ggsave(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".pdf", sep=""), img, width=4, height=7)
    

    DandJCount[check,"relLength"] = 0
  }
  
  cartegianProductDJ = expand.grid(Top.D.Gene = Dchain$v.name, Top.J.Gene = Jchain$v.name)
  
  completeDJ = merge(DandJCount, revDchain, by.x="Top.D.Gene", by.y="v.name", all=TRUE)
  completeDJ = merge(completeDJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all=TRUE)
  
  fltr = is.nan(completeDJ$relLength)
  if(any(fltr)){
	  completeDJ[fltr, "relLength"] = 1
  }

	clonaltype.in.replicates = inputdata
	clonaltype.in.replicates = clonaltype.in.replicates[clonaltype.in.replicates$Functionality %in% c("productive (see comment)","productive"),]
	clonaltype.in.replicates = clonaltype.in.replicates[!(is.na(clonaltype.in.replicates$ID) | is.na(clonaltype.in.replicates$Top.V.Gene) | is.na(clonaltype.in.replicates$Top.J.Gene)),]
	clonaltype = unlist(strsplit(clonaltype, ","))

	clonaltype.in.replicates$clonaltype = do.call(paste, c(clonaltype.in.replicates[clonaltype], sep = ":"))

	clonaltype.in.replicates = clonaltype.in.replicates[!duplicated(clonaltype.in.replicates$clonaltype),]

	clonaltype = clonaltype[-which(clonaltype == "Sample")]