#### The support.bioconductor.org editor has been updated to markdown! Please see more info at: Tutorial: Updated Support Site Editor

Question: Gviz: AlignmentsTrack from a bam file with non-UCSC chromosome names
1
3.0 years ago by

I have some bam files with chromosome names like "1" instead of "chr1". I run into problems when I want to plot data from ensembl along with coverage data from these bam files. The solution this far has been to use the command below to edit the bamfile before using it in my scripts:

samtools view -H test.bam | sed -e 's/SN:$$[0-9XY]$$/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - test.bam > test_chr.bam


I don't feel good doing that, so I would like to modify the function that reads data from the bam file somehow so that it returns the chromosome name prefixed with "chr". Is that possible? I tried editing the default import function, only editing the line giving seqnames, like this:

seqnames = if (is.null(reads$rname)) character() else paste("chr", reads$rname, sep="")

But that didn't lead me anywhere. Does anyone have suggestions?

gviz • 1.0k views
modified 2.1 years ago by Johannes Rainer1.4k • written 3.0 years ago by stianlagstad90
Answer: Gviz: AlignmentsTrack from a bam file with non-UCSC chromosome names
2
3.0 years ago by
Robert Ivanek590
Switzerland
Robert Ivanek590 wrote:

I think you are on the right track. The .import.bam.alignments function has a selection parameter (Granges object) which is constructed from the plotTracks arguments and limits the reads which are read from the BAM file. I assume your other tracks are using the "chr" in sequence names and therefore this selection object will point to "wrong" chromosome. Simple fix for that would be to add one additional line into the modified .import.bam.alignments function. Add this line after the  BamFile command (before creating the ScanBamParam object). And keep the modification you made before.

seqlevels(selection) <- sub("chr", "", seqlevels(selection))

That worked well, thank you! I guess I was too tired to consider how I could edit the selection at the time of posting :)

Is it possible to give a warning when the AlignmentsTrack doesn't find anything in the bamfile? Right now I have a boolean "nonUCSC" variable in my plotting function that I have to set to true if I know the bamfile has "1" and not "chr1" as seqname, to make my plotting function to use the alternative import function. I would like to notify the user that he/she can use this parameter if the AlignmentsTrack doesn't find anything.

I am facing a similar issue: my BAM files are also in non-UCSC format. Plotting the AlignmentTrack alone works, but combining it with the BiomartGeneRegionTrack fails because chromosome names don't match. Is there an easy way to change the chromosome names of the BiomartGeneRegionTrack object? There is a way to do this for ideograms: levels(itrack@bandTable$chrom) <- sub("^chr", "", levels(itrack@bandTable$chrom), ignore.case=T). I was wondering if there was a similar work-around for gene model tracks?

Answer: Gviz: AlignmentsTrack from a bam file with non-UCSC chromosome names
0
2.1 years ago by
Johannes Rainer1.4k
Italy
Johannes Rainer1.4k wrote:

Just another suggestion - if you have Ensembl based data, why not use ensembldb? You could use the dedicated getGeneRegionTrackForGviz to extract GRanges to pass to Gviz GeneRegionTrack method. Don't forget then to call options(ucscChromosomeNames = FALSE) too so that Gviz is not crying about the chromosome names not starting with "chr". Check the ensembldb vignette for more details and description.

Also, using EnsDb databases, since they match certain Ensembl releases, would allow you to ensure that the annotation you used for alignment or counting matches the annotation to plot the data.

cheers, jo