Question: Gviz: AlignmentsTrack from a bam file with non-UCSC chromosome names
gravatar for stianlagstad
2.8 years ago by
stianlagstad90 wrote:

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?

ADD COMMENTlink modified 22 months ago by Johannes Rainer1.3k • written 2.8 years ago by stianlagstad90
gravatar for Robert Ivanek
2.8 years ago by
Robert Ivanek530
Robert Ivanek530 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))
ADD COMMENTlink written 2.8 years ago by Robert Ivanek530

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

ADD REPLYlink written 2.8 years ago by stianlagstad90

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.

ADD REPLYlink written 2.8 years ago by stianlagstad90

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),`. I was wondering if there was a similar work-around for gene model tracks? 

ADD REPLYlink written 22 months ago by vladislav.kim10
gravatar for Johannes Rainer
22 months ago by
Johannes Rainer1.3k
Johannes Rainer1.3k 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

ADD COMMENTlink written 22 months ago by Johannes Rainer1.3k

Thank you! This works

ADD REPLYlink written 22 months ago by vladislav.kim10
Please log in to add an answer.


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