Question

Get Intron Seq

0

Entering edit mode

Lillian • 0

@c209ec01

Last seen 18 months ago

United States

Hello,

I'm wondering how to intron sequence by coordinates.

For example, I know that there might be an intron in chr20 start=50941210,end=50942030,

I'm wondering if I could use the BSgenome package in R to derive the DNAString object to get the sequence of intron, wondering if it is possible.

Thank you so much

intron BSgenome.Hsapiens.UCSC.hg38 BSgenome • 1.0k views

ADD COMMENT • link updated 18 months ago by ATpoint ★ 4.6k • written 18 months ago by Lillian • 0

score 0 · Answer 1 · 2023-06-20

0

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 2 days ago

United States

Yes, quite easily in fact.

> library(BSgenome.Hsapiens.UCSC.hg38)
> getSeq(Hsapiens, GRanges("chr20:50941210-50942030"))
DNAStringSet object of length 1:
    width seq
[1]   821 CTACAAAAAAAAAAGTGAATATATATATATATAT...TACTAATAAAGAAGTTGTAACACATGTACCTAC

You could wrap that in as.character if you want the entire sequence, or write to a FASTA file using writeXStringSet.

ADD COMMENT • link 18 months ago James W. MacDonald 67k

0

Entering edit mode

Hi!

Thank you so much for providing this.

I just have a question about it that I have tried to create a Granges object with the same coordinates as yours, however, why I feel like the sequences were different?

> gr <- GRanges(gr)
> gr
GRanges object with 1 range and 3 metadata columns:
    seqnames            ranges strand |            ENSG      pval       dPSI
       <Rle>         <IRanges>  <Rle> |     <character> <numeric>  <numeric>
  2    chr20 50941210-50942030      - | ENSG00000000419   0.17493 -0.0109049
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

> getSeq(BSgenome.Hsapiens.UCSC.hg38, gr)
DNAStringSet object of length 1:
    width seq                                                                                                       names               
[1]   821 GTAGGTACATGTGTTACAACTTCTTTATTAGTATAACTACTGGAAATGTAT...TATATATATATATATATATATATATATATATATTCACTTTTTTTTTTGTAG 2



 > getSeq(BSgenome.Hsapiens.UCSC.hg38, GRanges("chr20:50941210-50942030"))
DNAStringSet object of length 1:
    width seq
[1]   821 CTACAAAAAAAAAAGTGAATATATATATATATATATATATATATATATATATATAAATAAAA...CTTGCTATGAATACATTTCCAGTAGTTATACTAATAAAGAAGTTGTAACACATGTACCTAC

I just feel a bit confused because I use the same coordinates and BSgenome.

Thank you so much

ADD REPLY • link 18 months ago Lillian • 0

0

Entering edit mode

Your gr is stranded, the sequence is the reverse complement of James version due to strand-aware extraction. Unstranded defaults to positive strand.

ADD REPLY • link 18 months ago ATpoint ★ 4.6k