Thanks! I wrote a small script for that (see below); I hope this is helpful. I have also written another function (OverviewPlot
) to make heatmaps straight from a DNAcopy object. You can find it on our GitHub page. Please let us know if there is anything else we can help you with.
## Load segmented data
load("/PATH/TO/segment.Rdata")
## Replace chromosome names (23, 24) from DNAcopy object with "X" and "Y"
segment.CNA.object$output$chrom[segment.CNA.object$output$chrom == 23] <- "X"
segment.CNA.object$output$chrom[segment.CNA.object$output$chrom == 24] <- "Y"
## Split by factor
segment.CNA.object$output$ID <- as.factor(segment.CNA.object$output$ID)
segments <- split(segment.CNA.object$output, segment.CNA.object$output$ID)
library(GenomicRanges)
## Use the same GC/mappability file that you used in the CopywriteR analysis
load("PATH/TO/hg19_20kb/GC_mappability.rda")
seg.values <- lapply(segments, function(x, bins) {
## Create GRanges object
x <- x[, c("chrom", "loc.start", "loc.end", "seg.mean")]
colnames(x) <- c("Seqnames", "Start", "End", "Segmean")
x[, "Start"] <- ceiling(x[, "Start"])
x <- makeGRangesFromDataFrame(x, keep.extra.columns = TRUE)
## Find overlaps with original bins used for analysis
overlaps <- findOverlaps(x, bins)
## Retrieve and return segmentation values for bins
segment.values <- vector(length = length(bins))
segment.values[subjectHits(overlaps)] <- unname(unlist(mcols(x[queryHits(overlaps)])))
return(segment.values)
}, GC.mappa.grange)
## Reformatting
seg.values <- Reduce(function(x, y) {cbind(x, y)}, seg.values)
colnames(seg.values) <- levels(segment.CNA.object$output$ID)
bins <- as(GC.mappa.grange, "data.frame")[, c("seqnames", "start", "end")]
bins <- cbind(bins, paste0(bins$seqnames, ":", paste0(bins$start, "-", bins$end)))
colnames(bins) <- c("Chromosome", "Start", "End", "Feature")
seg.values <- cbind(bins, seg.values)
## Create tracking line for viewing options in IGV
igv.track.line <- "#track viewLimits=-3:3 graphType=heatmap color=255,0,0"
## Write output
write.table(igv.track.line, file = "PATH/TO/segment_values.igv", sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
suppressWarnings(write.table(seg.values, file = "PATH/TO/segment_values.igv", append = TRUE, sep = "\t", row.names = FALSE, quote = FALSE))
Thank you! I seem to be running into trouble at the seg.values part.
All 50 warnings are the same (as above). And all of the elements in seg.values are vectors consisting of all 0s ("segments" object has data though). Not sure where I went wrong... Any thoughts?
Seqlevels are the names of the chromosomes, and none of these are in common. Could it be that you have different chromosome notations (i.e., "chr1" vs "1") in your segment.Rdata file and your GC_mappability.rda file? You can check this by typing
segment.CNA.object$output$chrom
after loading segment.Rdata, andGC.mappa.grange
after loading the GC_mappability.rda file.PERFECT! That is exactly what happened! Fixed it with a
Thank you so much!