1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
|
### =========================================================================
### getTranscriptSeqs methods
### =========================================================================
setMethod("getTranscriptSeqs", c("GRangesList", "BSgenome"),
function(query, subject, ...)
{
extractTranscriptSeqs(subject, query)
}
)
setMethod("getTranscriptSeqs", c("GRangesList", "FaFile"),
function(query, subject, ...)
{
## order by exon rank, check for mixed chromosomes, mixed strand
txlist <- GenomicFeatures:::.makeUCSCTxListFromGRangesList(query,
decreasing.rank.on.minus.strand=FALSE)
## back to GRanges for getSeq method
widthEL <- elementLengths(query)
gr <- GRanges(
seqnames=Rle(txlist$chrom, widthEL),
ranges=IRanges(start=unlist(txlist$exonStarts),
end=unlist(txlist$exonEnds)),
strand=Rle(txlist$strand, widthEL))
seqs <- callGeneric(gr, subject, ...)
oneseq <- unlist(seqs)
widthView <- sum(width(query))
vws <- successiveViews(oneseq, widthView)
txseq <- as(vws, "DNAStringSet")
names(txseq) <- names(query)
txseq
}
)
setMethod("getTranscriptSeqs", c("GRanges", "FaFile"),
function(query, subject, ...)
{
fa <- open(subject)
idx <- scanFaIndex(fa)
seqmatch <- seqlevels(query) %in% seqlevels(idx)
if (any(seqmatch == FALSE)){
missingSeq <- seqlevels(subject)[seqmatch == FALSE]
stop(paste("sequence ", missingSeq,
" not found in 'query' \n", sep=""))
}
dna <- getSeq(fa, query)
close(fa)
dna
}
)
|
