search

ivanek / Gviz

This is the Gviz development repository. Gviz plots data and annotation information along genomic coordinates.
https://bioconductor.org/packages/Gviz/
Artistic License 2.0
75 stars 10 forks source link

yTicksAt not functional for DataTrack #63

Closed HeyLifeHD closed 2 years ago

HeyLifeHD commented 2 years ago

Dear Gviz authors,

I have been working on a large locus plot using Gviz and wanted to reduce the number of yTicks to make the plot more comprehensible.

I was able to set the yTicks using yTicksAt and numeric vector. This actually changed the displayPars of the individual DataTrack. However, when plotting all of the data using plotTracks and a list of different annotation and data tracks, the yTicks are set as usual (automatically calculated).

Thanks for your support,

Joschka

ivanek commented 2 years ago

Hi @HeyLifeHD,

Can you please provide an example of code?

Thanks, Robert

HeyLifeHD commented 2 years ago

Dear Robert,

This is the full code of the locus plot. The change of yAxisTicks did not work for any of the DataTracks.

`

locus plots

libraries

library(GenomicFeatures) library(rtracklayer) library(data.table) library(Gviz)

load data

annotations

anno_original_gtf <- makeTxDbFromGFF("/omics/groups/OE0219/internal/tinat/raw_data_repo/references/gencode.v29lift37.annotation.gtf")

anno_gtf <- makeTxDbFromGFF("/omics/groups/OE0219/internal/tinat/210726_shortRead_processing_deNovo_custom4/gffCompare.annotated.sorted_sub.gtf") anno <- readRDS("/omics/groups/OE0219/internal/tinat/210726_shortRead_processing_deNovo_custom4/gffCompare.annotated.sorted_repeat_anno.rds") anno_classi <- as.data.frame(data.table::fread("/omics/groups/OE0219/internal/tinat/210726_shortRead_processing_deNovo_custom4/gffCompare.mergedTranscripts.gtf.tmap"))

anno_original <- readRDS("/omics/groups/OE0219/internal/tinat/raw_data_repo/references/gencode.v29lift37.annotation.repeat_anno.rds")

repeats <- fread("/omics/groups/OE0219/internal/repguide/data/repeats/hg19_repeats.txt")

rnaseq -> merged tracks

treat1_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/treat1.bigWig") treat2_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/treat2939.bigWig") treat3_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/treat3.bigWig") treat4_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/treat4939.bigWig")

cage

treat1_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/treat1.bigWig") treat2_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/treat2939.bigWig") treat3_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/treat3.bigWig") treat4_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/treat4.bigWig")

nanopore

treat4_nanopore <- import.bw("/omics/groups/OE0219/internal/tinat/210709_nanopore_ashish_stringtie_assembly/pipeline-nanopore-ref-isoforms/bigwig/reads_aln_sorted.bw")

methylation

treat1_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150388_H2_treat1_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz") treat2_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150389_H2_treat2939_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz") treat3_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150386_H2_treat3_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz") treat4_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150387_H2_treat3_plus_treat2939_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz")

H3K9me3

treat1_9me3 <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150368_treat1_H3K9me3_ATCACG_L002_R1_complete_filtered.fastq.gz.bigWig") treat2_9me3 <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150384_treat2939_H3K9me3_TTAGGC_L002_R1_complete_filtered.fastq.gz.bigWig") treat3_9me3 <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150352_treat3_H3K9me3_CGATGT_L002_R1_complete_filtered.fastq.gz.bigWig") treat4_9me3 <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150360_treat4_H3K9me3_ACTTGA_L002_R1_complete_filtered.fastq.gz.bigWig")

H3K4me3

treat1_4me3 <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150366_treat1_H3K4me3_ATCACG_L001_R1_complete_filtered.fastq.gz.bigWig") treat2_4me3 <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150382_treat2939_H3K4me3_ACTTGA_L001_R1_complete_filtered.fastq.gz.bigWig") treat3_4me3 <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150350_treat3_H3K4me3_GCCAAT_L001_R1_complete_filtered.fastq.gz.bigWig") treat4_4me3 <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150358_treat4_H3K4me3_CTTGTA_L001_R1_complete_filtered.fastq.gz.bigWig")

H3K9ac

treat1_9ac <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150367_treat1_H3K9ac_CGATGT_L005_R1_complete_filtered.fastq.gz.bigWig") treat2_9ac <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150383_treat2939_H3K9ac_TGACCA_L005_R1_complete_filtered.fastq.gz.bigWig") treat3_9ac <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150351_treat3_H3K9ac_TTAGGC_L005_R1_complete_filtered.fastq.gz.bigWig") treat4_9ac <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150359_treat4_H3K9ac_GATCAG_L005_R1_complete_filtered.fastq.gz.bigWig")

H3K27ac

treat1_27ac <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150362_treat1_H3K27ac_ATCACG_L008_R1_complete_filtered.fastq.gz.bigWig") treat2_27ac <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150378_treat2939_H3K27ac_ACTTGA_L008_R1_complete_filtered.fastq.gz.bigWig") treat3_27ac <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150346_treat3_H3K27ac_GCCAAT_L008_R1_complete_filtered.fastq.gz.bigWig") treat4_27ac <-import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/GSM2150354_treat4_H3K27ac_CTTGTA_L008_R1_complete_filtered.fastq.gz.bigWig")

H3K14ac

chip.dir <- "/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/" treat1_14ac <- import.bw(file.path(chip.dir, paste0("GSM2597665_treat1-7.normalized", ".bigWig"))) treat2_14ac <- import.bw(file.path(chip.dir, paste0("GSM2597647_treat2-7.normalized", ".bigWig"))) treat3_14ac <- import.bw(file.path(chip.dir, paste0("GSM2597638_treat3-7.normalized", ".bigWig"))) treat4_14ac <- import.bw(file.path(chip.dir, paste0("GSM2597656_treat3+treat2-7.normalized", ".bigWig")))

H2BK5ac

treat1_5ac <- import.bw(file.path(chip.dir, paste0("GSM2597670_treat1-15.normalized", ".bigWig"))) treat2_5ac <- import.bw(file.path(chip.dir, paste0("GSM2597652_treat2-15.normalized", ".bigWig"))) treat3_5ac <- import.bw(file.path(chip.dir, paste0("GSM2597643_treat3-15.normalized", ".bigWig"))) treat4_5ac <- import.bw(file.path(chip.dir, paste0("GSM2597661_treat3+treat2-15.normalized", ".bigWig")))

prepare data

create granges from repeats file

repeats <- makeGRangesFromDataFrame(repeats, seqnames.field="genoName", start.field="genoStart", end.field="genoEnd", keep.extra.columns= TRUE) repeats$repName <- as.character(repeats$repName)

subset annos for symbol parsing

anno_original_transcript_symbol <- data.frame(row.names=anno_original[anno_original$type =="transcript",]$transcript_id,

SYMBOL= anno_original[anno_original$type =="transcript",]$transcript_name)

anno_transcript_symbol <- data.frame(row.names=anno[anno$type =="transcript",]$transcript_id, SYMBOL= anno[anno$type =="transcript",]$gene_name)

get class code information

anno_classi$class_code_simple <- NA anno_classi$class_code_simple <- ifelse(anno_classi$class_code == "=", "known", NA) anno_classi$class_code_simple <- ifelse(anno_classi$class_code %in% c("s", "x", "i", "y", "p", "u"), "non-chimeric (novel)", anno_classi$class_code_simple) anno_classi$class_code_simple <- ifelse(anno_classi$class_code %in% c("c", "k", "m", "n", "j", "e", "o"), "chimeric (novel)", anno_classi$class_code_simple) anno_transcript_class <- data.frame(row.names=anno_classi$qry_id, Class_code= anno_classi$class_code_simple)

replace with easy name for gviz function

anno_transcript_class$Class_code <- gsub("non-chimeric (novel)","NonChimeric",anno_transcript_class$Class_code,fixed=TRUE) anno_transcript_class$Class_code <- gsub("chimeric (novel)","Chimeric",anno_transcript_class$Class_code,fixed=TRUE) anno_transcript_class$Class_code <- gsub("known","Known",anno_transcript_class$Class_code,fixed=TRUE)

make granges of methylation

treat1_meth <- makeGRangesFromDataFrame(treat1_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(treat1_meth)<- treat1_meth$Smoothed_Methylation_Level_H2_treat1 score(treat1_meth)<- treat1_meth$Smoothed_Methylation_Level_H2_treat1 treat2_meth <- makeGRangesFromDataFrame(treat2_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(treat2_meth)<- treat2_meth$Smoothed_Methylation_Level_H2_treat2939 score(treat2_meth)<- treat2_meth$Smoothed_Methylation_Level_H2_treat2939 treat3_meth <- makeGRangesFromDataFrame(treat3_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(treat3_meth)<- treat3_meth$Smoothed_Methylation_Level_H2_treat3 score(treat3_meth)<- treat3_meth$Smoothed_Methylation_Level_H2_treat3 treat4_meth <- makeGRangesFromDataFrame(treat4_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(treat4_meth)<- treat4_meth$Smoothed_Methylation_Level_H2_treat3_plus_treat2939 score(treat4_meth)<- treat4_meth$Smoothed_Methylation_Level_H2_treat3_plus_treat2939

parameters

create directory

dir.create("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots")

region

ext <- 5000 fontSize <- 10

colors

col <- c("#00AFBB", "#E7B800", "#FC4E07","gray") names(col) <-c("treat3", "treat3andtreat2939", "treat1", "treat2939")

select region

temp <- "^DAPK1$" roi <- anno_original[grep(temp, anno_original$gene_name),]

roi_new <-GRanges( seqnames = unique(seqnames(roi)), ranges = IRanges(start = min(start(roi)), end = max(end(roi)))) roi_new$transcript_id <- unique(roi$gene_name)

plot regions

for (i in 1:length(roi_new)){

Define Region

lim <- c(start(roi_new[i,]), end(roi_new[i,]))
Chr<- as.character(seqnames(roi_new[i,]))
range<- GRanges(
seqnames = Chr,
ranges = IRanges(start = lim[1]-ext,
               end = lim[2]+ext))

#get ideogramm tracks
itrack <- IdeogramTrack(genome = "hg19", chromosome =Chr,fontcolor="black")

#genome axis track
getrack <- GenomeAxisTrack(fontcolor="black")

#gene annotation
#original
# grtrack_original <- GeneRegionTrack(anno_original_gtf, chromosome=Chr, transcriptAnnotation="symbol", showId=TRUE, 
#     geneSymbol=TRUE, name="Genecode",collapseTrack=TRUE, fill="darkgray",color="black", cex.feature = 0.5,
#     fontsize=fontSize-2,fontcolor.title="black", fontcolor.group="black")
# symbol(grtrack_original) <- as.character(anno_original_transcript_symbol[transcript(grtrack_original),])

#deNOvo
grtrack <- GeneRegionTrack(anno_gtf, chromosome=Chr, transcriptAnnotation="symbol", showId=TRUE, # fill="darkgray",
    geneSymbol=TRUE, name="deNovo",collapseTrack=TRUE,color="black", cex.feature = 0.75,
    fontsize=fontSize,fontcolor.title="black", fontcolor.group="black", Known="darkgray", Chimeric= "orange", NonChimeric="tomato3")

symbol(grtrack) <- as.character(anno_transcript_symbol[transcript(grtrack),])
feature(grtrack) <- anno_transcript_class[transcript(grtrack),]

#get annotation tracks of LTR12 repeats
anno_repeats <- AnnotationTrack(repeats[grep("LTR12", repeats$repName),], name="LTR12", 
    shape = "box",fill="darkgray",color="darkgray", fontsize=fontSize,fontcolor.title="black",fontcolor.feature ="black",
    chromosome=Chr,genome = "hg19", fill="darkgray", featureAnnotation="feature", cex.feature = 0.5)
feature(anno_repeats)<- repeats[grep("LTR12", repeats$repName),]$repName

#Data Tracks
#rnaseq
#get ylim
ylim_rnaseq <- c(-(max(c(max(score(treat1_rnaseq[treat1_rnaseq %over%range,])), max(score(treat2_rnaseq[treat2_rnaseq %over%range,])), 
    max(score(treat3_rnaseq[treat3_rnaseq %over%range,])), max(score(treat4_rnaseq[treat4_rnaseq %over%range,]))))/20), 
    max(c(max(score(treat1_rnaseq[treat1_rnaseq %over%range,])), max(score(treat2_rnaseq[treat2_rnaseq %over%range,])), 
    max(score(treat3_rnaseq[treat3_rnaseq %over%range,])), max(score(treat4_rnaseq[treat4_rnaseq %over%range,])))))    
yTicksAt_rnaseq <- c(0, plyr::round_any(max(ylim_rnaseq), accuracy=10, f = floor))
#prepare tracks
treat1_rnaseq_track <- DataTrack(range = treat1_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black", ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_rnaseq_track <- DataTrack(range = treat2_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\ntreat2939", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_rnaseq_track <- DataTrack(range = treat3_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_rnaseq_track <- DataTrack(range = treat4_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#cage
#get ylim
ylim_cage <- c(-(max(c(max(score(treat1_cage[treat1_cage %over%range,])), max(score(treat2_cage[treat2_cage %over%range,])), 
    max(score(treat3_cage[treat3_cage %over%range,])), max(score(treat4_cage[treat4_cage %over%range,]))))/20), max(c(max(score(treat1_cage[treat1_cage %over%range,])), max(score(treat2_cage[treat2_cage %over%range,])), 
    max(score(treat3_cage[treat3_cage %over%range,])), max(score(treat4_cage[treat4_cage %over%range,])))))
yTicksAt_cage <- c(0, plyr::round_any(max(ylim_cage), accuracy=10, f = floor))
#prepare tracks
treat1_cage_track <- DataTrack(range = treat1_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_cage_track <- DataTrack(range = treat2_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_cage_track <- DataTrack(range = treat3_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage,  yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_cage_track <- DataTrack(range = treat4_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#Methylation
treat1_meth_track <- DataTrack(range = treat1_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05), yTicksAt=c(0,1),
    type = c("histogram"), chromosome = Chr, name = "WGBS\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_meth_track <- DataTrack(range = treat2_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_meth_track <- DataTrack(range = treat3_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_meth_track <- DataTrack(range = treat4_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#H3K9me3
#get ylim
ylim_9me3 <- c(-(max(c(max(score(treat1_9me3[treat1_9me3 %over%range,])), max(score(treat2_9me3[treat2_9me3 %over%range,])), 
    max(score(treat3_9me3[treat3_9me3 %over%range,])), max(score(treat4_9me3[treat4_9me3 %over%range,]))))/20), max(c(max(score(treat1_9me3[treat1_9me3 %over%range,])), max(score(treat2_9me3[treat2_9me3 %over%range,])), 
    max(score(treat3_9me3[treat3_9me3 %over%range,])), max(score(treat4_9me3[treat4_9me3 %over%range,])))))
yTicksAt_9me3 <- c(0, plyr::round_any(max(ylim_9me3), accuracy=10, f = floor))
#prepare tracks
treat1_9me3_track <- DataTrack(range = treat1_9me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9me3,yTicksAt= yTicksAt_9me3,
    type = c("histogram"), chromosome = Chr, name = "H3K9me3\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_9me3_track <- DataTrack(range = treat2_9me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9me3,yTicksAt= yTicksAt_9me3,
    type = c("histogram"), chromosome = Chr, name = "H3K9me3\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_9me3_track <- DataTrack(range = treat3_9me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9me3,yTicksAt= yTicksAt_9me3,
    type = c("histogram"), chromosome = Chr, name = "H3K9me3\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_9me3_track <- DataTrack(range = treat4_9me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9me3,yTicksAt= yTicksAt_9me3,
    type = c("histogram"), chromosome = Chr, name = "H3K9me3\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#H3K4me3
#get ylim
ylim_4me3 <- c(-(max(c(max(score(treat1_4me3[treat1_4me3 %over%range,])), max(score(treat2_4me3[treat2_4me3 %over%range,])), 
    max(score(treat3_4me3[treat3_4me3 %over%range,])), max(score(treat4_4me3[treat4_4me3 %over%range,]))))/20), max(c(max(score(treat1_4me3[treat1_4me3 %over%range,])), max(score(treat2_4me3[treat2_4me3 %over%range,])), 
    max(score(treat3_4me3[treat3_4me3 %over%range,])), max(score(treat4_4me3[treat4_4me3 %over%range,])))))
yTicksAt_4me3 <- c(0, plyr::round_any(max(ylim_4me3), accuracy=10, f = floor))
#prepare tracks
treat1_4me3_track <- DataTrack(range = treat1_4me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_4me3, yTicksAt=yTicksAt_4me3, 
    type = c("histogram"), chromosome = Chr, name = "H3K4me3\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_4me3_track <- DataTrack(range = treat2_4me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_4me3, yTicksAt=yTicksAt_4me3, 
    type = c("histogram"), chromosome = Chr, name = "H3K4me3\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_4me3_track <- DataTrack(range = treat3_4me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_4me3, yTicksAt=yTicksAt_4me3, 
    type = c("histogram"), chromosome = Chr, name = "H3K4me3\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_4me3_track <- DataTrack(range = treat4_4me3, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_4me3, yTicksAt=yTicksAt_4me3, 
    type = c("histogram"), chromosome = Chr, name = "H3K4me3\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#H3K9ac
#get ylim
ylim_9ac <- c(-(max(c(max(score(treat1_9ac[treat1_9ac %over%range,])), max(score(treat2_9ac[treat2_9ac %over%range,])), 
    max(score(treat3_9ac[treat3_9ac %over%range,])), max(score(treat4_9ac[treat4_9ac %over%range,]))))/20), max(c(max(score(treat1_9ac[treat1_9ac %over%range,])), max(score(treat2_9ac[treat2_9ac %over%range,])), 
    max(score(treat3_9ac[treat3_9ac %over%range,])), max(score(treat4_9ac[treat4_9ac %over%range,])))))
yTicksAt_9ac <- c(0, plyr::round_any(max(ylim_9ac), accuracy=10, f = floor))
#prepare tracks
treat1_9ac_track <- DataTrack(range = treat1_9ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9ac, yTicksAt= yTicksAt_9ac,
    type = c("histogram"), chromosome = Chr, name = "H3K9ac\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_9ac_track <- DataTrack(range = treat2_9ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9ac, yTicksAt= yTicksAt_9ac,
    type = c("histogram"), chromosome = Chr, name = "H3K9ac\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_9ac_track <- DataTrack(range = treat3_9ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9ac, yTicksAt= yTicksAt_9ac,
    type = c("histogram"), chromosome = Chr, name = "H3K9ac\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_9ac_track <- DataTrack(range = treat4_9ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_9ac, yTicksAt= yTicksAt_9ac,
    type = c("histogram"), chromosome = Chr, name = "H3K9ac\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#H3K27ac
#get ylim
ylim_27ac <- c(-(max(c(max(score(treat1_27ac[treat1_27ac %over%range,])), max(score(treat2_27ac[treat2_27ac %over%range,])), 
    max(score(treat3_27ac[treat3_27ac %over%range,])), max(score(treat4_27ac[treat4_27ac %over%range,]))))/20), max(c(max(score(treat1_27ac[treat1_27ac %over%range,])), max(score(treat2_27ac[treat2_27ac %over%range,])), 
    max(score(treat3_27ac[treat3_27ac %over%range,])), max(score(treat4_27ac[treat4_27ac %over%range,])))))
yTicksAt_27ac <- c(0, plyr::round_any(max(ylim_27ac), accuracy=10, f = floor)/2, plyr::round_any(max(ylim_27ac), accuracy=10, f = floor))
#prepare tracks
treat1_27ac_track <- DataTrack(range = treat1_27ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_27ac, yTicksAt=yTicksAt_27ac,
    type = c("histogram"), chromosome = Chr, name = "H3K27ac\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_27ac_track <- DataTrack(range = treat2_27ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_27ac, yTicksAt=yTicksAt_27ac,
    type = c("histogram"), chromosome = Chr, name = "H3K27ac\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_27ac_track <- DataTrack(range = treat3_27ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_27ac, yTicksAt=yTicksAt_27ac,
    type = c("histogram"), chromosome = Chr, name = "H3K27ac\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_27ac_track <- DataTrack(range = treat4_27ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_27ac, yTicksAt=yTicksAt_27ac,
    type = c("histogram"), chromosome = Chr, name = "H3K27ac\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#H3K14ac
#get ylim
ylim_14ac <- c(-(max(c(max(score(treat1_14ac[treat1_14ac %over%range,])), max(score(treat2_14ac[treat2_14ac %over%range,])), 
    max(score(treat3_14ac[treat3_14ac %over%range,])), max(score(treat4_14ac[treat4_14ac %over%range,]))))/20), max(c(max(score(treat1_14ac[treat1_14ac %over%range,])), max(score(treat2_14ac[treat2_14ac %over%range,])), 
    max(score(treat3_14ac[treat3_14ac %over%range,])), max(score(treat4_14ac[treat4_14ac %over%range,])))))
yTicksAt_14ac <- c(0, plyr::round_any(max(ylim_14ac), accuracy=10, f = floor)/2, plyr::round_any(max(ylim_14ac), accuracy=10, f = floor))
#prepare tracks
treat1_14ac_track <- DataTrack(range = treat1_14ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_14ac, yTicksAt=yTicksAt_14ac,
    type = c("histogram"), chromosome = Chr, name = "H3K14ac\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_14ac_track <- DataTrack(range = treat2_14ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_14ac, yTicksAt=yTicksAt_14ac,
    type = c("histogram"), chromosome = Chr, name = "H3K14ac\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_14ac_track <- DataTrack(range = treat3_14ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_14ac, yTicksAt=yTicksAt_14ac,
    type = c("histogram"), chromosome = Chr, name = "H3K14ac\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_14ac_track <- DataTrack(range = treat4_14ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_14ac, yTicksAt=yTicksAt_14ac,
    type = c("histogram"), chromosome = Chr, name = "H3K14ac\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#H2BK5ac
#get ylim
ylim_5ac <- c(-(max(c(max(score(treat1_5ac[treat1_5ac %over%range,])), max(score(treat2_5ac[treat2_5ac %over%range,])), 
    max(score(treat3_5ac[treat3_5ac %over%range,])), max(score(treat4_5ac[treat4_5ac %over%range,]))))/20), max(c(max(score(treat1_5ac[treat1_5ac %over%range,])), max(score(treat2_5ac[treat2_5ac %over%range,])), 
    max(score(treat3_5ac[treat3_5ac %over%range,])), max(score(treat4_5ac[treat4_5ac %over%range,])))))
yTicksAt_5ac <- c(0, plyr::round_any(max(ylim_5ac), accuracy=10, f = floor)/2, plyr::round_any(max(ylim_5ac), accuracy=10, f = floor))
#prepare tracks
treat1_5ac_track <- DataTrack(range = treat1_5ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_5ac, yTicksAt=yTicksAt_5ac,
    type = c("histogram"), chromosome = Chr, name = "H2BK5ac\ntreat1", fill.histogram=col["treat1"],col.histogram=col["treat1"])
treat2_5ac_track <- DataTrack(range = treat2_5ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_5ac, yTicksAt=yTicksAt_5ac,
    type = c("histogram"), chromosome = Chr, name = "H2BK5ac\ntreat2", fill.histogram=col["treat2939"],col.histogram=col["treat2939"])
treat3_5ac_track <- DataTrack(range = treat3_5ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_5ac, yTicksAt=yTicksAt_5ac,
    type = c("histogram"), chromosome = Chr, name = "H2BK5ac\ntreat3", fill.histogram=col["treat3"],col.histogram=col["treat3"])
treat4_5ac_track <- DataTrack(range = treat4_5ac, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_5ac, yTicksAt=yTicksAt_5ac,
    type = c("histogram"), chromosome = Chr, name = "H2BK5ac\ntreat3+treat2", fill.histogram=col["treat3andtreat2939"],col.histogram=col["treat3andtreat2939"])

#Plot track
pdf(file.path("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots",paste0(roi_new[i,]$transcript_id, "_locus_","ext_",ext, "_v1.pdf")), width = 12, height = 24)
plotTracks(list(itrack ,getrack,grtrack, anno_repeats, # grtrack_original,
    treat1_rnaseq_track,treat2_rnaseq_track,treat3_rnaseq_track,treat4_rnaseq_track,
    treat1_cage_track,treat2_cage_track,treat3_cage_track,treat4_cage_track,
    treat1_meth_track,treat2_meth_track,treat3_meth_track,treat4_meth_track,
    treat1_9me3_track,treat2_9me3_track,treat3_9me3_track,treat4_9me3_track,
    treat1_4me3_track,treat2_4me3_track,treat3_4me3_track,treat4_4me3_track,
    treat1_9ac_track,treat2_9ac_track,treat3_9ac_track,treat4_9ac_track,
    treat1_27ac_track,treat2_27ac_track,treat3_27ac_track,treat4_27ac_track,
    treat1_14ac_track,treat2_14ac_track,treat3_14ac_track,treat4_14ac_track,
    treat1_5ac_track,treat2_5ac_track,treat3_5ac_track,treat4_5ac_track
    ),fontcolor.title="black",
    from =lim[1]-ext, to = lim[2]+ext, chromosome=Chr)
dev.off()

print(i)

}

`

Joschka

HeyLifeHD commented 2 years ago

Here is a shorter example. As you can see also for the methylation track were the yTicks are indicated as numeric vector, the visualization is not taking care of this.

locus plots

libraries

library(GenomicFeatures) library(rtracklayer) library(data.table) library(Gviz)

load data

annotations

anno_original_gtf <- makeTxDbFromGFF("/omics/groups/OE0219/internal/tinat/raw_data_repo/references/gencode.v29lift37.annotation.gtf")

anno_gtf <- makeTxDbFromGFF("/omics/groups/OE0219/internal/tinat/210726_shortRead_processing_deNovo_custom4/gffCompare.annotated.sorted_sub.gtf") anno <- readRDS("/omics/groups/OE0219/internal/tinat/210726_shortRead_processing_deNovo_custom4/gffCompare.annotated.sorted_repeat_anno.rds") anno_classi <- as.data.frame(data.table::fread("/omics/groups/OE0219/internal/tinat/210726_shortRead_processing_deNovo_custom4/gffCompare.mergedTranscripts.gtf.tmap"))

anno_original <- readRDS("/omics/groups/OE0219/internal/tinat/raw_data_repo/references/gencode.v29lift37.annotation.repeat_anno.rds")

repeats <- fread("/omics/groups/OE0219/internal/repguide/data/repeats/hg19_repeats.txt")

peptides

peptides_new <- readRDS(file.path("/omics/groups/OE0219/internal/tinat/integration/peptidomics/comparison_gene_expression", "peptides_list_new.rds"))

rnaseq -> merged tracks

DMSO_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/DMSO.normalized.bigWig") SB_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/SB939.normalized.bigWig") DAC_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/DAC.normalized.bigWig") DACSB_rnaseq <- import.bw("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/bigwig/DACSB939.normalized.bigWig")

cage

DMSO_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/DMSO.bigWig") SB_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/SB939.bigWig") DAC_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/DAC.bigWig") DACSB_cage <- import.bw("/omics/groups/OE0219/internal/tinat/raw_data_repo/CAGE/DACSB.bigWig")

nanopore

DACSB_nanopore <- import.bw("/omics/groups/OE0219/internal/tinat/210709_nanopore_ashish_stringtie_assembly/pipeline-nanopore-ref-isoforms/bigwig/reads_aln_sorted.bw")

methylation

DMSO_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/online/GSM2150388_H2_DMSO_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz") SB_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/online/GSM2150389_H2_SB939_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz") DAC_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/online/GSM2150386_H2_DAC_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz") DACSB_meth <- fread("/omics/groups/OE0219/internal/tinat/raw_data_repo/ChIP/online/GSM2150387_H2_DAC_plus_SB939_2lanes_merged.CG.ALL.call.gz.BSmooth.csv.gz")

prepare data

create granges from repeats file

repeats <- makeGRangesFromDataFrame(repeats, seqnames.field="genoName", start.field="genoStart", end.field="genoEnd", keep.extra.columns= TRUE) repeats$repName <- as.character(repeats$repName)

subset annos for symbol parsing

anno_original_transcript_symbol <- data.frame(row.names=anno_original[anno_original$type =="transcript",]$transcript_id,

SYMBOL= anno_original[anno_original$type =="transcript",]$transcript_name)

anno_transcript_symbol <- data.frame(row.names=anno[anno$type =="transcript",]$transcript_id, SYMBOL= anno[anno$type =="transcript",]$gene_name)

get class code information

anno_classi$class_code_simple <- NA anno_classi$class_code_simple <- ifelse(anno_classi$class_code == "=", "known", NA) anno_classi$class_code_simple <- ifelse(anno_classi$class_code %in% c("s", "x", "i", "y", "p", "u"), "non-chimeric (novel)", anno_classi$class_code_simple) anno_classi$class_code_simple <- ifelse(anno_classi$class_code %in% c("c", "k", "m", "n", "j", "e", "o"), "chimeric (novel)", anno_classi$class_code_simple) anno_transcript_class <- data.frame(row.names=anno_classi$qry_id, Class_code= anno_classi$class_code_simple)

replace with easy name for gviz function

anno_transcript_class$Class_code <- gsub("non-chimeric (novel)","NonChimeric",anno_transcript_class$Class_code,fixed=TRUE) anno_transcript_class$Class_code <- gsub("chimeric (novel)","Chimeric",anno_transcript_class$Class_code,fixed=TRUE) anno_transcript_class$Class_code <- gsub("known","Known",anno_transcript_class$Class_code,fixed=TRUE)

make granges of methylation

DMSO_meth <- makeGRangesFromDataFrame(DMSO_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(DMSO_meth)<- DMSO_meth$Smoothed_Methylation_Level_H2_DMSO score(DMSO_meth)<- DMSO_meth$Smoothed_Methylation_Level_H2_DMSO SB_meth <- makeGRangesFromDataFrame(SB_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(SB_meth)<- SB_meth$Smoothed_Methylation_Level_H2_SB939 score(SB_meth)<- SB_meth$Smoothed_Methylation_Level_H2_SB939 DAC_meth <- makeGRangesFromDataFrame(DAC_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(DAC_meth)<- DAC_meth$Smoothed_Methylation_Level_H2_DAC score(DAC_meth)<- DAC_meth$Smoothed_Methylation_Level_H2_DAC DACSB_meth <- makeGRangesFromDataFrame(DACSB_meth, seqnames.field="Chromosome", start.field="Start", end.field="Start",keep.extra.columns=TRUE) mcols(DACSB_meth)<- DACSB_meth$Smoothed_Methylation_Level_H2_DAC_plus_SB939 score(DACSB_meth)<- DACSB_meth$Smoothed_Methylation_Level_H2_DAC_plus_SB939

parameters

create directory

dir.create("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots/sub")

region

ext <- 1500 fontSize <- 10

get color code

library(RColorBrewer) col <- brewer.pal(12, "Paired") treat_col <- col[c(4,2,8,6)] names(treat_col)<-c("DMSO", "DAC", "SB939", "DACandSB939") class_col <- c("gray", col[c(9,10)]) names(class_col)<- c("known", "chimeric (novel)", "non-chimeric (novel)") ere_col <- col[c(1,12,5,7,3)] names(ere_col)<- c("LINE", "LTR", "no ERE", "other", "SINE") exon_col <- c("darkgray", "whitesmoke") names(exon_col)<- c("multi-exonic", "mono-exonic")

select region

temp <- "^FBP2$" roi <- anno_original[grep(temp, anno_original$gene_name),] roi_new <-GRanges( seqnames = unique(seqnames(roi)), ranges = IRanges(start = min(start(roi)), end = max(end(roi)))) roi_new$transcript_id <- unique(roi$transcript_id anno_transcript <- anno[anno$type=="transcript",]

subset peptide list based on our ORFs == Universe

peptides_new_ORF <- peptides_new[peptides_new$Species =="ORFs",]

get transcripts that give rise to uniquely presented peptides

peptides_new_ORF_noDMSO <- peptides_new_ORF[which(peptides_new_ORF$DMSO ==0),] peptides_new_ORF_noDMSO_split <- split(peptides_new_ORF_noDMSO, peptides_new_ORF_noDMSO$DAC_SB) lapply(peptides_new_ORF_noDMSO_split, length) names(peptides_new_ORF_noDMSO_split)<- c("1/3 DAC + SB939 unique", "2/3 DAC + SB939 unique", "3/3 DAC + SB939 unique") roi <- anno_transcript[anno_transcript$transcript_id %in% peptides_new_ORF_noDMSO_split$"3/3 DAC + SB939 unique"$transcript_id,] roi_new <- roi

plot regions

for (i in 1:length(roi_new)){

Define Region

lim <- c(start(roi_new[i,]), end(roi_new[i,]))
Chr<- as.character(seqnames(roi_new[i,]))
range<- GRanges(
seqnames = Chr,
ranges = IRanges(start = lim[1]-ext,
               end = lim[2]+ext))

#get ideogramm tracks
itrack <- IdeogramTrack(genome = "hg19", chromosome =Chr,fontcolor="black")

#genome axis track
getrack <- GenomeAxisTrack(fontcolor="black")

#deNOvo anno
grtrack <- GeneRegionTrack(anno_gtf, chromosome=Chr, transcriptAnnotation="symbol", showId=TRUE, # fill="darkgray",
    geneSymbol=TRUE, name="deNovo",collapseTrack=TRUE,color="black", cex.feature = 0.75,
    fontsize=fontSize,fontcolor.title="black", fontcolor.group="black", 
    Known=class_col["known"], Chimeric=class_col["chimeric (novel)"], NonChimeric=class_col["non-chimeric (novel)"])
symbol(grtrack) <- as.character(anno_transcript_symbol[transcript(grtrack),])
feature(grtrack) <- anno_transcript_class[transcript(grtrack),]

#get annotation tracks of LTR12 repeats
anno_repeats <- AnnotationTrack(repeats[grep("LTR12", repeats$repName),], name="LTR12", 
    shape = "box",fill="darkgray",color="darkgray", fontsize=fontSize,fontcolor.title="black",fontcolor.feature ="black",
    chromosome=Chr,genome = "hg19", fill="darkgray", featureAnnotation="feature", cex.feature = 0.5)
feature(anno_repeats)<- repeats[grep("LTR12", repeats$repName),]$repName

#Data Tracks
#rnaseq
#get ylim
ylim_rnaseq <- c(-(max(c(max(score(DMSO_rnaseq[DMSO_rnaseq %over%range,])), max(score(SB_rnaseq[SB_rnaseq %over%range,])), 
    max(score(DAC_rnaseq[DAC_rnaseq %over%range,])), max(score(DACSB_rnaseq[DACSB_rnaseq %over%range,]))))/20), 
    max(c(max(score(DMSO_rnaseq[DMSO_rnaseq %over%range,])), max(score(SB_rnaseq[SB_rnaseq %over%range,])), 
    max(score(DAC_rnaseq[DAC_rnaseq %over%range,])), max(score(DACSB_rnaseq[DACSB_rnaseq %over%range,])))))    
yTicksAt_rnaseq <- c(0, plyr::round_any(max(ylim_rnaseq), accuracy=10, f = floor))
#prepare tracks
DMSO_rnaseq_track <- DataTrack(range = DMSO_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black", ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nDMSO", fill.histogram=treat_col["DMSO"],col.histogram=treat_col["DMSO"])
SB_rnaseq_track <- DataTrack(range = SB_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nSB939", fill.histogram=treat_col["SB939"],col.histogram=treat_col["SB939"])
DAC_rnaseq_track <- DataTrack(range = DAC_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nDAC", fill.histogram=treat_col["DAC"],col.histogram=treat_col["DAC"])
DACSB_rnaseq_track <- DataTrack(range = DACSB_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#cage
#get ylim
ylim_cage <- c(-(max(c(max(score(DMSO_cage[DMSO_cage %over%range,])), max(score(SB_cage[SB_cage %over%range,])), 
    max(score(DAC_cage[DAC_cage %over%range,])), max(score(DACSB_cage[DACSB_cage %over%range,]))))/20), max(c(max(score(DMSO_cage[DMSO_cage %over%range,])), max(score(SB_cage[SB_cage %over%range,])), 
    max(score(DAC_cage[DAC_cage %over%range,])), max(score(DACSB_cage[DACSB_cage %over%range,])))))
yTicksAt_cage <- c(0, plyr::round_any(max(ylim_cage), accuracy=10, f = floor))
#prepare tracks
DMSO_cage_track <- DataTrack(range = DMSO_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nDMSO", fill.histogram=treat_col["DMSO"],col.histogram=treat_col["DMSO"])
SB_cage_track <- DataTrack(range = SB_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nSB", fill.histogram=treat_col["SB939"],col.histogram=treat_col["SB939"])
DAC_cage_track <- DataTrack(range = DAC_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage,  yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nDAC", fill.histogram=treat_col["DAC"],col.histogram=treat_col["DAC"])
DACSB_cage_track <- DataTrack(range = DACSB_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#Methylation
DMSO_meth_track <- DataTrack(range = DMSO_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05), yTicksAt=c(0,1),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nDMSO", fill.histogram=treat_col["DMSO"],col.histogram=treat_col["DMSO"])
SB_meth_track <- DataTrack(range = SB_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nSB", fill.histogram=treat_col["SB939"],col.histogram=treat_col["SB939"])
DAC_meth_track <- DataTrack(range = DAC_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nDAC", fill.histogram=treat_col["DAC"],col.histogram=treat_col["DAC"])
DACSB_meth_track <- DataTrack(range = DACSB_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#Plot track
pdf(file.path("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots/sub",paste0(roi_new[i,]$transcript_id, "_locus_","ext_",ext, "_v1.pdf")), width = 12, height = 8)
plotTracks(list(itrack ,getrack,grtrack, anno_repeats, # grtrack_original,
    DMSO_rnaseq_track,SB_rnaseq_track,DAC_rnaseq_track,DACSB_rnaseq_track,
    DMSO_cage_track,SB_cage_track,DAC_cage_track,DACSB_cage_track,
    DMSO_meth_track,SB_meth_track,DAC_meth_track,DACSB_meth_track
    ),fontcolor.title="black",
    from =lim[1]-ext, to = lim[2]+ext, chromosome=Chr)
dev.off()

print(i)

}

with nanopore data

dir.create("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots/sub_withNano") dir.create("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots/sub_withAlign")

select region

temp <- "^FBP2$" roi <- anno_original[grep(temp, anno_original$gene_name),] roi_new <-GRanges( seqnames = unique(seqnames(roi)), ranges = IRanges(start = min(start(roi)), end = max(end(roi)))) roi_new$transcript_id <- unique(roi$transcript_id)[2] anno_transcript <- anno[anno$type=="transcript",]

subset peptide list based on our ORFs == Universe

peptides_new_ORF <- peptides_new[peptides_new$Species =="ORFs",]

get transcripts that give rise to uniquely presented peptides

peptides_new_ORF_noDMSO <- peptides_new_ORF[which(peptides_new_ORF$DMSO ==0),] peptides_new_ORF_noDMSO_split <- split(peptides_new_ORF_noDMSO, peptides_new_ORF_noDMSO$DAC_SB) lapply(peptides_new_ORF_noDMSO_split, length) names(peptides_new_ORF_noDMSO_split)<- c("1/3 DAC + SB939 unique", "2/3 DAC + SB939 unique", "3/3 DAC + SB939 unique") roi <- anno_transcript[anno_transcript$transcript_id %in% peptides_new_ORF_noDMSO_split$"3/3 DAC + SB939 unique"$transcript_id,] roi_new <- roi

plot regions

for (i in 1:length(roi_new)){

Define Region

lim <- c(start(roi_new[i,]), end(roi_new[i,]))
Chr<- as.character(seqnames(roi_new[i,]))
range<- GRanges(
seqnames = Chr,
ranges = IRanges(start = lim[1]-ext,
               end = lim[2]+ext))

#get ideogramm tracks
itrack <- IdeogramTrack(genome = "hg19", chromosome =Chr,fontcolor="black")

#genome axis track
getrack <- GenomeAxisTrack(fontcolor="black")

#deNOvo anno
grtrack <- GeneRegionTrack(anno_gtf, chromosome=Chr, transcriptAnnotation="symbol", showId=TRUE, # fill="darkgray",
    geneSymbol=TRUE, name="deNovo",collapseTrack=TRUE,color="black", cex.feature = 0.75,
    fontsize=fontSize,fontcolor.title="black", fontcolor.group="black", 
    Known=class_col["known"], Chimeric=class_col["chimeric (novel)"], NonChimeric=class_col["non-chimeric (novel)"])
symbol(grtrack) <- as.character(anno_transcript_symbol[transcript(grtrack),])
feature(grtrack) <- anno_transcript_class[transcript(grtrack),]

#get annotation tracks of LTR12 repeats
anno_repeats <- AnnotationTrack(repeats[grep("LTR12", repeats$repName),], name="LTR12", 
    shape = "box",fill="darkgray",color="darkgray", fontsize=fontSize,fontcolor.title="black",fontcolor.feature ="black",
    chromosome=Chr,genome = "hg19", fill="darkgray", featureAnnotation="feature", cex.feature = 0.5)
feature(anno_repeats)<- repeats[grep("LTR12", repeats$repName),]$repName

#Data Tracks
#rnaseq
#get ylim
ylim_rnaseq <- c(-(max(c(max(score(DMSO_rnaseq[DMSO_rnaseq %over%range,])), max(score(SB_rnaseq[SB_rnaseq %over%range,])), 
    max(score(DAC_rnaseq[DAC_rnaseq %over%range,])), max(score(DACSB_rnaseq[DACSB_rnaseq %over%range,]))))/20), 
    max(c(max(score(DMSO_rnaseq[DMSO_rnaseq %over%range,])), max(score(SB_rnaseq[SB_rnaseq %over%range,])), 
    max(score(DAC_rnaseq[DAC_rnaseq %over%range,])), max(score(DACSB_rnaseq[DACSB_rnaseq %over%range,])))))    
yTicksAt_rnaseq <- c(0, plyr::round_any(max(ylim_rnaseq), accuracy=10, f = floor))
#prepare tracks
DMSO_rnaseq_track <- DataTrack(range = DMSO_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black", ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nDMSO", fill.histogram=treat_col["DMSO"],col.histogram=treat_col["DMSO"])
SB_rnaseq_track <- DataTrack(range = SB_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nSB939", fill.histogram=treat_col["SB939"],col.histogram=treat_col["SB939"])
DAC_rnaseq_track <- DataTrack(range = DAC_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nDAC", fill.histogram=treat_col["DAC"],col.histogram=treat_col["DAC"])
DACSB_rnaseq_track <- DataTrack(range = DACSB_rnaseq, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "RNAseq\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#nanopore
nano_track <- DataTrack(range = DACSB_nanopore, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",#ylim=ylim_rnaseq, yTicksAt=yTicksAt_rnaseq,
    type = c("histogram"), chromosome = Chr, name = "Nano\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#bam track for dacsb
DAC_SB_alignment <- AlignmentsTrack("/omics/groups/OE0219/internal/tinat/210712_shortRead_processing_knownRef/HISAT2/aligned_sorted/AS-641779-LR-57123.sorted.bam",
    isPaired = TRUE, type=c("sashimi"),col.sashimi=treat_col["DACandSB939"],
    fontsize=fontSize,fontcolor.title="black",col.axis="black",  sashimiFilter =  intronicParts(anno_gtf) ,sashimiFilterTolerance = 10000L)

#cage
#get ylim
ylim_cage <- c(-(max(c(max(score(DMSO_cage[DMSO_cage %over%range,])), max(score(SB_cage[SB_cage %over%range,])), 
    max(score(DAC_cage[DAC_cage %over%range,])), max(score(DACSB_cage[DACSB_cage %over%range,]))))/20), max(c(max(score(DMSO_cage[DMSO_cage %over%range,])), max(score(SB_cage[SB_cage %over%range,])), 
    max(score(DAC_cage[DAC_cage %over%range,])), max(score(DACSB_cage[DACSB_cage %over%range,])))))
yTicksAt_cage <- c(0, plyr::round_any(max(ylim_cage), accuracy=10, f = floor))
#prepare tracks
DMSO_cage_track <- DataTrack(range = DMSO_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nDMSO", fill.histogram=treat_col["DMSO"],col.histogram=treat_col["DMSO"])
SB_cage_track <- DataTrack(range = SB_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nSB", fill.histogram=treat_col["SB939"],col.histogram=treat_col["SB939"])
DAC_cage_track <- DataTrack(range = DAC_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage,  yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nDAC", fill.histogram=treat_col["DAC"],col.histogram=treat_col["DAC"])
DACSB_cage_track <- DataTrack(range = DACSB_cage, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=ylim_cage, yTicksAt=yTicksAt_cage,
    type = c("histogram"), chromosome = Chr, name = "CAGE\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#Methylation
DMSO_meth_track <- DataTrack(range = DMSO_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05), yTicksAt=c(0,1),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nDMSO", fill.histogram=treat_col["DMSO"],col.histogram=treat_col["DMSO"])
SB_meth_track <- DataTrack(range = SB_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nSB", fill.histogram=treat_col["SB939"],col.histogram=treat_col["SB939"])
DAC_meth_track <- DataTrack(range = DAC_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nDAC", fill.histogram=treat_col["DAC"],col.histogram=treat_col["DAC"])
DACSB_meth_track <- DataTrack(range = DACSB_meth, genome = "hg19",
    fontsize=fontSize,fontcolor.title="black",col.axis="black",ylim=c(-.05,1.05),
    type = c("histogram"), chromosome = Chr, name = "WGBS\nDAC+SB", fill.histogram=treat_col["DACandSB939"],col.histogram=treat_col["DACandSB939"])

#Plot track
pdf(file.path("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots/sub_withNano",paste0(roi_new[i,]$transcript_id, "_locus_","ext_",ext, "_v1.pdf")), width = 12, height = 8)
plotTracks(list(itrack ,getrack,grtrack, anno_repeats, # grtrack_original,
    DMSO_rnaseq_track,SB_rnaseq_track,DAC_rnaseq_track,DACSB_rnaseq_track,nano_track,
    DMSO_cage_track,SB_cage_track,DAC_cage_track,DACSB_cage_track,
    DMSO_meth_track,SB_meth_track,DAC_meth_track,DACSB_meth_track
    ),fontcolor.title="black",sizes=c(0.5,1,2,0.5,1,1,1,1,1,1,1,1,1,1,1,1,1),
    from =lim[1]-ext, to = lim[2]+ext, chromosome=Chr)
dev.off()

pdf(file.path("/omics/groups/OE0219/internal/tinat/210727_shortRead_processing_deNovo_custom4_quantification_analysis/locus_plots/sub_withAlign",paste0(roi_new[i,]$transcript_id, "_locus_","ext_",ext, "_v1.pdf")), width = 12, height = 8)
plotTracks(list(itrack ,getrack,grtrack, anno_repeats, # grtrack_original,
    DMSO_rnaseq_track,SB_rnaseq_track,DAC_rnaseq_track,DACSB_rnaseq_track,DAC_SB_alignment,
    DMSO_cage_track,SB_cage_track,DAC_cage_track,DACSB_cage_track,
    DMSO_meth_track,SB_meth_track,DAC_meth_track,DACSB_meth_track
    ),fontcolor.title="black",sizes=c(0.5,1,2,0.5,1,1,1,1,2.5,1,1,1,1,1,1,1,1),
    from =lim[1]-ext, to = lim[2]+ext, chromosome=Chr)
dev.off()
print(i)

}

ivanek commented 2 years ago

Hi @HeyLifeHD,

I think you are not using it in a correct way:

library(Gviz)

set.seed(0)
gr1 <- GRanges(seqnames = "chrI", ranges = IRanges(start = 1:10, width = 1), score = sample(x = 1:10 / 10, size = 10, replace = TRUE))
gr2 <- GRanges(seqnames = "chrI", ranges = IRanges(start = 1:10, width = 1), score = sample(x = 1:10 / 10, size = 10, replace = TRUE))
gr3 <- GRanges(seqnames = "chrI", ranges = IRanges(start = 1:10, width = 1), score = sample(x = 1:10 / 10, size = 10, replace = TRUE))

dt1 <- DataTrack(range = gr1, type = "l", col = "black")
dt2 <- DataTrack(range = gr2, type = "l", col = "blue")
dt3 <- DataTrack(range = gr3, type = "l", col = "red")

Default behavior - both ylim and yTicksAt is set automatically by Gviz

plotTracks(trackList = list(dt1, dt2, dt3), chromosome = "chrI")

Rplot01

Setting only yTicksAt

dt1 <- DataTrack(range = gr1, type = "l", col = "black", yTicksAt=seq(0.3,0.7,0.1))
plotTracks(trackList = list(dt1, dt2, dt3), chromosome = "chrI")

Rplot02

Setting both yTicksAt and ylim

dt1 <- DataTrack(range = gr1, type = "l", col = "black", yTicksAt=seq(0,1,0.1), ylim=c(0,1))
plotTracks(trackList = list(dt1, dt2, dt3), chromosome = "chrI")

Rplot03

Hope that helps. Closing it.

HeyLifeHD commented 2 years ago

Hi @ivanek,

Thanks for the response. I see how the code is working for your script. However, my ylim as well as yTicksAt both end up to be a two numeric long vector, as I want only the maximum and minimum to be indicated in the tracks. Still, the default ticks are used.

Best,

ivanek commented 2 years ago

Do you mean something like this? 

dt1 <- DataTrack(range = gr1, type = "l", col = "black", yTicksAt=c(0,1), ylim=c(0,1))
plotTracks(trackList = list(dt1, dt2, dt3), chromosome = "chrI")

Rplot04