Problem making plots #133

Gon1976 · 2023-07-27T14:06:58Z

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")}

JJchr.ratio.txt

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?

The text was updated successfully, but these errors were encountered:

Gon1976 · 2023-08-07T17:50:43Z

update, it's work but I only have plot for last chromosomes, from 554 to 570, I lost the plot from 003 to 553

Gon1976 · 2023-08-22T17:23:57Z

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Problem making plots #133

Problem making plots #133

Gon1976 commented Jul 27, 2023

Gon1976 commented Aug 7, 2023

Gon1976 commented Aug 22, 2023

Problem making plots #133

Problem making plots #133

Comments

Gon1976 commented Jul 27, 2023

Gon1976 commented Aug 7, 2023

Gon1976 commented Aug 22, 2023

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