[Patchwork-commits] r168 - in pkg: . TAPS TAPS/R TAPS/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 22 16:29:39 CEST 2013
Author: sebastian_d
Date: 2013-04-22 16:29:38 +0200 (Mon, 22 Apr 2013)
New Revision: 168
Added:
pkg/TAPS/
pkg/TAPS/.Rhistory
pkg/TAPS/DESCRIPTION
pkg/TAPS/NAMESPACE
pkg/TAPS/R/
pkg/TAPS/R/TAPS.r
pkg/TAPS/R/sysdata.rda
pkg/TAPS/R/zzz.r
pkg/TAPS/man/
pkg/TAPS/man/TAPS_call.Rd
pkg/TAPS/man/TAPS_plot.Rd
pkg/TAPS/man/TAPS_region.Rd
Log:
Initial addition of TAPS to Rforge
Added: pkg/TAPS/.Rhistory
===================================================================
--- pkg/TAPS/.Rhistory (rev 0)
+++ pkg/TAPS/.Rhistory 2013-04-22 14:29:38 UTC (rev 168)
@@ -0,0 +1,512 @@
+dupcor<-duplicateCorrelation(Data, design, block=block)
+dupcor$consensus.correlation
+cont.matrix<- makeContrasts(C_ATRAvsA=C_ATRA-C, IKK1_ATRAvsIKK1=IKK1_ATRA-IKK1, IKK2_ATRAvsIKK2=IKK2_ATRA-IKK2, levels=design)
+fit<-lmFit(Data,design, block=block, correlation=dupcor$consensus.correlation)
+fit2<-contrasts.fit(fit, cont.matrix)
+fit2<-eBayes(fit2)
+summary(decideTests(fit2))
+vennDiagram(fit2)
+results=decideTests(fit2)
+vennDiagram(results)
+cont.matrix<- makeContrasts(C_ATRAvsC=C_ATRA-C, IKK1_ATRAvsIKK1=IKK1_ATRA-IKK1, IKK2_ATRAvsIKK2=IKK2_ATRA-IKK2, levels=design)
+fit<-lmFit(Data,design, block=block, correlation=dupcor$consensus.correlation)
+fit2<-contrasts.fit(fit, cont.matrix)
+fit2<-eBayes(fit2)
+summary(decideTests(fit2))
+list1<-topTable(fit2, coef=1, number="all")
+list2<-topTable(fit2, coef=2, number="all")
+list3<-topTable(fit2, coef=3, number="all")
+write.table(list1, file="C_ATRAvsC.txt", sep="\t")
+write.table(list2, file="IKK1_ATRAvsIKK1.txt", sep="\t")
+write.table(list3, file="IKK2_ATRAvsIKK2.txt", sep="\t")
+vennDiagram(results)
+results=decideTests(fit2)
+vennDiagram(results)
+savehistory("E:/PROJEKT/PA126/hist120605.Rhistory")
+save.image("E:/PROJEKT/PA126/Analysis120605.RData")
+ls()
+design2=cbind(Control=c(1,1,1,1,0,0,0,0,0,0,0,0), IKK=c(0,0,0,0,1,1,1,1,1,1,1,1), ATRA=c(0,1,0,1,0,1,0,1,0,1,0,1))
+design2=cbind(Control=c(1,1,1,1,0,0,0,0,0,0,0,0), IKK=c(0,0,0,0,1,1,1,1,1,1,1,1))
+design3=cbind(Untreated=c(1,0,1,0,1,0,1,0,1,0,1,0), ATRA=c(0,1,0,1,0,1,0,1,0,1,0,1))
+design2
+design3
+designFinal=model.matrix(~design2+design3)
+designFinal
+colnames(design) <- c("Fem.Case","Male.Case",
+"UnTreated","ATRA")
+colnames(designFinal)=c("Control", "IKK", "Untreated", "ATRA")
+colnames(designFinal)=c("Control", "IKK", "Untreated", "ATRA")
+designFinal=model.matrix(~0+design2:design3)
+designFinal
+colnames(designFinal)=c("Control", "IKK", "Untreated", "ATRA")
+designFinal
+colnames(designFinal)=c("Control.Untreated", "IKK.Untreated", "Control.ATRA", "IKK.ATRA")
+designFinal
+contrastI=makeContrasts((IKK.ATRA-Control.ATRA)-(IKK.Untreated-Control.Untreated), levels=designFinal)
+library(limma)
+contrastI=makeContrasts((IKK.ATRA-Control.ATRA)-(IKK.Untreated-Control.Untreated), levels=designFinal)
+contrastI
+source('Y:/tillHanna/COLON/Gruppanalys_Colon_apr.r')
+setwd("E:/PROJEKT/PA126")
+Data<-read.table("RMAnorm_wA_main_data.txt", sep="\t")
+block=c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2)
+block
+block=c(1,1,2,2,1,1,2,2,1,1,2,2)
+design=cbind(C=c(1,0,1,0,0,0,0,0,0,0,0,0), C_ATRA=c(0,1,0,1,0,0,0,0,0,0,0,0), IKK1=c(0,0,0,0,1,0,1,0,0,0,0,0), IKK1_ATRA=c(0,0,0,0,0,1,0,1,0,0,0,0), IKK2=c(0,0,0,0,0,0,0,0,1,0,1,0), IKK2_ATRA=c(0,0,0,0,0,0,0,0,0,1,0,1))
+deisgn
+design
+dupcor<-duplicateCorrelation(Data, design, block=block)
+library(limma)
+dupcor<-duplicateCorrelation(Data, design, block=block)
+dupcor$consensus.correlation
+cont.matrix<- makeContrasts(C_ATRAvsA=C_ATRA-C, IKK1_ATRAvsIKK1=IKK1_ATRA-IKK1, IKK2_ATRAvsIKK2=IKK2_ATRA-IKK2, levels=design)
+fit<-lmFit(Data,design, block=block, correlation=dupcor$consensus.correlation)
+fit2<-contrasts.fit(fit, cont.matrix)
+fit2<-eBayes(fit2)
+summary(decideTests(fit2))
+vennDiagram(fit2)
+results=decideTests(fit2)
+vennDiagram(results)
+cont.matrix<- makeContrasts(C_ATRAvsC=C_ATRA-C, IKK1_ATRAvsIKK1=IKK1_ATRA-IKK1, IKK2_ATRAvsIKK2=IKK2_ATRA-IKK2, levels=design)
+fit<-lmFit(Data,design, block=block, correlation=dupcor$consensus.correlation)
+fit2<-contrasts.fit(fit, cont.matrix)
+fit2<-eBayes(fit2)
+summary(decideTests(fit2))
+list1<-topTable(fit2, coef=1, number="all")
+list2<-topTable(fit2, coef=2, number="all")
+list3<-topTable(fit2, coef=3, number="all")
+write.table(list1, file="C_ATRAvsC.txt", sep="\t")
+write.table(list2, file="IKK1_ATRAvsIKK1.txt", sep="\t")
+write.table(list3, file="IKK2_ATRAvsIKK2.txt", sep="\t")
+vennDiagram(results)
+results=decideTests(fit2)
+vennDiagram(results)
+savehistory("E:/PROJEKT/PA126/hist120605.Rhistory")
+dupcorF<-duplicateCorrelation(Data, designFinal, block=block)
+dupcorF$consensus.correlation
+fitF<-lmFit(Data,designFinal, block=block, correlation=dupcorF$consensus.correlation)
+fitF2<-contrasts.fit(fitF, contrastI)
+fitF2<-eBayes(fitF2)
+summary(decideTests(fitF2))
+topTable(fitF2)
+listF=topTable(fitF2, coef=1, number='all')
+write.table(listF, file="Interaction.txt", sep="\t")
+save.image("E:\\PROJEKT\\PA126\\Analysis120827.RData")
+contrastI
+setwd("E:/PROJEKT/PA071_AS_HGK/Analysis2012")
+Data=read.table("Results_RMAnorm_50prov_uAFFX_120830_data.txt", sep="\t")
+design=cbind(NK=c(1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), NT=c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), MDK=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0), MDT=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1))
+design=cbind(NK=c(1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), NT=c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), MDK=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0), MDT=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1))
+design
+library(limma)
+cont.matrix<- makeContrasts(MDKvNK_VS_MDTvNT=(MDK-NK)-(MDT-NT)), levels=design)
+cont.matrix<- makeContrasts(MDKvNK_VS_MDTvNT=(MDK-NK)-(MDT-NT), levels=design)
+cont.matrix
+fit<-lmFit(Data,design)
+fit2<-contrasts.fit(fit, cont.matrix)
+fit2<-eBayes(fit2)
+summary(decideTests(fit2))
+topTable(fit2)
+fitR<-lmFit(Data,design, method="robust")
+fitR2<-contrasts.fit(fitR, cont.matrix)
+fitR2<-eBayes(fitR2)
+summary(decideTests(fitR2))
+topTable(fitR2)
+list=topTable(fitR2, coef=1, number='all')
+write.table(list, file="MDKvNK_VS_MDTvNT", sep"\t")
+write.table(list, file="MDKvNK_VS_MDTvNT", sep="\t")
+write.table(list, file="MDKvNK_VS_MDTvNT.txt", sep="\t")
+DataA5=read.table("Results_RMAnorm_50prov_A5_uAFFX_120830_data.txt", sep="\t")
+fitR_A5<-lmFit(DataA5,design, method="robust")
+fitR_A5_2<-contrasts.fit(fitR_A5, cont.matrix)
+fitR_A5_2<-eBayes(fitR_A5_2)
+summary(decideTests(fitR_A5_2))
+topTable(fitR_A5_2)
+topTable(fitR_A5_2, number=30)
+listA5=topTable(fitR_A5_2, coef=1, number='all')
+write.table(listA5, file="MDKvNK_VS_MDTvNT_A5.txt", sep="\t")
+savehistory("E:/PROJEKT/PA071_AS_HGK/Analysis2012/hist120830.Rhistory")
+load("C:/Users/hagor/Desktop/Analysis120918.RData")
+#quartz(file='genome_dups.png',width=14,height=14,dpi=300,type='png')
+png(file='genome_dups_Recurrence.png', width=1000,height=1000)
+layout(matrix(1:1,nrow=1,byrow=T), widths=1,heights=1)
+## Recurrence
+plot_genome_dup('Recurrence',sampleData,sampleData$Recidiv=='recidiv' | sampleData$Stadium==4 )
+## Chromothripsis
+#plot_genome_dup('Chromothripsis',sampleData,sampleData$Fragmented!='')
+## BRAF 600
+#plot_genome_dup('BRAF 600',sampleData,sampleData$BRAF_600!='wt')
+## KRAS 12 13
+#plot_genome_dup('KRAS 12/13',sampleData,sampleData$KRAS_12_13!='wt')
+## KRAS 61
+#plot_genome_dup('KRAS 61',sampleData,sampleData$KRAS_61!='wt')
+## MSI-L
+#plot_genome_dup('MSI-L',sampleData,sampleData$MSI=='MSI-L')
+## MSI-H
+#plot_genome_dup('MSI-H',sampleData,sampleData$MSI=='MSI-H')
+## PI3K Exon 9
+#plot_genome_dup('PI3K-ex9',sampleData,sampleData$PI3K_ex9!='' & sampleData$PI3K_ex9!='wt')
+## PI3K Exon 20
+#plot_genome_dup('PI3K-ex20',sampleData,sampleData$PI3K_ex20!='' & sampleData$PI3K_ex20!='wt')
+dev.off()
+mss=sampleData$MSI!='MSI-H'
+dim(mss)
+mss
+sampleData[mss]
+mss[sampleData]
+png(file='genome_dups_Recurrence.png', width=1000,height=1000)
+layout(matrix(1:1,nrow=1,byrow=T), widths=1,heights=1)
+## Recurrence
+plot_genome_dup('Recurrence',sampleData$MSI!='MSI-H',sampleData$Recidiv=='recidiv' | sampleData$Stadium==4 )
+## Chromothripsis
+#plot_genome_dup('Chromothripsis',sampleData,sampleData$Fragmented!='')
+## BRAF 600
+#plot_genome_dup('BRAF 600',sampleData,sampleData$BRAF_600!='wt')
+## KRAS 12 13
+#plot_genome_dup('KRAS 12/13',sampleData,sampleData$KRAS_12_13!='wt')
+## KRAS 61
+#plot_genome_dup('KRAS 61',sampleData,sampleData$KRAS_61!='wt')
+## MSI-L
+#plot_genome_dup('MSI-L',sampleData,sampleData$MSI=='MSI-L')
+## MSI-H
+#plot_genome_dup('MSI-H',sampleData,sampleData$MSI=='MSI-H')
+## PI3K Exon 9
+#plot_genome_dup('PI3K-ex9',sampleData,sampleData$PI3K_ex9!='' & sampleData$PI3K_ex9!='wt')
+## PI3K Exon 20
+#plot_genome_dup('PI3K-ex20',sampleData,sampleData$PI3K_ex20!='' & sampleData$PI3K_ex20!='wt')
+dev.off()
+library(xps)
+setRepositories()
+install.packages("xps")
+library(xps)
+load("C:/Users/hagor/Dropbox/Colon_snp6/RESULTAT/Analysis121002.RData")
+ls()
+Model=glm(poor~loss1p, family=binomial)
+m <- summary(Model)
+m
+help(glm)
+save.image("C:/Users/hagor/Dropbox/Colon_snp6/RESULTAT/Analysis121002.RData")
+load("E:/PROJEKT/PA106_AS_HGK/AgeStudyOkt2012/Analysis121030.RData")
+block
+design
+dupcor$consensus.correlation
+fitB<-lmFit(Data,design, block=block, correlation=dupcor$consensus.correlation)
+help(topTable)
+load("C:/Users/hagor/Dropbox/Colon_snp6/Analysis121119.RData")
+setwd("C:/Users/hagor/Dropbox/Colon_snp6")
+source('TAPS_may2012.r')
+sampleData <- load.txt('allSamples.txt')
+load('ideogram.Rdata')
+chroms=ideogram
+#load('chroms.Rdata') #### Varning! En annan chroms har nog använts till frekvensplottarna...
+chroms$before <- temp <- 0
+chroms$Chromosome=chroms$chr
+for (i in 1:24) {
+chroms$before[chroms$chr==chrom_ucsc(i)] = temp
+temp=temp+ideogram$length[chroms$chr==chrom_ucsc(i)]
+}
+#load('ideogram.Rdata')
+nSamples=nrow(sampleData)
+#for (i in 1:ncol(sampleData)) sampleData[,i][sampleData[,i]==''] <- NA
+## Load and parse all samples
+samples <- chromwise <- meanCns <- meanImba <- maxfrag <- meanFrag <- maxfragchr <- fragmat <- LCn4mCn2 <- LCn6mCn2 <- LCn6mCn2 <-LTotal <- NULL
+for (i in 1:nSamples) {
+table <- load.txt(paste('allCNs/',sampleData$Sample[i],'_segmentCN.txt',sep=''))
+table$n <- i
+table$name <- as.character(sampleData$Sample[i])
+table$stage=sampleData$Stadium[i]
+table$rec=sampleData$Recidiv[i]
+table$kras1213=sampleData$KRAS_12_13[i]
+table$kras61=sampleData$KRAS_61[i]
+table$braf=sampleData$BRAF_600[i]
+table$msi=sampleData$MSI[i]
+table$nodes=sampleData$Antal[i]
+table$pi3k9=sampleData$PI3K_ex9[i]
+table$pi3k20=sampleData$PI3K_ex20[i]
+## Calculate degree of Cn 4/6 mCn 2 for each sample
+LCn6mCn2[i]=sum(table$lengthMB[table$Cn==6 & table$mCn==2 & !is.na(table$mCn)])
+LCn4mCn2[i]=sum(table$lengthMB[table$Cn==4 & table$mCn==2 & !is.na(table$mCn)])
+LCn46mCn2[i]=LCn4mCn2[i]+LCn6mCn2[i]
+LTotal[i]=sum(table$lengthMB)
+## Calculate chrom-arm-wise Cn-changes.
+changes=ideogram; changes$p=0; changes$q=0; changes$maxfrag=0; changes$pmed=NA; changes$qmed=NA; changes$plocal=NA; changes$qlocal=NA;
+for (c in as.character(changes$chr)) {
+ix=changes$chr==c
+# p arm: breakpoints
+temp=table[table$Chromosome==c & table$End<ideogram$mid[ideogram$chr==c],]
+temp=temp[!is.na(temp$imba) & temp$snps>100,]
+if (nrow(temp)>1) {
+changes$p[ix]=sum(abs(temp$imba[1:(nrow(temp)-1)] - temp$imba[2:(nrow(temp))])>0.15)
+} else changes$p[ix]=c(0)
+# median copy number
+if (nrow(temp)>0) changes$pmed[ix]=weightedMedian(temp$Cn,temp$probes)
+# Localized chromothripsis (20MB window)
+maxp <- max <- end <- 0; while (end<ideogram$start[ix]) {
+t=temp[temp$End<end & temp$Start<end-20e6,]
+if (nrow(t)>1) {
+breaks=sum(abs(t$imba[1:(nrow(t)-1)] - t$imba[2:(nrow(t))])>0.15)
+} else breaks=0
+if (breaks>max) {
+max=breaks; maxp=end-10e6
+}
+end=end+1e6
+}; changes$plocal[ix]=max
+# q arm: breakpoints
+temp=table[table$Chromosome==c & table$Start>ideogram$mid[ideogram$chr==c],]
+temp=temp[!is.na(temp$imba) & temp$snps>100,]
+if (nrow(temp)>1) {
+changes$q[ix]=sum(abs(temp$imba[1:(nrow(temp)-1)] - temp$imba[2:(nrow(temp))])>0.05)
+} else changes$q[ix]=c(0)
+# median copy number
+if (nrow(temp)>0) changes$qmed[ix]=weightedMedian(temp$Cn,temp$probes)
+# Localized chromothripsis (20MB window)
+maxp <- max <-0; end <- ideogram$end[ix]; while (end<ideogram$length[ix]) {
+t=temp[temp$End<end & temp$Start<end-20e6,]
+if (nrow(t)>1) {
+breaks=sum(abs(t$imba[1:(nrow(t)-1)] - t$imba[2:(nrow(t))])>0.05)
+} else breaks=0
+if (breaks>max) {
+max=breaks; maxp=end-10e6
+}
+end=end+1e6
+}; changes$qlocal[ix]=max
+}
+changes$p=changes$p/ (changes$start) * 100e6
+changes$q=changes$q/ (changes$length-changes$end) * 100e6
+changes$maxfrag=apply(changes[,6:7],1,max)
+ix <- as.numeric(deChrom_ucsc(table$Chromosome)) <= 22
+table$meanCn <- meanCns[i] <- round(weightedMean(table$Cn[ix], table$lengthMB[ix]),2)
+table$meanImba <- meanImba[i] <- sum(table$lengthMB[table$Cn != table$mCn*2],na.rm=T) / 3000#1 - weightedMean(table$mCn[ix] / (table$Cn[ix]-table$mCn[ix]), table$length[ix])
+maxfrag[i] <- max(changes$maxfrag)
+maxfragchr[i] <- as.character(changes$chr[changes$maxfrag==maxfrag[i]][1])
+meanFrag[i] <- mean(mean(changes$p,na.rm=T),mean(changes$q,na.rm=T))
+samples <- rbind(samples,table)
+chromwise[[i]] = changes
+}
+sampleData$meanCn=meanCns
+sampleData$meanImba=meanImba
+sampleData$maxFrag=maxfrag
+sampleData$maxFragChr=maxfragchr
+sampleData$meanFrag=meanFrag
+## Load and parse all samples
+samples <- chromwise <- meanCns <- meanImba <- maxfrag <- meanFrag <- maxfragchr <- fragmat <- LCn4mCn2 <- LCn6mCn2 <- LCn46mCn2 <-LTotal <- NULL
+for (i in 1:nSamples) {
+table <- load.txt(paste('allCNs/',sampleData$Sample[i],'_segmentCN.txt',sep=''))
+table$n <- i
+table$name <- as.character(sampleData$Sample[i])
+table$stage=sampleData$Stadium[i]
+table$rec=sampleData$Recidiv[i]
+table$kras1213=sampleData$KRAS_12_13[i]
+table$kras61=sampleData$KRAS_61[i]
+table$braf=sampleData$BRAF_600[i]
+table$msi=sampleData$MSI[i]
+table$nodes=sampleData$Antal[i]
+table$pi3k9=sampleData$PI3K_ex9[i]
+table$pi3k20=sampleData$PI3K_ex20[i]
+## Calculate degree of Cn 4/6 mCn 2 for each sample
+LCn6mCn2[i]=sum(table$lengthMB[table$Cn==6 & table$mCn==2 & !is.na(table$mCn)])
+LCn4mCn2[i]=sum(table$lengthMB[table$Cn==4 & table$mCn==2 & !is.na(table$mCn)])
+LCn46mCn2[i]=LCn4mCn2[i]+LCn6mCn2[i]
+LTotal[i]=sum(table$lengthMB)
+## Calculate chrom-arm-wise Cn-changes.
+changes=ideogram; changes$p=0; changes$q=0; changes$maxfrag=0; changes$pmed=NA; changes$qmed=NA; changes$plocal=NA; changes$qlocal=NA;
+for (c in as.character(changes$chr)) {
+ix=changes$chr==c
+# p arm: breakpoints
+temp=table[table$Chromosome==c & table$End<ideogram$mid[ideogram$chr==c],]
+temp=temp[!is.na(temp$imba) & temp$snps>100,]
+if (nrow(temp)>1) {
+changes$p[ix]=sum(abs(temp$imba[1:(nrow(temp)-1)] - temp$imba[2:(nrow(temp))])>0.15)
+} else changes$p[ix]=c(0)
+# median copy number
+if (nrow(temp)>0) changes$pmed[ix]=weightedMedian(temp$Cn,temp$probes)
+# Localized chromothripsis (20MB window)
+maxp <- max <- end <- 0; while (end<ideogram$start[ix]) {
+t=temp[temp$End<end & temp$Start<end-20e6,]
+if (nrow(t)>1) {
+breaks=sum(abs(t$imba[1:(nrow(t)-1)] - t$imba[2:(nrow(t))])>0.15)
+} else breaks=0
+if (breaks>max) {
+max=breaks; maxp=end-10e6
+}
+end=end+1e6
+}; changes$plocal[ix]=max
+# q arm: breakpoints
+temp=table[table$Chromosome==c & table$Start>ideogram$mid[ideogram$chr==c],]
+temp=temp[!is.na(temp$imba) & temp$snps>100,]
+if (nrow(temp)>1) {
+changes$q[ix]=sum(abs(temp$imba[1:(nrow(temp)-1)] - temp$imba[2:(nrow(temp))])>0.05)
+} else changes$q[ix]=c(0)
+# median copy number
+if (nrow(temp)>0) changes$qmed[ix]=weightedMedian(temp$Cn,temp$probes)
+# Localized chromothripsis (20MB window)
+maxp <- max <-0; end <- ideogram$end[ix]; while (end<ideogram$length[ix]) {
+t=temp[temp$End<end & temp$Start<end-20e6,]
+if (nrow(t)>1) {
+breaks=sum(abs(t$imba[1:(nrow(t)-1)] - t$imba[2:(nrow(t))])>0.05)
+} else breaks=0
+if (breaks>max) {
+max=breaks; maxp=end-10e6
+}
+end=end+1e6
+}; changes$qlocal[ix]=max
+}
+changes$p=changes$p/ (changes$start) * 100e6
+changes$q=changes$q/ (changes$length-changes$end) * 100e6
+changes$maxfrag=apply(changes[,6:7],1,max)
+ix <- as.numeric(deChrom_ucsc(table$Chromosome)) <= 22
+table$meanCn <- meanCns[i] <- round(weightedMean(table$Cn[ix], table$lengthMB[ix]),2)
+table$meanImba <- meanImba[i] <- sum(table$lengthMB[table$Cn != table$mCn*2],na.rm=T) / 3000#1 - weightedMean(table$mCn[ix] / (table$Cn[ix]-table$mCn[ix]), table$length[ix])
+maxfrag[i] <- max(changes$maxfrag)
+maxfragchr[i] <- as.character(changes$chr[changes$maxfrag==maxfrag[i]][1])
+meanFrag[i] <- mean(mean(changes$p,na.rm=T),mean(changes$q,na.rm=T))
+samples <- rbind(samples,table)
+chromwise[[i]] = changes
+}
+sampleData$meanCn=meanCns
+sampleData$meanImba=meanImba
+sampleData$maxFrag=maxfrag
+sampleData$maxFragChr=maxfragchr
+sampleData$meanFrag=meanFrag
+LCn46mCn2
+Ratio=LCn46mCn2/LTotal
+Ratio=LCn46mCn2/LTotal
+SelectedRatio=Ratio[samples$msi!='MSI-H',]
+Ratio=LCn46mCn2/LTotal
+SelectedRatio=Ratio[samples$msi!='MSI-H']
+dim(SelectedRatio)
+SelectedRatio
+Ratio=LCn46mCn2/LTotal
+SelectedRatio=Ratio[Ratio(samples$msi!='MSI-H')]
+ratio=LCn46mCn2/LTotal
+selected<- samples[samples$msi!='MSI-H',]
+selected
+View(sampleData)
+View(sampleData)
+View(sampleData)
+View(sampleData)
+sampleData$MSI!='MSI-H'
+ratio[sampleData$MSI!='MSI-H']
+ratio=(LCn46mCn2/LTotal)*100
+selected<- ratio[sampleData$MSI!='MSI-H']
+hist(selected)
+hist(selected,50)
+MSS&MSI-L=selected
+MSSandMSIL=selected
+hist(MSSandMSIL,50)
+dim(selected)
+selected
+hist(selected,50)
+plot(sampleData$meanCn[sampelData$MSI!='MSI-H'], selected)
+plot(sampleData$meanCn[sampleData$MSI!='MSI-H'], selected)
+PercentageTemp=(LCn46mCn2/LTotal)*100
+Percentage<- ratio[sampleData$MSI!='MSI-H']
+hist(Percentage,50)
+PercentageTemp=(LCn46mCn2/LTotal)*100
+Percentage_Cn4or6_mCn2<- ratio[sampleData$MSI!='MSI-H']
+hist(Percentage_Cn4or6_mCn2,50)
+MeanCN_MSSandMSILsamples=sampleData$meanCn[sampleData$MSI!='MSI-H']
+plot(MeanCN_MSSandMSILsample, Percentage_Cn4or6_mCn2)
+MeanCN_MSSandMSILsamples=sampleData$meanCn[sampleData$MSI!='MSI-H']
+plot(MeanCN_MSSandMSILsamples, Percentage_Cn4or6_mCn2)
+sum(sampleData$meanCn[sampleData$MSI!='MSI-H')
+sampleData$meanCn[sampleData$MSI!='MSI-H'
+]
+PercentageTemp=(LCn46mCn2/LTotal)*100
+Percentage_Cn4or6_mCn2<- ratio[sampleData$MSI!='MSI-H']
+hist(Percentage_Cn4or6_mCn2,100)
+save.image("C:/Users/hagor/Dropbox/Colon_snp6/Analysis121119.RData")
+savehistory("C:/Users/hagor/Desktop/hist121120.Rhistory")
+view(SampleData)
+samples
+SampleData
+sampleData
+sampleData[,1]
+sampleData[1,]
+ratio
+sampleData[1,1:3]
+cbind(sampleData[,2],ratio)
+sampleData[1,1:3]
+cbind(sampleData[,3],ratio)
+cbind(sampleData[,1:3],ratio)
+RatioSampleinfo=cbind(sampleData[,1:3],ratio)
+RatioSampleinfo=cbind(sampleData[,1:3],ratio,sampleData$meanCn)
+RatioSampleinfo
+write.table(RatioSampleInfo, file="RatiosampleInfo.txt", sep="\t")
+RatioSampleInfo=cbind(sampleData[,1:3],ratio,sampleData$meanCn)
+write.table(RatioSampleInfo, file="RatiosampleInfo.txt", sep="\t")
+RatioSampleinfo
+save.image("C:/Users/hagor/Desktop/DuplicationAnalysis.RData")
+load("E:/PROJEKT/PA176/Analysis130123.RData")
+ls()
+library(limma)
+help(vennDiagram)
+result1=decideTests(fitB1R2, p=0.2)
+result2=decideTests(fitB2R2, p=0.2)
+result=decideTests(fit_R2, p=0.2)
+result
+Result=cbind(result1, result2, result)
+vennDiagram(Result)
+summary(result1)
+summary(result2)
+summary(result)
+A=[1,2,3,4,5]
+A=c(1,2,3,4,5)
+A
+hiat(A)
+hist(A)
+help(hist)
+B=A>3
+B
+load("E:/PROJEKT/PA176/Analysis130123.RData")
+ls()
+result1=decideTests(fitB1R2, p=0.2)
+result2=decideTests(fitB2R2, p=0.2)
+result=decideTests(fit_R2, p=0.2)
+Result=cbind(result1, result2, result)
+ResultNew=cbind(Result[,1:3])
+dim(Result)
+dim(ResultNew)
+ResultNew=Result[,1:3]
+dim(ResultNew)
+ResultNew=Result[,1:2]
+dim(ResultNew)
+ResultNew=Result[,1:2]
+dim(Result)
+ResultNew=Result[,1:2]
+dim(ResultNew)
+vennDiagram(ResultNew)
+load("E:/PROJEKT/PA183/AnalysBatch1and2/Analysis130311.RData")
+ls()
+results=decideTests(fitNew2)
+summary(results)
+topTable(ditNew2)
+topTable(fitNew2)
+topTable(fitNew2, coef=1)
+topTable(fitNew2, coef=2)
+topTable(fitNew2, coef=2, number=20)
+topTable(fitNew2, coef=2, number=30)
+setRepositories()
+install.packages("TAPS")
+install.packages('TAPS')
+library(TAPS)
+setRepositories()
+install.packages('TAPS')
+setwd("C:/Users/hagor/Dropbox/TAPSdev")
+setRepositories()
+install.packages('TAPS')
+setwd("C:/Users/hagor/Dropbox/TAPSdev")
+setwd("C:/Users/hagor/Dropbox/TAPSdev/TAPS")
+install.packages('TAPS')
+install.packages("C:/Users/hagor/Dropbox/TAPSdev/TAPS")
+install.packages("C:/Users/hagor/Dropbox/TAPSdev")
+install.packages("C:/Users/hagor/Dropbox/TAPSdev/TAPS.zip", repos = NULL)
+library(TAPS)
+setRepositories()
+install.packages("C:/Users/hagor/Dropbox/TAPSdev/TAPS.zip", repos = NULL)
+library(TAPS)
+library('TAPS)
+''
+)
+)')'
+library('TAPS')
Added: pkg/TAPS/DESCRIPTION
===================================================================
--- pkg/TAPS/DESCRIPTION (rev 0)
+++ pkg/TAPS/DESCRIPTION 2013-04-22 14:29:38 UTC (rev 168)
@@ -0,0 +1,11 @@
+Package: TAPS
+Type: Package
+Title: Tumor Abberation Prediction Suite
+Version: 1.0
+Date: 2013-02-11
+Author: Markus Rasmussen, Hanna Goransson-Kultima
+Maintainer: Markus Rasmussen <Markus.Mayrhofer at medsci.uu.se>
+Description: Performs a allele-specific copy number analysis of array data.
+License: GPL-2
+Depends: R (>= 2.10), DNAcopy, stats, fields, affxparser
+Packaged: 2013-03-22 10:05:27 UTC; S_D
Added: pkg/TAPS/NAMESPACE
===================================================================
--- pkg/TAPS/NAMESPACE (rev 0)
+++ pkg/TAPS/NAMESPACE 2013-04-22 14:29:38 UTC (rev 168)
@@ -0,0 +1,5 @@
+export(TAPS_plot,
+ TAPS_call,
+ TAPS_region)
+
+
Added: pkg/TAPS/R/TAPS.r
===================================================================
--- pkg/TAPS/R/TAPS.r (rev 0)
+++ pkg/TAPS/R/TAPS.r 2013-04-22 14:29:38 UTC (rev 168)
@@ -0,0 +1,3090 @@
+
+TAPS_plot <- function(directory=NULL,#xlim=c(-1,2),ylim=c(0,1),
+ bin=400) {
+#Automatically check, and if needed install, packages stats and fields
+
+#list.of.packages <- c("stats", "fields")
+#new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
+#if(length(new.packages)) install.packages(new.packages)
+
+if (is.null(directory))
+ {
+ cat("No directory supplied, using working directory.")
+ directory = "."
+ #cat("You have not assigned a directory containing one or more folders of samples for TAPS_plot to execute. \n")
+ #cat("Example: \"/user/mysamples/\" or, to run it in your current working directory, \".\" \n")
+ #directory = readline("Please supply such a directory now: ")
+ }
+
+## This function takes a directory as input, then builds short-segment TAPS scatter plots for each sample (subdirectory) in the directory.
+ setwd(directory)
+ subs <- getSubdirs()
+ if (is.null(subs)) { ## check samples = subdirectories or a single sample = current directory
+ subs=thisSubdir()
+ setwd('..')
+ }
+
+ # create SampleData file if there is none.
+ if (length(grep('SampleData.txt',dir()))==0) save.txt(data.frame(Sample=subs,cn2='',delta='',loh='',completed.analysis='no'),file='SampleData.txt')
+
+ for (i in 1:length(subs)) try( {
+ setwd(subs[i])
+ name <- subs[i]
+ cat(' ..loading', subs[i])
+
+ if(length(grep("*.cyhd.cychp",dir()))==1) ##cyhd sample
+ {
+ ##read CYCHP file from ChAS with logratio info##
+ ################################################
+ library(affxparser)
+
+ name=list.files(".",pattern="*.cychp")
+ temp=readCcg(name)
+ tempLog2=temp$dataGroups$ProbeSets$dataSets$CopyNumber$table[,1:4]
+ #Value=NULL
+ Start=tempLog2$Position
+ End=tempLog2$Position
+ Value=tempLog2$Log2Ratio
+ tempChr=tempLog2$Chromosome
+ nrows=dim(tempLog2)[1]
+ Chromosome=rep(NA,nrows)
+ for (i in 1:nrows) {
+ Chromosome[i]=paste('chr', tempChr[i], sep="")
+ }
+
+
+ Log2=as.data.frame(cbind(Chromosome, Start, End, Value))
+
+ colnames(Log2)=c("Chromosome","Start","End","Value")
+
+ Log2$Chromosome = as.character(Log2$Chromosome)
+ Log2$Chromosome[Log2$Chromosome=="chr24"] = "chrX"
+ Log2$Chromosome[Log2$Chromosome=="chr25"] = "chrY"
+ Log2$Chromosome = as.factor(Log2$Chromosome)
+
+ Log2$Start = as.integer(as.character(Log2$Start))
+ Log2$End = as.integer(as.character(Log2$End))
+ Log2$Value = as.numeric(as.character(Log2$Value))
+
+ ##read txt file from ChAS with snp info##
+ #########################################
+
+ name=list.files(".",pattern="*.cyhd.txt")
+ tempalf=read.table(file=name, header=FALSE, sep="\t", skip=12)
+ SignalA=tempalf$V4
+ SignalB=tempalf$V5
+ tempChr=tempalf$V8
+ Start=tempalf$V9
+ End=tempalf$V9
+ Chromosome=Value=rep(NA,nrows)
+ nrows=dim(tempalf)[1]
+ for (i in 1:nrows) {
+ Value[i]=SignalB[i]/(SignalB[i]+SignalA[i])
+ Chromosome[i]=paste('chr', tempChr[i], sep="")
+ }
+
+ alf=as.data.frame(cbind(Chromosome, Start, End, Value))
+
+ colnames(alf)=c("Chromosome","Start","End","Value")
+
+ alf$Start = as.integer(as.character(alf$Start))
+ alf$End = as.integer(as.character(alf$End))
+ alf$Value = as.numeric(as.character(alf$Value))
+ }
+ else
+ {
+ Log2 <- readLog2() ## Log-R
+ alf <- readAlf() ## Allele Frequency
+ }
+
+ Log2=Log2[!is.nan(Log2$Value),]
+ Log2=Log2[!is.na(Log2$Value),]
+ alf=alf[!is.nan(alf$Value),]
+ alf=alf[!is.na(alf$Value),]
+
+ segments <- readSegments() ## segments if available (CBS recommended)
+
+ cat(' ..processing')
+ if (is.null(segments)) { ## segmentation using DNA-copy if needed (must then be installed)
+ segments <- segment_DNAcopy(Log2)
+ save.txt(segments,'_segments.txt')
+ }
+
+ segments$Value <- segments$Value-median(Log2$Value) ## Median-centering
+ Log2$Value <- Log2$Value-median(Log2$Value) ## Median-centering
+
+
+ allRegions <- makeRegions(Log2, alf, segments) ## Calculates necessary data for segments (all functions are in this file)
+ save(allRegions,file='allRegions.Rdata')
+ regs <- regsFromSegs(Log2,alf,segments,bin=bin,min=5) ## Calculates the same data for shortened segments, very useful for visualization
+ save(regs,file='shortRegions.Rdata')
+
+ #Save for TAPS_region()
+ save(regs,Log2,alf,segments,file="TAPS_plot_output.Rdata")
+
+ #Clear warnings generated previously so hopefully I can see what is actually causing the program to fail.
+ #assign("last.warning", NULL, envir = baseenv())
+
+ cat('..plotting.\n')
+
+ OverviewPlot(regs$chr,regs$start,regs$end,regs$logs,regs$scores,ideogram=NULL,
+ as.character(Log2$Chromosome),Log2$Start,Log2$Value,as.character(alf$Chromosome),alf$Start,alf$Value,
+ name=name,xlim=c(-1,1.5),ylim=c(0,1))
+
+ ## Chromosome-wise plots for manual analysis
+ regions=allRegions$regions
+
+ karyotype_chroms(regs$chr,regs$start,regs$end,regs$logs,regs$scores,ideogram=NULL,
+ as.character(Log2$Chromosome),Log2$Start,Log2$Value,as.character(alf$Chromosome),alf$Start,
+ alf$Value,name=name,xlim=c(-1,2),ylim=c(0,1))
+
+ cat('..done\n')
+ setwd('..')
+ }, silent=T)
+}
+###
+
+###
+TAPS_call <- function(directory=NULL,xlim=c(-1,1),ylim=c(0,1),minseg=1,maxCn=12) {
+## TAPS_call outputs the total and minor allele copy numbers of all segments as a text file, and as images for visual confirmation.
+## sampleInfo_TAPS.txt must be present in each sample folder. If TAPS_plot could not make a good guess of the Log-R of copy number 2
+## and the Log-R difference to a deletion, you must interpret the scatter plots and edit sampleInfo_TAPS.txt.
+ if (is.null(directory))
+ {
+ cat("No directory supplied, using working directory.")
+ directory = "."
+ #cat("You have not assigned a directory containing one or more folders of samples for TAPS_call to execute. \n")
+ #cat("Example: \"/user/mysamples/\" or, to run it in your current working directory, \".\" \n")
+ #directory = readline("Please supply such a directory now: ")
+ }
+
+
+ setwd(directory)
+ #subs <- getSubdirs()
+
+ sampleData <- load.txt('SampleData.txt')
+ subs=as.character(sampleData$Sample)
+ #save.txt(sampleData,file='sampleData_old.csv')
+
+ if (is.null(subs)) {
+ subs=thisSubdir()
+ setwd('..')
+ }
+ for (i in 1:length(subs)) if (sampleData$completed.analysis[i]=='no') {
+ setwd(subs[i])
+ name <- subs[i]
+ sampleInfo <- sampleData[sampleData$Sample==subs[i],]
+ if (nrow(sampleInfo)==1) {
+
+ cat(' ..loading', subs[i])
+ Log2 <- readLog2()
+ alf <- readAlf(localDir)
+ segments <- readSegments()
+
+ segments$Value <- segments$Value-median(Log2$Value)
+ Log2$Value <- Log2$Value-median(Log2$Value)
+
+ cat(' ..processing.\n')
+
+ load('allRegions.Rdata') ## These were prepared in TAPS_plot
+ #allRegions <- makeRegions(Log2, alf, segments)
+
+ t <- findCNs(Log2,alf,allRegions,dmin=0.9,maxCn=maxCn,ceiling=1,sampleInfo=sampleInfo) ## estimates the Log-R and Allelic Imbalance Ration of all variants up to maxCn
+
+ u <- setCNs(allRegions,t$int,t$ai,maxCn) ## Assigns copy number variant for all segments
+ allRegions$regions <- u$regions
+ save.txt(u$regions,file=paste(name,'_segmentCN.txt',sep='')) ## adjacent segments with idendical copy number are merged (except over centromere) and all are saved to a text file
+ regions=allRegions$regions
+ save(t,regions,file="regions_t.Rdata")
+
+ karyotype_check(regions$Chromosome,regions$Start,regions$End,regions$log2,regions$imba,regions$Cn,regions$mCn,t,ideogram=NULL,name=name,xlim=xlim,ylim=ylim)
+
+ karyotype_chromsCN(regions$Chromosome,regions$Start,regions$End,regions$log2,
+ regions$imba,regions$Cn,regions$mCn,ideogram=NULL,
+ Log2$Chromosome,Log2$Start,Log2$Value,alf$Chromosome,
+ alf$Start,alf$Value,t,name=name,xlim=c(-1,1),ylim=c(0,1))
+
+ cat('..done\n')
+ sampleData$completed.analysis[i] <- ''
+ } else cat('Skipped',name,'\n')
+
+ setwd('..')
+ }
+ save.txt(sampleData,file='SampleData_new.csv')
+}
+###
+regsFromSegs <- function (Log2,alf, segments, bin=200,min=1) {
+## This function builds short segments and calcualtes their average Log-R and Allelic Imbalance Ratio.
+ rownames(Log2)=1:nrows(Log2)
+ rownames(alf)=1:nrows(alf)
+ regs=list('chr'=NULL,'start'=NULL,'end'=NULL,'logs'=NULL,'scores'=NULL,'het'=NULL,'hom'=NULL,'probes'=NULL,'snps'=NULL,'key1'=rep(NA,nrow(Log2)),'key2'=rep(NA,nrow(alf)))
+ n=nrow(segments)
+ s_check=NULL
+for (c in 1:n) { ## for every segment
+ tlog=Log2[Log2$Chromosome==segments$Chromosome[c],] ## Log-R on this chromosome
+ talf=alf[alf$Chromosome==segments$Chromosome[c],] ## Allele Freq on this chromosome
+ tlog=tlog[(tlog$Start>=segments$Start[c])&(tlog$Start<segments$End[c]),] ## Log-R on this segment
+ talf=talf[(talf$Start>=segments$Start[c])&(talf$Start<segments$End[c]),] ## Allele Freq on this segment
+
+ tsnps=nrow(talf) ## number of snps and probes in this segment
+ tprobes=nrow(tlog)
+ tnregs=max(trunc(tsnps/bin),1) ## Further split into how many (a minimum of 1) subsegments
+ tcuts=segments$Start[c] ## The first start pos
+ tlength=segments$End[c]-segments$Start[c] ## Length of this whole segment
+ for (i in 1:tnregs) tcuts = c(tcuts, tcuts[1]+round(i/tnregs*tlength)) ## Break the segment at these positions
+
+ for (r in 1:(tnregs)) { ## build the subsegments
+ regs$chr=c(regs$chr,as.character(segments$Chromosome[c])) ## Chromosome
+ s_=tcuts[r] ## Start
+ e_=tcuts[r+1] ## End
+ thisalf=talf[(talf$Start>=s_)&(talf$Start<=e_),] ## get the Log-R values
+ thislog=tlog[(tlog$Start>=s_)&(tlog$Start<=e_),] ## and the allele frequency
+ regs$key1[as.integer(rownames(thislog))]=length(regs$chr) ## store their positions for fast access during plotting
+ regs$key2[as.integer(rownames(thisalf))]=length(regs$chr) ## --"--
+ regs$logs=c( regs$logs, mean(thislog$Value) ) ## store average log ratio of this segment
+ regs$probes=c(regs$probes,nrow(thislog)) ## store number of probes
+ regs$snps=c(regs$snps,nrow(thisalf)) ## store number of bi-allelic probes (SNPs)
+ #regs$or_seg=c(regs$or_seg,c)
+ regs$start=c(regs$start,s_) ## store start and end positions
+ regs$end=c(regs$end,e_)
+
+ if (nrow(thisalf)>min) { ## Time to calculate Allelic Imbalance Ratio (if enough SNPs)
+ t1=sort( abs(thisalf$Value-0.5) ) ## distance from middle (het) in the allele freq pattern, ascending
+ if (length(unique(t1))>3) { ## do not attempt clustering with too few snps
+ xx=NULL
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/patchwork -r 168
More information about the Patchwork-commits
mailing list