Mercurial > repos > davidvanzessen > complete_immunerepertoire_igg
changeset 10:b8db36cfe6ad draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Mon, 12 Jan 2015 11:07:58 -0500 |
| parents | 7dbc9ebcefce |
| children | 1b356cbde361 |
| files | RScript.r |
| diffstat | 1 files changed, 55 insertions(+), 53 deletions(-) [+] |
line wrap: on
line diff
--- a/RScript.r Mon Jan 12 05:35:42 2015 -0500 +++ b/RScript.r Mon Jan 12 11:07:58 2015 -0500 @@ -31,6 +31,7 @@ outfile = args[2] #path to output file outdir = args[3] #path to output folder (html/images/data) clonaltype = args[4] #clonaltype definition, or 'none' for no unique filtering +ct = unlist(strsplit(clonaltype, ",")) species = args[5] #human or mouse locus = args[6] # IGH, IGK, IGL, TRB, TRA, TRG or TRD filterproductive = ifelse(args[7] == "yes", T, F) #should unproductive sequences be filtered out? (yes/no) @@ -42,13 +43,15 @@ setwd(outdir) # remove weird rows -inputdata = inputdata[inputdata$Sample != "",] +inputdata = inputdata[inputdata$Sample != "",] #remove the allele from the V,D and J genes inputdata$Top.V.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.V.Gene) inputdata$Top.D.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.D.Gene) inputdata$Top.J.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.J.Gene) + inputdata$clonaltype = 1:nrow(inputdata) + PRODF = inputdata UNPROD = inputdata if(filterproductive){ @@ -63,10 +66,9 @@ #remove duplicates based on the clonaltype if(clonaltype != "none"){ + clonaltype = paste(clonaltype, ",Sample", sep="") #add sample column to clonaltype, unique within samples PRODF$clonaltype = do.call(paste, c(PRODF[unlist(strsplit(clonaltype, ","))], sep = ":")) PRODF = PRODF[!duplicated(PRODF$clonaltype), ] - UNPROD$clonaltype = do.call(paste, c(UNPROD[unlist(strsplit(clonaltype, ","))], sep = ":")) - UNPROD = UNPROD[!duplicated(UNPROD$clonaltype), ] } PRODF$freq = 1 @@ -96,9 +98,8 @@ inputdata.dt = data.table(inputdata) #for speed -ct = unlist(strsplit(clonaltype, ",")) if(clonaltype == "none"){ - ct = c("ID") + ct = c("clonaltype") } inputdata.dt$samples_replicates = paste(inputdata.dt$Sample, inputdata.dt$Replicate, sep="_") @@ -471,6 +472,7 @@ { clonalityFrame = inputdata if(clonaltype != "none"){ + clonalityFrame$clonaltype = do.call(paste, c(clonalityFrame[unlist(strsplit(clonaltype, ","))], sep = ":")) clonalityFrame$ReplicateConcat = paste(clonalityFrame$clonaltype, clonalityFrame$Sample, clonalityFrame$Replicate, sep = ":") clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ] } @@ -515,7 +517,7 @@ ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) lapply(ReplicateSplit, FUN=ReplicatePrint) - ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(Reads), ReadsSquaredSum=sum(squared)), by=c("Sample")]) + ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")]) clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T) @@ -558,54 +560,54 @@ if(all(imgtcolumns %in% colnames(inputdata))) { newData = data.frame(data.table(PRODF)[,list(unique=.N, - VH.DEL=mean(X3V.REGION.trimmed.nt.nb, na.rm=T), - P1=mean(P3V.nt.nb, na.rm=T), - N1=mean(N1.REGION.nt.nb, na.rm=T), - P2=mean(P5D.nt.nb, na.rm=T), - DEL.DH=mean(X5D.REGION.trimmed.nt.nb, na.rm=T), - DH.DEL=mean(X3D.REGION.trimmed.nt.nb, na.rm=T), - P3=mean(P3D.nt.nb, na.rm=T), - N2=mean(N2.REGION.nt.nb, na.rm=T), - P4=mean(P5J.nt.nb, na.rm=T), - DEL.JH=mean(X5J.REGION.trimmed.nt.nb, na.rm=T), - Total.Del=( mean(X3V.REGION.trimmed.nt.nb, na.rm=T) + - mean(X5D.REGION.trimmed.nt.nb, na.rm=T) + - mean(X3D.REGION.trimmed.nt.nb, na.rm=T) + - mean(X5J.REGION.trimmed.nt.nb, na.rm=T)), - - Total.N=( mean(N1.REGION.nt.nb, na.rm=T) + - mean(N2.REGION.nt.nb, na.rm=T)), - - Total.P=( mean(P3V.nt.nb, na.rm=T) + - mean(P5D.nt.nb, na.rm=T) + - mean(P3D.nt.nb, na.rm=T) + - mean(P5J.nt.nb, na.rm=T))), - by=c("Sample")]) + VH.DEL=mean(X3V.REGION.trimmed.nt.nb, na.rm=T), + P1=mean(P3V.nt.nb, na.rm=T), + N1=mean(N1.REGION.nt.nb, na.rm=T), + P2=mean(P5D.nt.nb, na.rm=T), + DEL.DH=mean(X5D.REGION.trimmed.nt.nb, na.rm=T), + DH.DEL=mean(X3D.REGION.trimmed.nt.nb, na.rm=T), + P3=mean(P3D.nt.nb, na.rm=T), + N2=mean(N2.REGION.nt.nb, na.rm=T), + P4=mean(P5J.nt.nb, na.rm=T), + DEL.JH=mean(X5J.REGION.trimmed.nt.nb, na.rm=T), + Total.Del=( mean(X3V.REGION.trimmed.nt.nb, na.rm=T) + + mean(X5D.REGION.trimmed.nt.nb, na.rm=T) + + mean(X3D.REGION.trimmed.nt.nb, na.rm=T) + + mean(X5J.REGION.trimmed.nt.nb, na.rm=T)), + + Total.N=( mean(N1.REGION.nt.nb, na.rm=T) + + mean(N2.REGION.nt.nb, na.rm=T)), + + Total.P=( mean(P3V.nt.nb, na.rm=T) + + mean(P5D.nt.nb, na.rm=T) + + mean(P3D.nt.nb, na.rm=T) + + mean(P5J.nt.nb, na.rm=T))), + by=c("Sample")]) write.table(newData, "junctionAnalysisProd.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F) - newData = data.frame(data.table(UNPROD)[,list(unique=.N, - VH.DEL=mean(X3V.REGION.trimmed.nt.nb, na.rm=T), - P1=mean(P3V.nt.nb, na.rm=T), - N1=mean(N1.REGION.nt.nb, na.rm=T), - P2=mean(P5D.nt.nb, na.rm=T), - DEL.DH=mean(X5D.REGION.trimmed.nt.nb, na.rm=T), - DH.DEL=mean(X3D.REGION.trimmed.nt.nb, na.rm=T), - P3=mean(P3D.nt.nb, na.rm=T), - N2=mean(N2.REGION.nt.nb, na.rm=T), - P4=mean(P5J.nt.nb, na.rm=T), - DEL.JH=mean(X5J.REGION.trimmed.nt.nb, na.rm=T), - Total.Del=( mean(X3V.REGION.trimmed.nt.nb, na.rm=T) + - mean(X5D.REGION.trimmed.nt.nb, na.rm=T) + - mean(X3D.REGION.trimmed.nt.nb, na.rm=T) + - mean(X5J.REGION.trimmed.nt.nb, na.rm=T)), - - Total.N=( mean(N1.REGION.nt.nb, na.rm=T) + - mean(N2.REGION.nt.nb, na.rm=T)), - - Total.P=( mean(P3V.nt.nb, na.rm=T) + - mean(P5D.nt.nb, na.rm=T) + - mean(P3D.nt.nb, na.rm=T) + - mean(P5J.nt.nb, na.rm=T))), - by=c("Sample")]) + newData = data.frame(data.table(UNPROD)[,list(unique=.N, + VH.DEL=mean(X3V.REGION.trimmed.nt.nb, na.rm=T), + P1=mean(P3V.nt.nb, na.rm=T), + N1=mean(N1.REGION.nt.nb, na.rm=T), + P2=mean(P5D.nt.nb, na.rm=T), + DEL.DH=mean(X5D.REGION.trimmed.nt.nb, na.rm=T), + DH.DEL=mean(X3D.REGION.trimmed.nt.nb, na.rm=T), + P3=mean(P3D.nt.nb, na.rm=T), + N2=mean(N2.REGION.nt.nb, na.rm=T), + P4=mean(P5J.nt.nb, na.rm=T), + DEL.JH=mean(X5J.REGION.trimmed.nt.nb, na.rm=T), + Total.Del=( mean(X3V.REGION.trimmed.nt.nb, na.rm=T) + + mean(X5D.REGION.trimmed.nt.nb, na.rm=T) + + mean(X3D.REGION.trimmed.nt.nb, na.rm=T) + + mean(X5J.REGION.trimmed.nt.nb, na.rm=T)), + + Total.N=( mean(N1.REGION.nt.nb, na.rm=T) + + mean(N2.REGION.nt.nb, na.rm=T)), + + Total.P=( mean(P3V.nt.nb, na.rm=T) + + mean(P5D.nt.nb, na.rm=T) + + mean(P3D.nt.nb, na.rm=T) + + mean(P5J.nt.nb, na.rm=T))), + by=c("Sample")]) write.table(newData, "junctionAnalysisUnProd.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F) }