diff --git a/.DS_Store b/.DS_Store index 126ff7a..0e2a99c 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/DNAfunction.R b/DNAfunction.R index 8ff487a..a776d12 100755 --- a/DNAfunction.R +++ b/DNAfunction.R @@ -54,6 +54,7 @@ ASCATout<-function(ASCATdata,sample,chromosome){ ACN=data.frame(nMajor=segment$nMajor,nMinor=segment$nMinor) segment$nMajor=apply(ACN,1,max) segment$nMinor=apply(ACN,1,min) + segment=data.frame(chr=segment$chr,startpos=segment$startpos,endpos=segment$endpos,nMajor=segment$nMajor,nMinor=segment$nMinor) alpha=ascat.output$aberrantcellfraction ploidy=ascat.output$ploidy DNAout=list(alpha=alpha,segment=segment,ploidy=ploidy) @@ -647,8 +648,16 @@ getMajorMinorCN <- function(state, symmetric = TRUE){ } return(list(majorCN = majorCN, minorCN = minorCN)) } -TITANout<-function(DNAinput){ +TITANout<-function(DNAinput,chromosome){ + DNAinput=DNAinput[DNAinput$nfrac>=0.2&DNAinput$nfrac<=0.8,] + rlength=nchar(as.character(DNAinput$ref)) + alength=nchar(as.character(DNAinput$alt)) + DNAinput=DNAinput[rlength==1&alength==1,] titandata=data.frame(chr=DNAinput$chr,pos=DNAinput$pos,ref=DNAinput$ref,refCount=DNAinput$refNumT,Nref=DNAinput$alt,NrefCount=DNAinput$altNumT) + chro=intersect(chromosome,names(table(titandata$chr))[table(titandata$chr)>1]) + index=match(as.character(titandata$chr),chro) + titandata=titandata[!is.na(index),] + DNAinput=DNAinput[!is.na(index),] write.table(titandata,"TITAN.input",col.names=TRUE,row.names=FALSE,sep="\t",quote=FALSE) numClusters <- 2 params <- loadDefaultParameters(copyNumber=5,numberClonalClusters=numClusters) @@ -663,10 +672,11 @@ TITANout<-function(DNAinput){ results$AllelicRatio=as.numeric(results$AllelicRatio) results$LogRatio=as.numeric(results$LogRatio) segs <- outputTitanSegments(results, id = "test", convergeParams,filename=NULL) + segment=data.frame(chr=segs$Chromosome,stratpos=segs$Start_Position.bp.,endpos=segs$End_Position.bp.,nMajor=segs$MajorCN,nMinor=segs$MinorCN) ploidy <- tail(convergeParams$phi, 1) normal <- tail(convergeParams$n, 1) - mean(as.numeric(segs$Cellular_Frequency[!is.na(segs$Cellular_Frequency)])*(1-normal)) - DNAout=list(segment=segs,ploidy=ploidy,alpha=1-normal) + #mean(as.numeric(segs$Cellular_Frequency[!is.na(segs$Cellular_Frequency)])*(1-normal)) + DNAout=list(segment=segment,ploidy=ploidy,alpha=1-normal) return(DNAout) } @@ -705,7 +715,7 @@ DNACNA<-function(DNAinput,sample,tempth,chromosome,i){ return(ASCATres) }else if (i==2){ registerDoMC(cores = 4) - TITANres<-try(TITANout(DNAinput=DNAinput),silent=TRUE) + TITANres<-try(TITANout(DNAinput=DNAinput,chromosome=chromosome),silent=TRUE) return (TITANres) }else if (i ==3){ FACETSres<-try(FACETSout(DNAinput=DNAinput),silent=TRUE) @@ -844,9 +854,9 @@ DNArun<-function(SNPinput,somaticinput,sample,temppath){ times=1000 res=foreach(i = 1:4) %dopar% DNACNA(DNAinput=DNAinput,sample=sample,tempth=temppath,chromosome=chromosome,i) print.noquote("Somatic mutation copy number") - somaticres<-foreach(i =1:4) %dopar% somaticCN(res=res,i,somaticdata=somaticdata,method=methods) + somaticres<-foreach(i =1:length(res)) %dopar% somaticCN(res=res,i,somaticdata=somaticdata,method=methods) somaticindex=0 - for (i in 1:4){ + for (i in 1:length(somaticres)){ if (!is.null(somaticres[[i]])){ if (dim(somaticres[[i]])[1]>0){ somaticindex=1 @@ -855,15 +865,15 @@ DNArun<-function(SNPinput,somaticinput,sample,temppath){ } } if (somaticindex!=0){ - randomdata<-foreach(i=1:4) %dopar% randomSample(data=data.frame(Tsum=somaticres[[i]]$Tsum,tfrac=somaticres[[i]]$tfrac),times=times) - randomout<-foreach(i =1:4) %dopar% randomCN(res,i,randomdata,somaticres,method=methods) - aveDIF<-foreach(i=1:4) %dopar% randomDIF(randomout,somaticres,i) - realDIF<-foreach(i=1:4) %dopar% realDiff(somaticres,i) - realsub<-foreach(i=1:4) %dopar% realSub(somaticres,i) - randomsub<-foreach(i = 1:4) %dopar% subcloneP(x=randomout[[i]]) + randomdata<-foreach(i=1:length(somaticres)) %dopar% randomSample(data=data.frame(Tsum=somaticres[[i]]$Tsum,tfrac=somaticres[[i]]$tfrac),times=times) + randomout<-foreach(i =1:length(somaticres)) %dopar% randomCN(res,i,randomdata,somaticres,method=methods) + aveDIF<-foreach(i=1:length(somaticres)) %dopar% randomDIF(randomout,somaticres,i) + realDIF<-foreach(i=1:length(somaticres)) %dopar% realDiff(somaticres,i) + realsub<-foreach(i=1:length(somaticres)) %dopar% realSub(somaticres,i) + randomsub<-foreach(i = 1:length(somaticres)) %dopar% subcloneP(x=randomout[[i]]) PDl=list() PSl=list() - for (i in 1:4){ + for (i in 1:length(realDIF)){ if (!is.null(realDIF[[i]])){ PDl[[i]]=1-length(aveDIF[[i]][aveDIF[[i]]>=realDIF[[i]]])/length(aveDIF[[i]]) PSl[[i]]=length(randomsub[[i]][randomsub[[i]]=segment$start&DNAinput$pos<=segment$end,] + ll=dim(DNAinput)[1] + gdouble=dim(subdata)[1] + }else{ + ll=sum(as.numeric(segment$end-segment$start)) + gdouble=sum(as.numeric(segment$end[segment$TCN>2]-segment$start[segment$TCN>2])) + } + + return(gdouble/ll) + +} DACRE<-function(resout){ germ=resout[resout$type=="Germline",] germ1=do.call(rbind,apply(germ[,8:10],1,altalle)) diff --git a/Texomer.R b/Texomer.R index 3917364..a2799ab 100644 --- a/Texomer.R +++ b/Texomer.R @@ -40,7 +40,12 @@ if (length(args)==6){ heterogeneity=DNAout$heterogeneity DNApurity=DNAout$alpha ploidy=DNAout$ploidy - summaryres=rbind(c("DNApurity",DNApurity),c("Heterogeneity",heterogeneity),c("Ploidy",ploidy)) + segment=DNAout$segment + colnames(segment)=c("chr","start","end","Dmajor","Dminor") + somatic=DNAout$somatic + WGDs=WGD(segment,SNPinput) + summaryres=rbind(c("DNApurity",DNApurity),c("Heterogeneity",heterogeneity),c("Ploidy",ploidy),c("WGD",WGDs)) + print.noquote(paste("The initial method = ",DNAout$method,sep="")) if (RNAinput=="NA"){ segment=DNAout$segment somatic=DNAout$somatic diff --git a/Texomer.py b/Texomer.py index 1ea6c9c..021626f 100755 --- a/Texomer.py +++ b/Texomer.py @@ -270,7 +270,7 @@ def main(): elif options.Input == "BAM": if options.Tumor and options.Normal: if not options.reference: - print "Please input the location of refereence sequence through -f." + print "Please input the location of reference sequence through -f." sys.exit(1) else: ref=options.reference