freec_CNVPlot.R 5.93 KB
Newer Older
1
library('optparse')
Wai Yi Leung's avatar
Wai Yi Leung committed
2 3
library('naturalsort')

4 5
# Script taken from  http://bioinfo-out.curie.fr/projects/freec/tutorial.html and modified for biopet

6
option_list <- list(
7
    make_option(c("-m", "--mappability"), dest="mappability"),
8 9 10 11
    make_option(c("-p", "--ploidy"), default=2, type="integer", dest="ploidy"),
    make_option(c("-i", "--input"), dest="input"),
    make_option(c("-o", "--output"), dest="output")
    )
12

13 14
parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
opt = parse_args(parser)
15 16


17 18 19 20 21 22 23 24 25 26 27
#
# Load mappability track
#

mappabilityFile <- opt$mappability
mappabilityTrack <- read.table(mappabilityFile, header=FALSE, col.names=c("chrom", "start", "end", "score"))

mappabilityTrack$Start <- mappabilityTrack$start+1
mappabilityTrack$Chromosome <- gsub("chr", "", mappabilityTrack$chrom)


28 29 30 31
#
# Load Data
#

32

33 34 35 36 37 38 39 40 41 42 43
dataTable <- read.table( opt$input , header=TRUE)
input_ratio <- data.frame(dataTable)

chromosomes <- naturalsort(levels(input_ratio$Chromosome))
input_ratio$Chromosome <- factor(input_ratio$Chromosome, levels=chromosomes, ordered=T)

sorted_ratio <- input_ratio[order(input_ratio$Chromosome),]
ratio <- merge(sorted_ratio, mappabilityTrack, sort=TRUE)
ratio <- ratio[order(ratio$Chromosome, ratio$Start),]


44 45 46 47
ploidy <- opt$ploidy
ppi <- 300
plot_margins <- c(3,4,1,2)+0.1
label_positions <- c(2,0.5,0)
48 49 50 51 52


maxLevelToPlot <- 3
for (i in c(1:length(ratio$Ratio))) {
	if (ratio$Ratio[i]>maxLevelToPlot) {
Wai Yi Leung's avatar
Wai Yi Leung committed
53
		ratio$Ratio[i]=maxLevelToPlot
54 55 56
	}
}

57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
#
# Plot the graphs per chromosome
#

for (i in chromosomes) {

    png(filename = paste(opt$output, ".", i,".png",sep=""), width = 4 * ppi, height = 2.5 * ppi,
        res=ppi, bg = "white")
    par(mfrow = c(1,1))
    par(mar=plot_margins)
    par(mgp=label_positions)


    tt <- which(ratio$Chromosome==i)
    if (length(tt)>0) {
        plot(ratio$Start[tt],
                ratio$Ratio[tt]*ploidy,
                ylim = c(0,maxLevelToPlot*ploidy),
                xlab = paste ("position, chr",i),
                ylab = "normalized CN",
                pch = ".",
                col = colors()[88])

        title(outer=TRUE)
        tt <- which(ratio$Chromosome==i  & ratio$CopyNumber>ploidy )
        points(ratio$Start[tt],ratio$Ratio[tt]*ploidy,pch = ".",col = colors()[136])

        tt <- which(ratio$Chromosome==i  & ratio$Ratio==maxLevelToPlot & ratio$CopyNumber>ploidy)
        points(ratio$Start[tt],ratio$Ratio[tt]*ploidy,pch = ".",col = colors()[136],cex=4)

        tt <- which(ratio$Chromosome==i  & ratio$CopyNumber<ploidy & ratio$CopyNumber!= -1)
        points(ratio$Start[tt],ratio$Ratio[tt]*ploidy,pch = ".",col = colors()[461])
        tt <- which(ratio$Chromosome==i)

        #UNCOMMENT HERE TO SEE THE PREDICTED COPY NUMBER LEVEL:
        #points(ratio$Start[tt],ratio$CopyNumber[tt], pch = ".", col = colors()[24],cex=4)
    }
    #tt <- which(ratio$Chromosome==i)

	#UNCOMMENT HERE TO SEE THE EVALUATED MEDIAN LEVEL PER SEGMENT:
	#points(ratio$Start[tt],ratio$MedianRatio[tt]*ploidy, pch = ".", col = colors()[463],cex=4)

    dev.off()
}










111
png(filename = paste(opt$output, ".png",sep=""), width = 16 * ppi, height = 10 * ppi,
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
    res=ppi, bg = "white")
par(mfrow = c(6,4))
par(mar=plot_margins)
par(mgp=label_positions)

for (i in chromosomes) {
    tt <- which(ratio$Chromosome==i)
    if (length(tt)>0) {
        plot(ratio$Start[tt],
                ratio$Ratio[tt]*ploidy,
                ylim = c(0,maxLevelToPlot*ploidy),
                xlab = paste ("position, chr",i),
                ylab = "normalized CN",
                pch = ".",
                col = colors()[88])

        tt <- which(ratio$Chromosome==i  & ratio$CopyNumber>ploidy )
        points(ratio$Start[tt],ratio$Ratio[tt]*ploidy,pch = ".",col = colors()[136])

        tt <- which(ratio$Chromosome==i  & ratio$Ratio==maxLevelToPlot & ratio$CopyNumber>ploidy)
        points(ratio$Start[tt],ratio$Ratio[tt]*ploidy,pch = ".",col = colors()[136],cex=4)

        tt <- which(ratio$Chromosome==i  & ratio$CopyNumber<ploidy & ratio$CopyNumber!= -1)
        points(ratio$Start[tt],ratio$Ratio[tt]*ploidy,pch = ".",col = colors()[461])
        tt <- which(ratio$Chromosome==i)

        #UNCOMMENT HERE TO SEE THE PREDICTED COPY NUMBER LEVEL:
        #points(ratio$Start[tt],ratio$CopyNumber[tt], pch = ".", col = colors()[24],cex=4)
140
	}
141 142
	#tt <- which(ratio$Chromosome==i)

143 144
	#UNCOMMENT HERE TO SEE THE EVALUATED MEDIAN LEVEL PER SEGMENT:
	#points(ratio$Start[tt],ratio$MedianRatio[tt]*ploidy, pch = ".", col = colors()[463],cex=4)
145

146 147
}

148 149 150 151
dev.off()



Wai Yi Leung's avatar
Wai Yi Leung committed
152 153 154 155
# Export the whole genome graph

png(filename = paste(opt$output, ".wg.png",sep=""), width = 16 * ppi, height = 10 * ppi,
res=ppi, bg = "white")
156 157 158 159 160

plot_margins <- c(3,4,2,2)+0.1
label_positions <- c(2,0.5,0)

par(mfrow = c(1,1))
Wai Yi Leung's avatar
Wai Yi Leung committed
161 162 163 164 165 166 167 168 169 170 171 172 173
par(mar=plot_margins)
par(mgp=label_positions)
par(xaxs="i", yaxs="i")


maxLevelToPlot <- 3
for (i in c(1:length(ratio$Ratio))) {
    if (ratio$Ratio[i]>maxLevelToPlot) {
        ratio$Ratio[i]=maxLevelToPlot
    }
}

for (i in c(1:length(ratio$Start))) {
174
    ratio$Position[i] = (i-1) *5000 +1
Wai Yi Leung's avatar
Wai Yi Leung committed
175 176
}

177 178 179 180 181 182 183

plotRatioLT <- 0.10

filteredSet <- ratio[ ratio$score > plotRatioLT, ]

plot(filteredSet$Position,
filteredSet$Ratio*ploidy,
Wai Yi Leung's avatar
Wai Yi Leung committed
184 185 186 187 188 189 190 191
ylim = c(0,maxLevelToPlot*ploidy),
xlab = paste ("Chr. on genome"),
ylab = "normalized CN",
pch = ".",
col = colors()[88])


title(outer=TRUE)
192 193
tt <- which(filteredSet$CopyNumber>ploidy)
points(filteredSet$Position[tt],filteredSet$Ratio[tt]*ploidy,pch = ".",col = colors()[136])
Wai Yi Leung's avatar
Wai Yi Leung committed
194

195 196 197 198 199
tt <- which(filteredSet$Ratio==maxLevelToPlot & filteredSet$CopyNumber>ploidy)
points(filteredSet$Position[tt],filteredSet$Ratio[tt]*ploidy,pch = ".",col = colors()[136],cex=4)

tt <- which(filteredSet$CopyNumber<ploidy & filteredSet$CopyNumber!= -1)
points(filteredSet$Position[tt],filteredSet$Ratio[tt]*ploidy,pch = ".",col = colors()[461], bg="black")
Wai Yi Leung's avatar
Wai Yi Leung committed
200 201 202


for (chrom in chromosomes) {
203 204 205
    tt <- which(filteredSet$Chromosome == chrom)
    print(filteredSet[tt[1],])
    xpos <- filteredSet$Position[tt][1]
Wai Yi Leung's avatar
Wai Yi Leung committed
206
    abline(v=xpos, col="grey")
207 208
    axis(3, at=xpos, labels=chrom , las=2)
}
Wai Yi Leung's avatar
Wai Yi Leung committed
209 210


211
dev.off()