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

#背景

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

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

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

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

では、そもそも両方の鎖で遺伝子が重なっている領域がどこにどれだけあるか。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)
       element seqnames  start    end width strand tx_id   tx_name
6687   YPR200C   chrXVI 939279 939671   393      -  6663   YPR200C
6688   YPR201W   chrXVI 939922 941136  1215      +  6408   YPR201W
6689   YPR202W   chrXVI 943032 943896   865      +  6409   YPR202W
6690   YPR202W   chrXVI 943880 944188   309      +  6410   YPR203W
6691 YPR204C-A   chrXVI 946856 947338   483      -  6664 YPR204C-A
6692   YPR204W   chrXVI 944603 947701  3099      +  6411   YPR204W

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

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)
GRanges with 528 ranges and 2 metadata columns:
            seqnames           ranges strand   |     tx_id     tx_name
               <Rle>        <IRanges>  <Rle>   | <integer> <character>
    YAL004W     chrI [140760, 141407]      +   |        38     YAL004W
    YAL016W     chrI [124879, 126786]      +   |        32     YAL016W
  YAL019W-A     chrI [114250, 114819]      +   |        29   YAL019W-A
    YAL027W     chrI [ 94687,  95472]      +   |        28     YAL027W
  YAL031W-A     chrI [ 84669,  84977]      +   |        25   YAL031W-A
  YAL034W-A     chrI [ 79718,  80587]      +   |        23   YAL034W-A
    YAL038W     chrI [ 71786,  73288]      +   |        20     YAL038W
    YAL042W     chrI [ 61316,  62563]      +   |        18     YAL042W
  YAL044W-A     chrI [ 57518,  57850]      +   |        17   YAL044W-A
        ...      ...              ...    ... ...       ...         ...
    YPR137W   chrXVI [802359, 804080]      +   |      6371     YPR137W
    YPR143W   chrXVI [818323, 819075]      +   |      6373     YPR143W
    YPR150W   chrXVI [830999, 831520]      +   |      6376     YPR150W
    YPR160W   chrXVI [861306, 864014]      +   |      6384     YPR160W
  YPR170W-B   chrXVI [883239, 883457]      +   |      6390   YPR169W-A
  YPR170W-B   chrXVI [883239, 883595]      +   |      6391   YPR170W-B
    YPR178W   chrXVI [892332, 893729]      +   |      6396     YPR178W
    YPR198W   chrXVI [934034, 935665]      +   |      6407     YPR198W
    YPR204W   chrXVI [944603, 947701]      +   |      6411     YPR204W
  ---
  seqlengths:
      chrI   chrII  chrIII   chrIV    chrV ... chrXIII  chrXIV   chrXV  chrXVI    chrM
    230218  813184  316620 1531933  576874 ...  924431  784333 1091291  948066   85779

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

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

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

result

> subsetByOverlaps(query=sacCer3.cgene.GRange, subject=sacCer3.wgene.GRange,
                   ignore.strand=TRUE)
GRanges with 526 ranges and 2 metadata columns:
            seqnames           ranges strand   |     tx_id     tx_name
               <Rle>        <IRanges>  <Rle>   | <integer> <character>
    YAL005C     chrI [139503, 141431]      -   |       103     YAL005C
  YAL016C-A     chrI [124755, 125069]      -   |        98   YAL016C-A
    YAL020C     chrI [113614, 114615]      -   |        95     YAL020C
    YAL026C     chrI [ 95386,  95823]      -   |        88   YAL026C-A
    YAL031C     chrI [ 84749,  87031]      -   |        86     YAL031C
  YAL034C-B     chrI [ 79489,  79842]      -   |        83   YAL034C-B
  YAL037C-B     chrI [ 72326,  73300]      -   |        80   YAL037C-B
  YAL042C-A     chrI [ 61231,  61608]      -   |        77   YAL042C-A
    YAL045C     chrI [ 57488,  57796]      -   |        74     YAL045C
        ...      ...              ...    ... ...       ...         ...
    YPR130C   chrXVI [793374, 793781]      -   |      6620     YPR130C
    YPR136C   chrXVI [802325, 802837]      -   |      6623     YPR136C
    YPR142C   chrXVI [818130, 818693]      -   |      6629     YPR142C
    YPR151C   chrXVI [831055, 831675]      -   |      6635     YPR151C
  YPR160C-A   chrXVI [862577, 862861]      -   |      6642   YPR160C-A
    YPR170C   chrXVI [882983, 883318]      -   |      6648     YPR170C
    YPR177C   chrXVI [892313, 892684]      -   |      6652     YPR177C
    YPR197C   chrXVI [933898, 934461]      -   |      6661     YPR197C
  YPR204C-A   chrXVI [946856, 947338]      -   |      6664   YPR204C-A
  ---
  seqlengths:
      chrI   chrII  chrIII   chrIV    chrV ... chrXIII  chrXIV   chrXV  chrXVI    chrM
    230218  813184  316620 1531933  576874 ...  924431  784333 1091291  948066   85779