How to draw the chromosome specific coverage plot using GFF3 file

Question

[ajaybioinfo]

Dear all,

Good evening, I am working on the genome sequencing of a plant, now my guide want a make an image of gene coverage, repeat coverage and SNP coverage specific to chromosomes in a single picture. for your reference please see this link http://www.nature.com/ng/journal/v48/n6/fig_tab/ng.3565_F1.html. my guide wants that I should make picture like this.

I have genome in chromosomes in FASTA format, gene prediction result in GFF3 format, SNP result in vcf format and repeat prediction result in GFF3 format.

please guide me and tell how I can make image like this, As I am struggling since 12 days without any success.

Thanks in advance and waiting for your guidance and suggestion.

Warm Regards Ajay

[lakhujanivijay]

What all have you tried in the past 18 days?

[ajaybioinfo]

Dear Vijay I have tried in R to generate required image but could not understand how I can generate the coverage graph,.. have you any idea how I can generate those plot if yes please tell.

Thanks in Advance Ajay

Answer 1

[bernatgel]

Hi @ajaybioinfo,

It won't be exactly the same image, but you can get something similar with the R package karyoploteR.

1 - You can create a genome.bed file with the length of your chromosomes and use it to create a plot as in the Custom Genome tutorial. If you want to plot only one chromosome at a time, use the parameter "chromosomes".

2 - You can load the GFF's with the import function from rtracklayer as shown here and then use the kpPlotDensity function to plot the density along the genome.

3 - You can the load the VCF file with the readVcf function from VariantAnnotation, extract a GRanges with the position of the variants with rowRanges and finally use kpPlotDensity again to plot the SNP density.

Hope this helps

Answer 2

[ajaybioinfo]

Thanks sir for your suggestion

I read your tutorial and able to draw the " ideogram using a custom genomes tutorial using below code

---

library(karyoploteR) custom.genome <- toGRanges(data.frame(chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20"), start=c(1, 1), end=c(15232102, 13796269, 16490533, 21168693, 11781399, 11560298, 11574172, 11330313, 18393972, 10058126, 14275844, 12202370, 15640443, 9516842, 10244494, 18252290, 10284815, 15009555, 23318555, 13640176))) custom.genome <- toGRanges("mygenome1.txt") AK <- plotKaryotype(genome=custom.genome)

---

But I dont understand how I can draw the coverage plot of the Gene and Repeat. Both the data is in the below format

chr start end name annotation 1 2345 3862 gene1 LRR1 1 5678 5764 gene2 LRR2 1 9011 7666 gene3 LRR3 1 12344 9568 gene4 LRR4 1 15677 11470 gene5 LRR5 1 19010 13372 gene6 LRR6 1 22343 15274 gene7 LRR7

Please guide me how I can draw the gene density plot which you have shown in your tutorial

Please help me as we are in the final phase of manuscript material preparation ,I hope you can help student like me as I am struggling since 6 days with no success.

Waiting for your kind reply and suggestion. I would be a great help if you can tell the steps along with some reference code so that I can draw this Image.

Thanking You

[bernatgel]

Hi @ajaybioinfo,

I've replied to the email you sent me.

Here's the reply

If you want to add some tracks with density data in the format you have given, you will only need to load that data into R and plot it with kpPlotDensity.

Please follow these steps:

1 - Load all your data into R. You can use toGRanges function to directly read the files if they are in the order you specified (1st column chr, 2nd start, 3rd end, and then others optionally). If your data file is called "genes.txt" do it with

my.genes <- toGRanges("genes.txt")

This will read the data and create a GRanges object as needed. If you have any problem you can read the regioneR vignette

2 - Once you have loaded all your data, you will have to plot it using kpPlotDensity. You can find more information on this function the tutorial page. If you have, for example, 2 GRanges objects: "my.genes" and "my.repeats" you can plot them with:

#First prepare your custom genome as you have already done
#Create the karyoplot
AK <- plotKaryotype(genome=custom.genome)

#Plot the genes (on the bottom half of the plotting region as specified by r0 and r1)
kpPlotDensity(AK, data=my.genes, r0=0, r1=0.5)

#Then plot the repeats on above the genes
kpPlotDensity(AK, data=my.repeats, r0=0.5, r1=1)

In the kpPlotDensity Tutorial page and in the Gene Density example you can find more information on how to customize your density plot (change the window width for density calculation), or change the colors, etc...

In the Data Positioning page you can find more information on the r0 and r1 argument and how to use them to change the positioning of your data.

[ajaybioinfo]

Respected Sir,

Thanks for your reply and suggestion, as per your suggestion I have run the commands and successfully loaded the gene file whose density I need to plot as per your suggested command (my.genes <- toGRanges("myband.txt")

But When I tried to generate the plot using command (kpPlotDensity(AK, data=my.gene, r0=0, r1=0.5)

I am getting below Error:

Error in tileGenome(stats::setNames(karyoplot$chromosome.lengths, karyoplot$chromosomes), : 'seqlengths' has names that are NA or the empty string.

now I am again stuck, how to solve this error,

Command used by me:

---

library(karyoploteR) library(GenomeInfoDb) library(GenomicRanges) custom.genome <- toGRanges(data.frame(chr=c("1"), start=c(1, 1), end=c(15232102))) AK <- plotKaryotype(genome=custom.genome) my.genes <- toGRanges("myband.txt") kpPlotDensity(AK, data=my.gene, r0=0, r1=0.5)

---

I am requesting you, please suggest me correction and command so that I can draw this image and finished my assigned task.

waiting for your valuable suggestions Ajay

[bernatgel]

Hi @ajaybioinfo,

In your code, when you create the custom genome, you have two starts and only one chromosome. This creates an invalid genome (with two chromosomes of the same name) and prevent kpPlotDensity from working.

This should work

library(karyoploteR) 
custom.genome <- toGRanges(data.frame(chr=c("1"), start=c(1), end=c(15232102)))
AK <- plotKaryotype(genome=custom.genome) 
my.genes <- toGRanges("myband.txt")
my.genes <- createRandomRegions(nregions = 1000, genome=custom.genome)
kpPlotDensity(AK, data=my.genes, r0=0, r1=0.5)

Hope this helps and sorry for the late response.