Entering edit mode
Prashantha Hebbar
▴
260
@prashantha-hebbar-3526
Last seen 5.0 years ago
Dear Friends,
My following geneRange object has all chromosomal information in one
GRange object.
library(VariantAnnotation)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
txg = transcriptsBy(txdb,by="gene")
> geneRange <- sort(unlist(range(txg[as.character(myGenes$geneId)])))
> geneRange GRanges with 11 ranges and 0 metadata columns: seqnames
ranges strand <rle> <iranges> <rle> 4067 chr8 [
56792386, 56925006] + 26579 chr11 [ 69061622, 69156450]
+ 595 chr11 [ 69455873, 69469242] + 2017 chr11 [ 70244612,
70282690] + 10413 chr11 [101981192, 102104154] + 330
chr11 [102188181, 102210135] + 329 chr11 [102217966,
102249401] + 9965 chr11 [ 69513006, 69519106] - 2249
chr11 [ 69587797, 69590171] - 2248 chr11 [ 69624736,
69634192] - 47 chr17 [ 40023179, 40075272] - ---
I have to get SNP and genotype information for the same from 1000
genome. Below code does that. However, each time I have to create ftp
url and connect to 1000 Genome resource.
thats painstaking.
i=0
for(chr in as.vector(seqnames(geneRange))){
baseUrl <- paste("ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/phase1/
analysis_results/integrated_call_sets/ALL.chr",chr,".integrated_phase1
_v3.20101123.snps_indels_svs.genotypes.vcf.gz",sep="")
i=i+1
print(baseUrl)
inGene <- snpsInGene(geneRange[i],baseUrl,geno="GT",info=NA)
print(inGene)
}
If I have my geneRange information as per chromosome in GRangeList
form, my script becomes more efficient. because, I can retrive
genotypes from each chromosomes in one go.
So, May I know a way to convert geneRange GRange object to GRangeList
Object as per chromosome number?
Thanking you in anticipation.
Regards,
Prashantha
[[alternative HTML version deleted]]