search

BoevaLab / FREEC

Control-FREEC: Copy number and genotype annotation in whole genome and whole exome sequencing data
147 stars 49 forks source link

Problem making plots #133

Open Gon1976 opened 1 year ago

Gon1976 commented 1 year ago

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 JJchr ratio txt log2

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?

Gon1976 commented 1 year ago

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

Gon1976 commented 1 year ago

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