Hello..I have a list of .CEL files from HTA 2.0 chips. I created an ExpressionSet object by reading them with the R package "oligo". This object has data on 70523 features (genes) and 48 samples. What is best way to annotate these genes so that I can create a "GenomicRanges" object with the following attributes: Chromosome, Transcript_start, Transcript end ?
> library(pd.hta.2.0)
> load(paste0(path.package("pd.hta.2.0"), "/extdata/netaffxTranscript.rda"))
> annot <- pData(netaffxTranscript)
> annot <- annot[!annot$seqname %in% c("", NA, "NA"),]
> annot$strand[annot$strand == ""] <- "*"
> pd.hta.gr <- GRanges(annot$seqname, IRanges(annot$start, annot$stop), annot$strand)
EDIT: added names
> namespd.hta.gr) <- annot$transcriptclusterid
> pd.hta.gr
GRanges object with 67758 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
TC01000001.hg.1 chr1 [ 11869, 14409] +
TC01000002.hg.1 chr1 [ 29554, 31109] +
TC01000003.hg.1 chr1 [ 69091, 70008] +
TC01000004.hg.1 chr1 [160446, 161525] +
TC01000005.hg.1 chr1 [317811, 328581] +
... ... ... ...
TCUn_gl000212000002.hg.0 chrUn_gl000212 [ 65184, 66597] +
TCUn_gl000219000001.hg.0 chrUn_gl000219 [ 45397, 57412] +
TCUn_gl000222000004.hg.0 chrUn_gl000222 [122791, 122814] -
TCUn_gl000223000002.hg.0 chrUn_gl000223 [ 351, 377] -
TCUn_gl000229000001.hg.0 chrUn_gl000229 [ 16405, 16424] +
-------
seqinfo: 60 sequences from an unspecified genome; no seqlengths
Sincere thanks for your reply!…But this gives the annotation for 67758 transcripts, while I have the expression data on 70523 transcripts. How can I get the annotation for the remaining transcripts?
Thank you very much!
Hello James,
I have a similar situation. I have read all my 30 CEL file into R and Normalized them using Oligo package. How can I use pd.hta.2.0 package to annotate my data ? I have been searching for codes. From what I have seen some people have tried to build their on library but I am not sure if is that a necessary step to annotate my data?
thanks in advance.
Best,
Hoss
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
The short answer is that you can't. If you re-read the code in my previous answer, you will note this line:
> annot <- annot[!annot$seqname %in% c("", NA, "NA"),]which removes all the data where the seqname (chromosome) is missing. It turns out there are 2995 probesets for which there is no chromosome in the annotation file. If you don't have a chromosome, how will you add it to the GRanges object? If we look a bit closer at these probesets, they don't have start/end information either:
> annot.noseqname <- annot[annot$seqname == "", 1:6] > head(annot.noseqname) transcriptclusterid probesetid seqname strand start stop 2824546_st 2824546_st 2824546_st 0 0 2824549_st 2824549_st 2824549_st 0 0 2824551_st 2824551_st 2824551_st 0 0 2824554_st 2824554_st 2824554_st 0 0 2827992_st 2827992_st 2827992_st 0 0 2827995_st 2827995_st 2827995_st 0 0 > table(annot.noseqname$start) 0 2995 > table(annot.noseqname$stop) 0 2995And if we go to netaffx to see what these probesets are supposed to measure, we either return no results, or we get redirected to some random page. Try querying for 2824546_st at netaffx, or better yet, one of these PSR probesets:
> tail(annot.noseqname) transcriptclusterid probesetid PSR6_ssto_hap7002446.hg.1 PSR6_ssto_hap7002446.hg.1 PSR6_ssto_hap7002447.hg.1 PSR6_ssto_hap7002447.hg.1 PSR6_ssto_hap7002454.hg.1 PSR6_ssto_hap7002454.hg.1 PSR6_ssto_hap7002455.hg.1 PSR6_ssto_hap7002455.hg.1 PSR6_ssto_hap7002456.hg.1 PSR6_ssto_hap7002456.hg.1 PSR6_ssto_hap7002459.hg.1 PSR6_ssto_hap7002459.hg.1 seqname strand start stop PSR6_ssto_hap7002446.hg.1 0 0 PSR6_ssto_hap7002447.hg.1 0 0 PSR6_ssto_hap7002454.hg.1 0 0 PSR6_ssto_hap7002455.hg.1 0 0 PSR6_ssto_hap7002456.hg.1 0 0 PSR6_ssto_hap7002459.hg.1 0 0So long story short, you cannot get the position information for all of the probesets, because some of them have no such information.
The 2995 probesets in question are various types of controls, and are not considered main content for gene level analysis.