VOOZH about

URL: https://qiita.com/wakuteka/items/10027edccc6c2e244cd2

⇱ DNA二重らせんの両側にそれぞれ別の遺伝子が存在する箇所をRで抽出する #Bioconductor - Qiita

👁 Image
8

Go to list of users who liked

6

Share on X(Twitter)

Share on Facebook

Add to Hatena Bookmark

More than 5 years have passed since last update.

@wakutekain👁 Image
株式会社LIFULL

DNA二重らせんの両側にそれぞれ別の遺伝子が存在する箇所をRで抽出する

8
Last updated at Posted at 2013-01-16

背景

DNAは二重らせんであり、ATGCの塩基からなる二つの鎖がたがいに相補鎖となっている。（これらを＋鎖と－鎖と呼んだり、二重らせん構造の発見者の名をとってWatson鎖とCrick鎖と呼んだりする。）
遺伝子はこの二つの鎖のどちらか片側にコードされているが、まれに同じ塩基位置の両方の鎖に異なる遺伝子が重なってコードされていることがある。例えばRNA-Seqによって遺伝子の発現量を網羅的に知りたい場合、このような箇所では貼り付いたリードがどちらの遺伝子から得られたものかわからない。

そんな時、鎖特異的なRNA-Seqを行えば、リードが貼り付いた位置で両方の鎖に遺伝子があってもどちらの遺伝子から得られたリードであるか推測することができる。

下の図の中段は鎖特異的に得られたリードのマッピング結果を表していて、青いリードは＋鎖にコード（右向きに表示：→）されている遺伝子に、赤いリードは－鎖にコード（左向きに表示：←）されている遺伝子に貼り付いている。

👁 YAL004W

準備（パッケージ、サンプルデータの読み込み）

では、そもそも両方の鎖で遺伝子が重なっている領域がどこにどれだけあるか。Bioconductorの
GenomicRangesパッケージを使って、それを調べてみる。遺伝子リストは出芽酵母（Saccharomyces cerevisiae）のものを用いた。

code
# パッケージの読み込みlibrary(GenomicRanges)library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)# 遺伝子名とその開始位置・終了位置などのリストをGRangeList形式で読み込むsacCer3.gene.GRangesList<-transcriptsBy(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene,by="gene")

ここで得られた遺伝子リストはGRangeList形式になっているので、中身をデータフレーム形式で閲覧したいときはas()関数を使う。

code
sacCer3.gene.dataframe<-as(sacCer3.gene.GRangesList,"data.frame")
result
>tail(sacCer3.gene.dataframe)elementseqnamesstartendwidthstrandtx_idtx_name6687YPR200CchrXVI939279939671393-6663YPR200C6688YPR201WchrXVI9399229411361215+6408YPR201W6689YPR202WchrXVI943032943896865+6409YPR202W6690YPR202WchrXVI943880944188309+6410YPR203W6691YPR204C-AchrXVI946856947338483-6664YPR204C-A6692YPR204WchrXVI9446039477013099+6411YPR204W

遺伝子リストを＋鎖、－鎖に分割して重なりを見る

GRangeList形式で遺伝子リストが読み込めたら、unlistした後で個々のstrandごとに分解する。

＋,－それぞれのGRangeオブジェクトを作成する
sacCer3.gene.GRange<-unlist(sacCer3.gene.GRangesList)sacCer3.wgene.GRange<-sacCer3.gene.GRange[strand(sacCer3.gene.GRange)=="+",]sacCer3.cgene.GRange<-sacCer3.gene.GRange[strand(sacCer3.gene.GRange)=="-",]

ここで、subsetByOverlaps(query=..., subject=...)関数にignore.strand=TRUEオプションを与えて使うと、queryリスト中からsubjectで与えたリストと領域が重なっている遺伝子が抽出できる。

result
>subsetByOverlaps(query=sacCer3.wgene.GRange,subject=sacCer3.cgene.GRange,ignore.strand=TRUE)GRangeswith528rangesand2metadatacolumns:seqnamesrangesstrand|tx_idtx_name<Rle><IRanges><Rle>|<integer><character>YAL004WchrI[140760,141407]+|38YAL004WYAL016WchrI[124879,126786]+|32YAL016WYAL019W-AchrI[114250,114819]+|29YAL019W-AYAL027WchrI[94687,95472]+|28YAL027WYAL031W-AchrI[84669,84977]+|25YAL031W-AYAL034W-AchrI[79718,80587]+|23YAL034W-AYAL038WchrI[71786,73288]+|20YAL038WYAL042WchrI[61316,62563]+|18YAL042WYAL044W-AchrI[57518,57850]+|17YAL044W-A.....................YPR137WchrXVI[802359,804080]+|6371YPR137WYPR143WchrXVI[818323,819075]+|6373YPR143WYPR150WchrXVI[830999,831520]+|6376YPR150WYPR160WchrXVI[861306,864014]+|6384YPR160WYPR170W-BchrXVI[883239,883457]+|6390YPR169W-AYPR170W-BchrXVI[883239,883595]+|6391YPR170W-BYPR178WchrXVI[892332,893729]+|6396YPR178WYPR198WchrXVI[934034,935665]+|6407YPR198WYPR204WchrXVI[944603,947701]+|6411YPR204W---seqlengths:chrIchrIIchrIIIchrIVchrV...chrXIIIchrXIVchrXVchrXVIchrM2302188131843166201531933576874...924431784333109129194806685779

リストの一番最初にある YAL004W 遺伝子の周辺をIGVで可視化したのが最初にあげた図である。

👁 YAL004W

これを見ると、 YAL004W 遺伝子は＋鎖にコード（左向きに表示）されているのに、青いリードがほとんど貼り付いておらず、この位置にあるリードの多くは－鎖にコードされた YAL005C 遺伝子由来であると推測できる。

なお、遺伝子領域が重なっている「－鎖の」遺伝子を表示させたいときは、さきほどのsubsetbyOverlaps()関数のqueryとsubjectを逆にしてやれば良い。

result
>subsetByOverlaps(query=sacCer3.cgene.GRange,subject=sacCer3.wgene.GRange,ignore.strand=TRUE)GRangeswith526rangesand2metadatacolumns:seqnamesrangesstrand|tx_idtx_name<Rle><IRanges><Rle>|<integer><character>YAL005CchrI[139503,141431]-|103YAL005CYAL016C-AchrI[124755,125069]-|98YAL016C-AYAL020CchrI[113614,114615]-|95YAL020CYAL026CchrI[95386,95823]-|88YAL026C-AYAL031CchrI[84749,87031]-|86YAL031CYAL034C-BchrI[79489,79842]-|83YAL034C-BYAL037C-BchrI[72326,73300]-|80YAL037C-BYAL042C-AchrI[61231,61608]-|77YAL042C-AYAL045CchrI[57488,57796]-|74YAL045C.....................YPR130CchrXVI[793374,793781]-|6620YPR130CYPR136CchrXVI[802325,802837]-|6623YPR136CYPR142CchrXVI[818130,818693]-|6629YPR142CYPR151CchrXVI[831055,831675]-|6635YPR151CYPR160C-AchrXVI[862577,862861]-|6642YPR160C-AYPR170CchrXVI[882983,883318]-|6648YPR170CYPR177CchrXVI[892313,892684]-|6652YPR177CYPR197CchrXVI[933898,934461]-|6661YPR197CYPR204C-AchrXVI[946856,947338]-|6664YPR204C-A---seqlengths:chrIchrIIchrIIIchrIVchrV...chrXIIIchrXIVchrXVchrXVIchrM2302188131843166201531933576874...924431784333109129194806685779
8

Go to list of users who liked

6
0

Go to list of comments

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
8

Go to list of users who liked

6