Problem making plots
Gon1976 opened this issue · 2 comments
Hi
I have a problem making plots.
I change my chromosomes names as chr001..2..3...
I change the script to include my numbers but the plot only show the last chromosome.
This the change a made in the script (including the exact chromosomes names:
maxLevelToPlot <- 3
ratio$Ratio[ratio$Ratio>maxLevelToPlot]=maxLevelToPlot
for (i in c('003', '006', '008', '009', '013', '015', '016', '019', '023', '024', '025', '027', '028', '029', '030', '033', '034', '035', '036', '037', '038', '039', '041', '042', '043', '047', '048', '049', '050', '051', '052', '053', '055', '056', '057', '058', '059', '061', '064', '065', '067', '068', '069', '072', '073', '074', '075', '076', '077', '079', '080', '081', '083', '085', '086', '088', '089', '090', '091', '092', '093', '095', '098', '099', '100', '101', '102', '103', '104', '105', '107', '108', '110', '111', '112', '113', '114', '117', '118', '119', '123', '124', '126', '127', '128', '129', '133', '134', '135', '137', '138', '139', '140', '145', '149', '151', '152', '153', '154', '155', '156', '159', '160', '162', '163', '164', '166', '167', '168', '169', '171', '172', '173', '178', '179', '182', '183', '184', '185', '186', '187', '188', '189', '190', '193', '195', '196', '197', '198', '199', '200', '201', '203', '204', '206', '207', '208', '209', '212', '215', '216', '217', '218', '220', '221', '222', '224', '225', '226', '227', '230', '231', '232', '233', '235', '236', '239', '241', '242', '248', '250', '251', '252', '253', '254', '255', '256', '257', '258', '259', '260', '261', '262', '263', '264', '265', '266', '267', '268', '269', '270', '271', '272', '273', '274', '275', '276', '277', '278', '279', '280', '281', '282', '283', '284', '285', '286', '287', '288', '289', '290', '291', '292', '293', '294', '295', '296', '297', '298', '299', '300', '301', '302', '303', '304', '305', '306', '307', '308', '309', '310', '311', '312', '313', '314', '315', '316', '317', '318', '319', '320', '321', '322', '323', '324', '325', '326', '327', '328', '329', '330', '331', '332', '333', '334', '335', '336', '337', '338', '339', '340', '341', '342', '343', '344', '345', '346', '347', '348', '349', '350', '351', '352', '353', '354', '355', '356', '357', '358', '359', '360', '361', '362', '363', '364', '365', '366', '367', '368', '369', '370', '371', '372', '373', '374', '375', '376', '377', '378', '379', '380', '381', '382', '383', '384', '385', '386', '387', '388', '389', '390', '391', '392', '393', '394', '395', '396', '397', '398', '399', '400', '401', '402', '403', '404', '405', '406', '407', '408', '409', '410', '411', '412', '413', '414', '415', '416', '417', '418', '419', '420', '421', '422', '423', '424', '425', '426', '427', '428', '429', '430', '431', '432', '433', '434', '435', '436', '437', '438', '439', '440', '441', '442', '443', '444', '445', '446', '447', '448', '449', '450', '451', '452', '453', '454', '455', '456', '457', '458', '459', '460', '461', '462', '463', '464', '465', '466', '467', '468', '469', '470', '471', '472', '473', '474', '475', '476', '477', '478', '479', '480', '481', '482', '483', '484', '485', '486', '487', '488', '489', '490', '491', '492', '493', '494', '495', '496', '497', '498', '499', '500', '501', '502', '503', '504', '505', '506', '507', '508', '509', '510', '511', '512', '513', '514', '515', '516', '517', '518', '519', '520', '521', '522', '523', '524', '525', '526', '527', '528', '529', '530', '531', '532', '533', '534', '535', '536', '537', '538', '539', '540', '541', '542', '543', '544', '545', '546', '547', '548', '549', '550', '551', '552', '553', '554', '555', '556', '557', '558', '559', '560', '561', '562', '563', '564', '565', '566', '567', '568', '569', '570')) {
tt <- which(ratio$Chromosome==i)
if (length(tt)>0) {
plot(ratio$Start[tt],ratio$Ratio[tt]ploidy,ylim = c(0,maxLevelToPlotploidy),xlab = paste ("position, chr",i),ylab = "normalized copy number profile",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)
}
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()
} else {cat ("WARNING: To get a .png image with copy number profile, you can provide as input a file with suffix 'ratio.txt'\n")}
I am attaching my ratio file and the observed graph. As you can see, the same plot is repeated 13 times (instead of one plot for each chromosome. Can you help ? I need to change other part of the script?
update, it's work but I only have plot for last chromosomes, from 554 to 570, I lost the plot from 003 to 553
Update, my new R code is:
#!/usr/bin/env Rscript
One way to run this script is:
cat makeGraph.R | R --slave --args <_ratio.txt> [<_BAF.txt>]
Ploidy value will be inferred from the ratio file
args <- commandArgs()
BAFfileInd = 0;
ratioFileInd = 0;
#find which argument is Ratio.txt and which BAF.txt:
for (i in c(1:length(args))) {
if (length(grep("ratio.txt", args[i]))) {
ratioFileInd = i;
}
if (length(grep("BAF", args[i]))) {
BAFfileInd = i;
}
}
#------------------------------------------------------
#plot .png for the _ratio.txt file:
if (ratioFileInd) {
#read the file and get ploidy value:
ratio <-read.table(args[ratioFileInd], header=TRUE);
ratio<-data.frame(ratio)
ploidy = median (ratio$CopyNumber[which(ratio$MedianRatio>0.8 & ratio$MedianRatio<1.2)], na.rm = T)
cat (c("INFO: Selected ploidy:", ploidy, "\n"))
#------------------------------------------------------
#Plotting in the log scale:
offset = 0.01
png(filename = paste(args[ratioFileInd],".log2.png",sep = ""), width = 1180, height = 1180,
units = "px", pointsize = 20, bg = "white", res = NA)
plot(1:10)
op <- par(mfrow = c(5,5))
for (i in c('003', '006', '008', '009', '013', '015', '016', '019', '023', '024', '025', '027', '028', '029', '030', '033', '034', '035', '036', '037', '038', '039', '041', '042', '043', '047', '048', '049', '050', '051', '052', '053', '055', '056', '057', '058', '059', '061', '064', '065', '067', '068', '069', '072', '073', '074', '075', '076', '077', '079', '080', '081', '083', '085', '086', '088', '089', '090', '091', '092', '093', '095', '098', '099', '100', '101', '102', '103', '104', '105', '107', '108', '110', '111', '112', '113', '114', '117', '118', '119', '123', '124', '126', '127', '128', '129', '133', '134', '135', '137', '138', '139', '140', '145', '149', '151', '152', '153', '154', '155', '156', '159', '160', '162', '163', '164', '166', '167', '168', '169', '171', '172', '173', '178', '179', '182', '183', '184',
'185', '186', '187', '188', '189', '190', '193', '195', '196', '197', '198', '199', '200', '201', '203', '204', '206', '207', '208', '209', '212', '215', '216', '217', '218', '220', '221', '222', '224', '225', '226', '227', '230', '231', '232', '233', '235', '236', '239', '241', '242', '248', '250', '251', '252', '253', '254', '255', '256', '257', '258', '259', '260', '261', '262', '263', '264', '265', '266', '267', '268', '269', '270', '271', '272', '273', '274', '275', '276', '277', '278', '279', '280', '281', '282', '283', '284', '285', '286', '287', '288', '289', '290', '291', '292', '293', '294', '295', '296', '297', '298', '299', '300', '301', '302', '303', '304', '305', '306', '307', '308', '309', '310', '311', '312', '313', '314', '315', '316', '317', '318', '319', '320', '321', '322', '323', '324', '325', '326', '327',
'328', '329', '330', '331', '332', '333', '334', '335', '336', '337', '338', '339', '340', '341', '342', '343', '344', '345', '346', '347', '348', '349', '350', '351', '352', '353', '354', '355', '356', '357', '358', '359', '360', '361', '362', '363', '364', '365', '366', '367', '368', '369', '370', '371', '372', '373', '374', '375', '376', '377', '378', '379', '380', '381', '382', '383', '384', '385', '386', '387', '388', '389', '390', '391', '392', '393', '394', '395', '396', '397', '398', '399', '400', '401', '402', '403', '404', '405', '406', '407', '408', '409', '410', '411', '412', '413', '414', '415', '416', '417', '418', '419', '420', '421', '422', '423', '424', '425', '426', '427', '428', '429', '430', '431', '432', '433', '434', '435', '436', '437', '438', '439', '440', '441', '442', '443', '444', '445', '446', '447',
'448', '449', '450', '451', '452', '453', '454', '455', '456', '457', '458', '459', '460', '461', '462', '463', '464', '465', '466', '467', '468', '469', '470', '471', '472', '473', '474', '475', '476', '477', '478', '479', '480', '481', '482', '483', '484', '485', '486', '487', '488', '489', '490', '491', '492', '493', '494', '495', '496', '497', '498', '499', '500', '501', '502', '503', '504', '505', '506', '507', '508', '509', '510', '511', '512', '513', '514', '515', '516', '517', '518', '519', '520', '521', '522', '523', '524', '525', '526', '527', '528', '529', '530', '531', '532', '533', '534', '535', '536', '537', '538', '539', '540', '541', '542', '543', '544', '545', '546', '547', '548', '549', '550', '551', '552', '553', '554', '555', '556', '557', '558', '559', '560', '561', '562', '563', '564', '565', '566', '567',
'568', '569', '570')) {
tt <- which(ratio$Chromosome==i)
if (length(tt)>0) {
plot(ratio$Start[tt],log2(ratio$Ratio[tt]+offset),xlab = paste ("position, chr",i),ylab = "normalized copy number profile (log2)",pch = ".",col = colors()[88])
tt <- which(ratio$Chromosome==i & ratio$CopyNumber>ploidy )
points(ratio$Start[tt],log2(ratio$Ratio[tt]+offset),pch = ".",col = colors()[136])
tt <- which(ratio$Chromosome==i & ratio$CopyNumber<ploidy & ratio$CopyNumber!= -1)
points(ratio$Start[tt],log2(ratio$Ratio[tt]+offset),pch = ".",col = colors()[461])
tt <- which(ratio$Chromosome==i)
#UNCOMMENT HERE TO SEE THE PREDICTED COPY NUMBER LEVEL:
#points(ratio$Start[tt],log2(ratio$CopyNumber[tt]/ploidy+offset), 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],log2(ratio$MedianRatio[tt]+offset), pch = ".", col = colors()[463],cex=4)
}
dev.off()
#------------------------------------------------------
#Plotting in raw ratio values:
png(filename = paste(args[ratioFileInd],".png",sep = ""), width = 1180, height = 1180,
units = "px", pointsize = 20, bg = "white", res = NA)
plot(1:10)
op <- par(mfrow = c(5,5))
#replace high values of ratio with value "maxLevelToPlot":
maxLevelToPlot <- 3
ratio$Ratio[ratio$Ratio>maxLevelToPlot]=maxLevelToPlot
for (i in c('003', '006', '008', '009', '013', '015', '016', '019', '023', '024', '025', '027', '028', '029', '030', '033', '034', '035', '036', '037', '038', '039', '041', '042', '043', '047', '048', '049', '050', '051', '052', '053', '055', '056', '057', '058', '059', '061', '064', '065', '067', '068', '069', '072', '073', '074', '075', '076', '077', '079', '080', '081', '083', '085', '086', '088', '089', '090', '091', '092', '093', '095', '098', '099', '100', '101', '102', '103', '104', '105', '107', '108', '110', '111', '112', '113', '114', '117', '118', '119', '123', '124', '126', '127', '128', '129', '133', '134', '135', '137', '138', '139', '140', '145', '149', '151', '152', '153', '154', '155', '156', '159', '160', '162', '163', '164', '166', '167', '168', '169', '171', '172', '173', '178', '179', '182', '183', '184',
'185', '186', '187', '188', '189', '190', '193', '195', '196', '197', '198', '199', '200', '201', '203', '204', '206', '207', '208', '209', '212', '215', '216', '217', '218', '220', '221', '222', '224', '225', '226', '227', '230', '231', '232', '233', '235', '236', '239', '241', '242', '248', '250', '251', '252', '253', '254', '255', '256', '257', '258', '259', '260', '261', '262', '263', '264', '265', '266', '267', '268', '269', '270', '271', '272', '273', '274', '275', '276', '277', '278', '279', '280', '281', '282', '283', '284', '285', '286', '287', '288', '289', '290', '291', '292', '293', '294', '295', '296', '297', '298', '299', '300', '301', '302', '303', '304', '305', '306', '307', '308', '309', '310', '311', '312', '313', '314', '315', '316', '317', '318', '319', '320', '321', '322', '323', '324', '325', '326', '327',
'328', '329', '330', '331', '332', '333', '334', '335', '336', '337', '338', '339', '340', '341', '342', '343', '344', '345', '346', '347', '348', '349', '350', '351', '352', '353', '354', '355', '356', '357', '358', '359', '360', '361', '362', '363', '364', '365', '366', '367', '368', '369', '370', '371', '372', '373', '374', '375', '376', '377', '378', '379', '380', '381', '382', '383', '384', '385', '386', '387', '388', '389', '390', '391', '392', '393', '394', '395', '396', '397', '398', '399', '400', '401', '402', '403', '404', '405', '406', '407', '408', '409', '410', '411', '412', '413', '414', '415', '416', '417', '418', '419', '420', '421', '422', '423', '424', '425', '426', '427', '428', '429', '430', '431', '432', '433', '434', '435', '436', '437', '438', '439', '440', '441', '442', '443', '444', '445', '446', '447',
'448', '449', '450', '451', '452', '453', '454', '455', '456', '457', '458', '459', '460', '461', '462', '463', '464', '465', '466', '467', '468', '469', '470', '471', '472', '473', '474', '475', '476', '477', '478', '479', '480', '481', '482', '483', '484', '485', '486', '487', '488', '489', '490', '491', '492', '493', '494', '495', '496', '497', '498', '499', '500', '501', '502', '503', '504', '505', '506', '507', '508', '509', '510', '511', '512', '513', '514', '515', '516', '517', '518', '519', '520', '521', '522', '523', '524', '525', '526', '527', '528', '529', '530', '531', '532', '533', '534', '535', '536', '537', '538', '539', '540', '541', '542', '543', '544', '545', '546', '547', '548', '549', '550', '551', '552', '553', '554', '555', '556', '557', '558', '559', '560', '561', '562', '563', '564', '565', '566', '567',
'568', '569', '570')) {
tt <- which(ratio$Chromosome==i)
if (length(tt)>0) {
plot(ratio$Start[tt],ratio$Ratio[tt]ploidy,ylim = c(0,maxLevelToPlotploidy),xlab = paste ("position, chr",i),ylab = "normalized copy number profile",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)
}
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()
} else {cat ("WARNING: To get a .png image with copy number profile, you can provide as input a file with suffix 'ratio.txt'\n")}
#------------------------------------------------------
#plot .png for the _BAF.txt file:
if (BAFfileInd) {
BAF <-read.table(args[BAFfileInd], header=TRUE);
BAF<-data.frame(BAF)
png(filename = paste(args[BAFfileInd],".png",sep = ""), width = 1180, height = 1180,
units = "px", pointsize = 20, bg = "white", res = NA)
plot(1:10)
op <- par(mfrow = c(5,5))
for (i in c(1:22,'X','Y')) {
tt <- which(BAF$Chromosome==i)
if (length(tt)>0){
lBAF <-BAF[tt,]
plot(lBAF$Position,lBAF$BAF,ylim = c(-0.1,1.1),xlab = paste ("position, chr",i),ylab = "BAF",pch = ".",col = colors()[1])
tt <- which(lBAF$A==0.5)
points(lBAF$Position[tt],lBAF$BAF[tt],pch = ".",col = colors()[92])
tt <- which(lBAF$A!=0.5 & lBAF$A>=0)
points(lBAF$Position[tt],lBAF$BAF[tt],pch = ".",col = colors()[62])
tt <- 1
pres <- 1
if (length(lBAF$A)>4) {
for (j in c(2:(length(lBAF$A)-pres-1))) {
if (lBAF$A[j]==lBAF$A[j+pres]) {
tt[length(tt)+1] <- j
}
}
points(lBAF$Position[tt],lBAF$A[tt],pch = ".",col = colors()[24],cex=4)
points(lBAF$Position[tt],lBAF$B[tt],pch = ".",col = colors()[24],cex=4)
}
tt <- 1
pres <- 1
if (length(lBAF$FittedA)>4) {
for (j in c(2:(length(lBAF$FittedA)-pres-1))) {
if (lBAF$FittedA[j]==lBAF$FittedA[j+pres]) {
tt[length(tt)+1] <- j
}
}
points(lBAF$Position[tt],lBAF$FittedA[tt],pch = ".",col = colors()[463],cex=4)
points(lBAF$Position[tt],lBAF$FittedB[tt],pch = ".",col = colors()[463],cex=4)
}
}
}
dev.off()
} else {cat ("WARNING: To get a .png image with BAF profile, you can provide as input a file with suffix 'BAF.txt'\n")}
And my chromo names from the ratio file are:
Chromosome
003
004
006
007
....
....
....
568
569
570