Search
Question: Givz Doesn't plot coverage genetracks that are aligned against ensemble ref
0
gravatar for ta_awwad
29 days ago by
ta_awwad0
ta_awwad0 wrote:

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

ADD COMMENTlink modified 29 days ago • written 29 days ago by ta_awwad0
0
gravatar for florian.hahne@novartis.com
29 days ago by
Switzerland

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

ADD COMMENTlink written 29 days ago by florian.hahne@novartis.com1.5k
0
gravatar for ta_awwad
29 days ago by
ta_awwad0
ta_awwad0 wrote:

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 COMMENTlink modified 29 days ago • written 29 days ago by ta_awwad0

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 REPLYlink written 29 days ago by Robert Ivanek350

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 REPLYlink written 28 days ago by ta_awwad0

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 REPLYlink written 28 days ago by florian.hahne@novartis.com1.5k

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 REPLYlink written 28 days ago by ta_awwad0

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 REPLYlink written 28 days ago by Robert Ivanek350

yes I did .. still not solved

ADD REPLYlink written 28 days ago by ta_awwad0

 

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 REPLYlink written 28 days ago by ta_awwad0

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 REPLYlink modified 28 days ago • written 28 days ago by ta_awwad0

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

options(ucscChromosomeNames=FALSE)
ADD REPLYlink written 28 days ago by Robert Ivanek350
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 299 users visited in the last hour