Question

Givz Doesn't plot coverage genetracks that are aligned against ensemble ref

0

Entering edit mode

ta_awwad ▴ 10

@ta_awwad-11382

Last seen 6 weeks ago

Frankfurt am Main

Hello everybody,

I am using Gviz to generate snapshots of my RNA-seq experiment. I did the alignment of my bam against ensemble reference genome.

I am using the following script:

from <- 57290000
to <- 57400000
chr <- 'chr8'
gen <- 'mm10'
mut <- "mutant.bw"
wt <- "WT.bw"
mutant <- DataTrack(range = mut, genome = gen,chromosome = chr,type = "h", col="#ff0000", ylim=c(0, 1), name="Muatnt",col.sampleNames="black", background.title="white",col.axis="black", col.title="black")

wildtype <- DataTrack(range = wt, genome = gen,chromosome = chr,type = "h", col="#4e5b6e", ylim=c(0, 1), name="Wildtype", col.sampleNames="black", background.title="white",col.axis="black", col.title="black")

ensGenes <- UcscTrack(genome= gen, chromosome= chr, 
                   track= "GENCODE VM11 (Ensembl 86)",  from= from, to= to, 
                   trackType= "GeneRegionTrack",rstarts= "exonStarts",
                   rends= "exonEnds", gene= "name", 
                   symbol= "name2", transcript= "name", strand= "strand", 
                   fill= "#960000", name= "Ensembl Genes")

axTrack <- GenomeAxisTrack(col="#000000", fontcolor="#000000")

idxTrack <- IdeogramTrack(genome= gen, chromosome=chr, 
                          fontcolor= "#000000", outline=TRUE)

plotTracks(list(idxTrack, axTrack, ensGenes ,mutant, wildtype),from=from, to=to, showTitle=T, showId=T, fontsize = 18)

however, I am getting only the ensemble gene track and my bigwig track is empty.

any ideas how to solve this?

thanks much

gviz • 1.8k views

ADD COMMENT • link 6.5 years ago ta_awwad ▴ 10

score 0 · Answer 1 · 2017-10-23

Are you sure that the chromosome names match? Sometimes you will only have the numeric part (e.g., 8) in your BAM file. A bit hard to say what is going wrong without seeing the actual data.

Also, have you considered using the somewhat more specialised AlignmentsTrack class? It will be way more efficient for you BAM alignments, plus you get the added benefit of the individual read alignments if you need them

score 0 · Answer 2 · 2017-10-23

0

Entering edit mode

ta_awwad ▴ 10

@ta_awwad-11382

Last seen 6 weeks ago

Frankfurt am Main

Thanks much Florian for your answer. actually the problem seems to be chromosome names ..ensembl doesn't add "chr" thing to chromosome identifier .. this might be the problem. I did transform my bam to bigwig with STAR aligner, this means that most probably the chromosome names did not change ..

do you have any idea how to overcome this point?

ADD COMMENT • link 6.5 years ago ta_awwad ▴ 10

0

Entering edit mode

You can either use gene model based on ensembl or simply replace the seqlevels in your ensGenes.

seqlevels(ranges(ensGenes)) <- sub("^chr", "", seqlevels(ranges(ensGenes)))

ADD REPLY • link 6.5 years ago Robert Ivanek ▴ 730

0

Entering edit mode

thanks Robert .. I did this but still not working .. tried on IGV it is working fine I could see my tracks... any further suggestion?

ADD REPLY • link 6.5 years ago ta_awwad ▴ 10

0

Entering edit mode

Obviously you will need to use the numerical chromosome name now also for plotting:

plotTracks(list(idxTrack, axTrack, ensGenes ,mutant, wildtype),from=from, to=to, showTitle=T, showId=T, fontsize = 18, chromosome = 8)

IGV is irrelevant, you are not changing your BAM file here. If you need more help, please copy past a few rows from your BAM file in order for us to see what it looks like.

ADD REPLY • link 6.5 years ago florian.hahne@novartis.com ★ 1.6k

0

Entering edit mode

BAM file is binary .. and it is really huge I can't load it to r .. however I shared here below the head of my wiggle file which I have generated from the bam files.

ADD REPLY • link 6.5 years ago ta_awwad ▴ 10

0

Entering edit mode

Have you also switched off UCSC chromosome naming convention?

options(ucscChromosomeNames=FALSE)

Otherwise it is for me hard to help without seeing the data and complete code. Can you share it?

ADD REPLY • link 6.5 years ago Robert Ivanek ▴ 730

0

Entering edit mode

yes I did .. still not solved

ADD REPLY • link 6.5 years ago ta_awwad ▴ 10

0

Entering edit mode

this is the head of my wig file:

variableStep chrom=1

3202662 0.00909

3202663 0.00909

3202664 0.00909

3202665 0.00909

3202666 0.00909

3202667 0.00909

3202668 0.00909

3202669 0.00909

3202670 0.00909

3202671 0.00909

3202672 0.00909

3202673 0.00909

3202674 0.00909

3202675 0.00909

3202676 0.00909

3202677 0.00909

3202678 0.00909

3202679 0.00909

3202680 0.00909

3202681 0.00909

3202682 0.00909

3202683 0.00909

ADD REPLY • link 6.5 years ago ta_awwad ▴ 10

0

Entering edit mode

I also tried to generate my gene model using ensemble gif file:

myGenes <- '/gtf/Mus_musculus.GRCm38.90.gtf'

txdb <- makeTxDbFromGFF (myGenes, format= "gtf", dataSource= "Ensembl", organism= "Mus musculus")

txTrack <- GeneRegionTrack(txdb, name="gene", showId=TRUE)

but this error message popped up:

Error in FUN(X[[i]], ...) : Invalid chromosome identifier 'X'
Please consider setting options(ucscChromosomeNames=FALSE) to allow for arbitrary chromosome identifiers

when I put ucscChromosomeNames=F I am still not able to see my track!!!!???

here is the head of my gtf file:

#!genome-build GRCm38.p5

#!genome-version GRCm38

#!genome-date 2012-01

#!genome-build-accession NCBI:GCA_000001635.7

#!genebuild-last-updated 2017-06

1       havana  gene    3073253 3074322 .       +       .       gene_id "ENSMUSG

00000102693"; gene_version "1"; gene_name "4933401J01Rik"; gene_source "havana";

 gene_biotype "TEC";

1       havana  transcript      3073253 3074322 .       +       .       gene_id 

"ENSMUSG00000102693"; gene_version "1"; transcript_id "ENSMUST00000193812"; tran

script_version "1"; gene_name "4933401J01Rik"; gene_source "havana"; gene_biotyp

e "TEC"; transcript_name "4933401J01Rik-201"; transcript_source "havana"; transc

ript_biotype "TEC"; tag "basic"; transcript_support_level "NA";

1       havana  exon    3073253 3074322 .       +       .       gene_id "ENSMUSG

00000102693"; gene_version "1"; transcript_id "ENSMUST00000193812"; transcript_v

ersion "1"; exon_number "1"; gene_name "4933401J01Rik"; gene_source "havana"; ge

ne_biotype "TEC"; transcript_name "4933401J01Rik-201"; transcript_source "havana

"; transcript_biotype "TEC"; exon_id "ENSMUSE00001343744"; exon_version "1"; tag

 "basic"; transcript_support_level "NA";

1       ensembl gene    3102016 3102125 .       +       .       gene_id "ENSMUSG

00000064842"; gene_version "1"; gene_name "Gm26206"; gene_source "ensembl"; gene

_biotype "snRNA";

ADD REPLY • link 6.5 years ago ta_awwad ▴ 10

0

Entering edit mode

You need to set ucscChromosomeNames to FALSE. BTW. what is the meaning of "!!!!???" is this a command?

options(ucscChromosomeNames=FALSE)

ADD REPLY • link 6.5 years ago Robert Ivanek ▴ 730