lawremi / ggbio

Grid and ggplot2 based visualization for biological data
111 stars 24 forks source link

autoplot with vcf object #89

Open MichelMoser opened 7 years ago

MichelMoser commented 7 years ago

Hi,

I just started using ggbio in order to plot gene tracks along with ggplots. When i try to add variant information as a new track i have problems creating the track:

Code used:

vcf <- readVcf(file = "INF37-39.SNPs.vcf", genome="Peaxi162")

p2 <- autoplot(vcf[seqnames(vcf)=="Peaxi162Scf00003"], type = "fixed") + xlim(4000, 8000)
Error in `[[<-`(`*tmp*`, name, value = c("C", "C", "C", "G", "C", "G",  : 
  11913 elements in value to replace 11893 elements

So it looks like there is a nonmatching number of elements between whatever it need to replace.

I would like to produce something like in the plot below:

image

which was created using similar code according to http://girke.bioinformatics.ucr.edu/CSHL_RNAseq/mydoc/mydoc_Rgraphics_7/:

library(rtracklayer); library(GenomicFeatures); library(Rsamtools); library(GenomicAlignments); library(VariantAnnotation)
ga <- readGAlignments("./data/SRR064167.fastq.bam", use.names=TRUE, param=ScanBamParam(which=GRanges("Chr5", IRanges(4000, 8000))))
p1 <- autoplot(ga, geom = "rect")
p2 <- autoplot(ga, geom = "line", stat = "coverage")
vcf <- readVcf(file="data/varianttools_gnsap.vcf", genome="ATH1")
p3 <- autoplot(vcf[seqnames(vcf)=="Chr5"], type = "fixed") + xlim(4000, 8000) + theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y=element_blank())
txdb <- makeTxDbFromGFF(file="./data/TAIR10_GFF3_trunc.gff", format="gff3")
p4 <- autoplot(txdb, which=GRanges("Chr5", IRanges(4000, 8000)), names.expr = "gene_id")
tracks(Reads=p1, Coverage=p2, Variant=p3, Transcripts=p4, heights = c(0.3, 0.2, 0.1, 0.35)) + ylab("")

Any hints would be appreciated. Thank you, Michel

lawremi commented 7 years ago

It would help to have access to the file (or just a subset that reprodues the error).