thr3ads.net - R help - [R] Whole genome searching of 100bp "D" sequence [Apr 2011]

If this information is useful, please help other people find it:
Share via:

Amos Folarin

2011-Apr-15 18:53 UTC

[R] Whole genome searching of 100bp "D" sequence

Hi,
I was wondering I'm going about this in the correct way. I need to test if
there are coding sequences or exons in hg19 which match a string of 100bp
"D" i.e. [A,G or T]. However I'm getting a strange result.

I get a hit on chr7, using the 100bp search however when I search with 60bp
sequence of "D" I don't get any hits.


library("BSgenome")
library("Biostrings")
library("BSgenome.Hsapiens.UCSC.hg19")
library("biomaRt")
library("GenomicFeatures")




#extract the alignments which match real genes
#(add this onto stefans script … req to buid the gr object using the regions
from the BSgenome alignment of the 100bp seqs.)

txdb <- makeTranscriptDbFromUCSC(genome = "hg19", tablename =
"knownGene")
## do once locally & save



#query.plus <-
DNAString("DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD")
# 100bp C free sequence
query.plus <-
DNAString("DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD")
#
60bp C free sequence
query.minus <- reverseComplement(query.plus)


chrList <- c("chr1", "chr2", "chr3",
"chr4", "chr5", "chr6", "chr7",
"chr8",
"chr9", "chr10", "chr11", "chr12",
"chr13", "chr14", "chr15", "chr16",
"chr17", "chr18", "chr19", "chr20",
"chr21", "chr22", "chrX", "chrY")

#access/group the by subset of annotation
annotGr <- exons(txdb)
#annotGr <- cds(txdb)

wholeGenomeMatch <- function()
{
    for(i in chrList)
    {
    #matches on plus strand
    matches.plus.strand <- matchPattern(query.plus, Hsapiens[[i]],
fixed=FALSE, max.mismatch=5)
    mp <- as.matrix(matches.plus.strand)

    #matches on minus strand
    matches.minus.strand <- matchPattern(query.minus, Hsapiens[[i]],
fixed=FALSE, max.mismatch=5)
    mm <- as.matrix(matches.minus.strand)

    OL.p <- NULL
    if(nrow(mp) > 0)
    {
    ##### MATCH THE POSITIVE STRAND ######
    #need to get start and end positions (end = start + (length-1))
    gr <- GRanges(
    seqnames = rep(i,nrow(mp)),
    ranges = IRanges(start = mp[[1]],
    end = mp[[1]]+(mp[,2]-1)),
    strand = rep("+", nrow(mp)))

    #OL <- findOverlaps(query=gr, subject=annotGr)
    OL.p <- annotGr[(!is.na(match(annotGr, gr)))]
    #as.data.frame(OL) ## view the results

    #tdata.p <- annotGr[unique(queryHits(OL)),]
    #tdata.p <- as.data.frame(tdata.p, row.names = NULL, optional = FALSE)

    #write.table(tdata.p, paste(i,"plus.txt"))
    cat( paste(i,"plus.txt\n"))
    }

    OL.m <- NULL
    if(nrow(mm) > 0)
    {
    ##### MATCH THE MINUS STRAND ########
    #need to get start and end positions (end = start + (length-1))
    gr <- GRanges(
    seqnames = rep(i,nrow(mm)),
    ranges = IRanges(start = mm[[1]],
    end = mm[[1]]+(mm[,2]-1)),
    strand = rep("-", nrow(mm)))

    #OL <- findOverlaps(query=gr, subject=annotGr)
    OL.m <- annotGr[(!is.na(match(annotGr, gr)))]
    #tdata.m <- annotGr[unique(queryHits(OL)),]
    #tdata.m <- as.data.frame(tdata.m, row.names = NULL, optional = FALSE)

    #write.table(tdata.m, paste(i,"minus.txt"))

    cat( paste(i,"minus.txt\n"))

    #write.table(rbind(tdata.p, tdata.p), paste(i, ".txt"))
    }

    ##write all matching results into one file
    write.table(rbind(as.data.frame(OL.p), as.data.frame(OL.m)),
file="60bp_v_exons_max-mismatch=5.txt", append=TRUE)

    }

}

wholeGenomeMatch()




#in the first instance with the 100bp search I get 1 hit in the output file:
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"1" "chr7" 420815 422845 2031 "+" 97277
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"

#with the 60 bp query string my output file reads no hits
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"
"seqnames" "start" "end" "width"
"strand" "exon_id"



Your help is much appreciated,
Amos

	[[alternative HTML version deleted]]

Ben Bolker

2011-Apr-15 22:14 UTC

head link

[R] Whole genome searching of 100bp "D" sequence

Amos Folarin <amosfolarin <at> gmail.com> writes:
> I was wondering I'm going about this in the correct way. I need to test
if
> there are coding sequences or exons in hg19 which match a string of 100bp
> "D" i.e. [A,G or T]. However I'm getting a strange result.
> 
> I get a hit on chr7, using the 100bp search however when I search with 60bp
> sequence of "D" I don't get any hits.
> 
  Try posting this to the Bioconductor e-mail list, where you will
have a much better likelihood of getting useful advice. Only a few
hard-core bioinformatics people on this list.  (Do mention that you
posted here first.)

  Ben Bolker

Reasonably Related Threads

Search for more apparently analagous threads

R help - Apr 2011 - Whole genome searching of 100bp "D" sequence

[R] Whole genome searching of 100bp "D" sequence

[R] Whole genome searching of 100bp "D" sequence

Reasonably Related Threads