Question: How can I create a GmapGenome that doesn't start at index 1?
0
4 months ago by
Kyle Johnsen20
United States/Brigham Young University
Kyle Johnsen20 wrote:

Let me explain: I want to use a small "reference genome" of only a gene, for use in examples. But I need it to still have the coordinates that correspond to the actual genome coordinates. But when I create a GmapGenome using this truncated fasta, it begins at index 1, which prevents me from using it as wanted.

The header of the fasta file: >18 dna:chromosome chromosome:GRCz11:18:5213338:5227420:1

This is how I build the genome: refGenome <- gmapR::GmapGenome('slc24a5.fa', create=TRUE)

I want to be able to reference the genome using coordinates between 5213338 and 5227420, instead of 1 and 14083 (the length of the gene). Is this possible?

To clarify, what I need to do in the end is call variants (with VariantTools) using BAM files (aligned to genomic coordinates).

gmapr • 97 views
modified 4 months ago • written 4 months ago by Kyle Johnsen20
Answer: How can I create a GmapGenome that doesn't start at index 1?
1
4 months ago by
United States
Michael Lawrence11k wrote:

You'll have to map the resulting coordinates using e.g. mapToTranscripts(). See the GenomicAlignments package.

Ok, in my case, I want to call variants in this region, but the BAM files are in genomic coordinates. If I understand right, mapToTranscripts() would help if I had the BAM files in transcript coordinates, called variants, and then converted to genomic coordinates. Is that right, and if so, is that the only way?

Sorry, I meant mapFromTranscripts. You could call variants in your restricted genome space and then map to the true genome using that function.