comparison report_clonality/RScript.r @ 64:55f18bf19d72 draft

Uploaded

63:bd6fb6c03948

sample.count$perc_prod = round(sample.count$Productive / sample.count$All * 100)
sample.count = merge(sample.count, prod.unique.sample.count, by="Sample", all.x=T)
sample.count$perc_prod_un = round(sample.count$Productive_unique / sample.count$All * 100)

sample.count = merge(sample.count , unprod.sample.count, by="Sample", all.x=T)
print(sample.count)
sample.count$perc_unprod = round(sample.count$Unproductive / sample.count$All * 100)
sample.count = merge(sample.count, unprod.unique.sample.count, by="Sample", all.x=T)
sample.count$perc_unprod_un = round(sample.count$Unproductive_unique / sample.count$All * 100)

print(sample.count)

#then sample/replicate
rep.count = merge(input.rep.count, prod.rep.count, by=c("Sample", "Replicate"), all.x=T)
rep.count$perc_prod = round(rep.count$Productive / rep.count$All * 100)
rep.count = merge(rep.count, prod.unique.rep.count, by=c("Sample", "Replicate"), all.x=T)

# ---------------------- calculating the clonality score ----------------------

if("Replicate" %in% colnames(inputdata)) #can only calculate clonality score when replicate information is available
{
  print("Report Clonality - Clonality")
  if(clonality_method == "boyd"){
    samples = split(clonalityFrame, clonalityFrame$Sample, drop=T)
   
    for (sample in samples){
      res = data.frame(paste=character(0))

      coincidence.table = data.frame(table(res$type))
      colnames(coincidence.table) = c("Coincidence Type",  "Raw Coincidence Freq")
      write.table(coincidence.table, file=paste("lymphclon_coincidences_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T)
    }
  } else {
    write.table(clonalityFrame, "clonalityComplete.csv", sep=",",quote=F,row.names=F,col.names=T)
      
    clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")])
    clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
    clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
    clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
    clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
    

	  print(paste("Removing", sum(fltr), "sequences without a identified D"))
	  PRODF = PRODF[!fltr,]
  }
  
  
  print(names(PRODF))
  #ensure certain columns are in the data (files generated with older versions of IMGT Loader)
  col.checks = c("N3.REGION.nt.nb", "N4.REGION.nt.nb")
  for(col.check in col.checks){
	  if(!(col.check %in% names(PRODF))){
		  print(paste(col.check, "not found adding new column"))
		  if(nrow(PRODF) > 0){ #because R is anoying...
			PRODF[,col.check] = 0

		  } else {
			UNPROD = cbind(UNPROD, data.frame(N3.REGION.nt.nb=numeric(0), N4.REGION.nt.nb=numeric(0)))
		  }
	  }
  }
  
  print(names(PRODF))
  
  num_median = function(x, na.rm=T) { as.numeric(median(x, na.rm=na.rm)) }
  
  newData = data.frame(data.table(PRODF)[,list(unique=.N, 
                                               VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),

64:55f18bf19d72

sample.count$perc_prod = round(sample.count$Productive / sample.count$All * 100)
sample.count = merge(sample.count, prod.unique.sample.count, by="Sample", all.x=T)
sample.count$perc_prod_un = round(sample.count$Productive_unique / sample.count$All * 100)

sample.count = merge(sample.count , unprod.sample.count, by="Sample", all.x=T)
sample.count$perc_unprod = round(sample.count$Unproductive / sample.count$All * 100)
sample.count = merge(sample.count, unprod.unique.sample.count, by="Sample", all.x=T)
sample.count$perc_unprod_un = round(sample.count$Unproductive_unique / sample.count$All * 100)

#then sample/replicate
rep.count = merge(input.rep.count, prod.rep.count, by=c("Sample", "Replicate"), all.x=T)
rep.count$perc_prod = round(rep.count$Productive / rep.count$All * 100)
rep.count = merge(rep.count, prod.unique.rep.count, by=c("Sample", "Replicate"), all.x=T)

# ---------------------- calculating the clonality score ----------------------

if("Replicate" %in% colnames(inputdata)) #can only calculate clonality score when replicate information is available
{
  print("Report Clonality - Clonality")
  write.table(clonalityFrame, "clonalityComplete.csv", sep=",",quote=F,row.names=F,col.names=T)
  if(clonality_method == "boyd"){
    samples = split(clonalityFrame, clonalityFrame$Sample, drop=T)
   
    for (sample in samples){
      res = data.frame(paste=character(0))

      coincidence.table = data.frame(table(res$type))
      colnames(coincidence.table) = c("Coincidence Type",  "Raw Coincidence Freq")
      write.table(coincidence.table, file=paste("lymphclon_coincidences_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T)
    }
  } else {
    clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")])
    
    #write files for every coincidence group of >1
    samples = unique(clonalFreq$Sample)
    for(sample in samples){
		clonalFreqSample = clonalFreq[clonalFreq$Sample == sample,]
		if(max(clonalFreqSample$Type) > 1){
			for(i in 2:max(clonalFreqSample$Type)){
				clonalFreqSample = clonalFreqSample[clonalFreqSample$Type == i,]
				print(clonalFreqSample[clonalFreqSample$Type == i,])
				PRODF.sub = PRODF[PRODF$clonaltype %in% clonalFreqSample$clonaltype,]
				PRODF.sub = PRODF.sub[order(PRODF.sub$clonaltype),]
				write.table(PRODF.sub, file=paste("coincidences_", sample, "_", i, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
			}
		}
	}
	
    clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
    clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
    clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
    clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
    

	  print(paste("Removing", sum(fltr), "sequences without a identified D"))
	  PRODF = PRODF[!fltr,]
  }
  
  
  #ensure certain columns are in the data (files generated with older versions of IMGT Loader)
  col.checks = c("N.REGION.nt.nb", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "N3.REGION.nt.nb", "N4.REGION.nt.nb")
  for(col.check in col.checks){
	  if(!(col.check %in% names(PRODF))){
		  print(paste(col.check, "not found adding new column"))
		  if(nrow(PRODF) > 0){ #because R is anoying...
			PRODF[,col.check] = 0

		  } else {
			UNPROD = cbind(UNPROD, data.frame(N3.REGION.nt.nb=numeric(0), N4.REGION.nt.nb=numeric(0)))
		  }
	  }
  }
  
  num_median = function(x, na.rm=T) { as.numeric(median(x, na.rm=na.rm)) }
  
  newData = data.frame(data.table(PRODF)[,list(unique=.N, 
                                               VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),