Consensus CDSの領域をGRangesオブジェクトにする

Bioconductor - TxDb.Hsapiens.UCSC.hg19.knownGeneは便利なのですが、ある遺伝子に対応するtranscriptが多すぎてどうしたものかと困ってました。で、UCSCのゲノムブラウザで見るとConsensus CDS(CCDS)のアノテーションもあるのですが、このBioconductorのパッケージからそれを取得する方法が分かりませんでした。正確にはmakeTranscriptDbFromUCSC(genome="hg19", tablename="ccdsGene")などとしてCCDSの領域は取得はできるんですが、他のアノテーション(Gene IDなど)が一切消えるので全く使いものにならない感じです。

仕方がないのでCCDS ProjectのFTPサーバから元データを引っ張ってきてそれをパースしてアノテーションのあるCCDSのデータを作ったのでちょろっと共有しておきます。

まずはCCDS ProjectのFTPサーバからCCDSの領域とアノテーションを収めたファイルを取得します (ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/)。 reference genomeのビルドによっていくつかありますので、自分の目的に合ったものを取得しましょう。私はftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/archive/Hs37.3/CCDS.current.txtを取りました。

次にこのデータをRがパースしやすいようにPerlで整形します。元データはこんな感じのでデータになってます。

CCDS.current.txt:

#chromosome  nc_accession    gene    gene_id ccds_id ccds_status cds_strand  cds_from    cds_to  cds_locations   match_type
1   NC_000001.8 LINC00115   79854   CCDS1.1 Withdrawn   -   801942  802433  [801942-802433] Identical
1   NC_000001.10    SAMD11  148398  CCDS2.2 Public  +   861321  879532  [861321-861392, 865534-865715, 866418-866468, 871151-871275, 874419-874508, 874654-874839, 876523-876685, 877515-877630, 877789-877867, 877938-878437, 878632-878756, 879077-879187, 879287-879532] Identical
1   NC_000001.10    NOC2L   26155   CCDS3.1 Public  -   880073  894619  [880073-880179, 880436-880525, 880897-881032, 881552-881665, 881781-881924, 883510-883611, 883869-883982, 886506-886617, 887379-887518, 887791-887979, 888554-888667, 889161-889271, 889383-889461, 891302-891392, 891474-891594, 892273-892404, 892478-892652, 894308-894460, 894594-894619]   Identical
1   NC_000001.10    PLEKHN1 84069   CCDS4.1 Public  +   901911  909954  [901911-901993, 902083-902182, 905656-905802, 905900-905980, 906065-906137, 906258-906385, 906492-906587, 906703-906783, 907454-907529, 907667-907803, 908240-908389, 908565-908705, 908879-909019, 909212-909430, 909695-909743, 909821-909954]    Identical
1   NC_000001.10    HES4    57801   CCDS5.1 Public  -   934438  935352  [934438-934811, 934905-934992, 935071-935166, 935245-935352]    Identical
1   NC_000001.10    ISG15   9636    CCDS6.1 Public  +   948953  949857  [948953-948955, 949363-949857]  Identical
1   NC_000001.10    C1orf159    54991   CCDS7.2 Public  -   1018272 1026922 [1018272-1018366, 1019732-1019762, 1019860-1019885, 1021257-1021391, 1022518-1022583, 1022881-1022976, 1025732-1025807, 1026851-1026922]    Identical
1   NC_000001.10    TTLL10  254173  CCDS8.1 Public  +   1115433 1120521 [1115433-1115719, 1115862-1115980, 1116110-1116239, 1117120-1117194, 1117740-1117825, 1118255-1118426, 1119299-1119470, 1120348-1120521]    Identical
1   NC_000001.10    TNFRSF18    8784    CCDS9.1 Public  -   1138970 1141950 [1138970-1139339, 1139778-1139865, 1140749-1140871, 1141764-1141950]    Identical

CDS領域の範囲が一行に収められてます。 RマスターならRでもサクッといけるんでしょうが、私は以下のPerl スクリプトで事前にTSVファイルに変換しました。 (追記: CCDSの入力ファイルが0-based coordinateだったので1-basedに修正)

parse_ccds.pl:

#!/usr/bin/env perl

use v5.10;
use strict;
use warnings;

say "chromosome\tnc_accession\tgene\tgene_id\tccds_id\tccds_status\tcds_strand\tcds_from\tcds_to\tmatch_type\tstart\tstop";

while (<>) {
    next if /^#/;
    chomp;
    my ($chromosome, $nc_accession, $gene, $gene_id, $ccds_id, $ccds_status, $cds_strand, $cds_from, $cds_to, $cds_locations, $match_type) = split /\t/;
    my @cds_ranges = map { my ($start, $stop) = split /-/; {start => $start+1, stop => $stop+1} } (split /,\s/, (substr $cds_locations, 1, -1));
    foreach (@cds_ranges) {
        say "$chromosome\t$nc_accession\t$gene\t$gene_id\t$ccds_id\t$ccds_status\t$cds_strand\t$cds_from\t$cds_to\t$match_type\t$_->{start}\t$_->{stop}";
    }
}

そして、tsvファイルを吐きます。

cat CCDS.current.txt | perl parse_ccds.pl > CCDS.current.tsv

あとはRで読み込むだけです。幸い、makeGRangesFromDataFrameを使うことで、読み込んだデータフレームを簡単にGRangesオブジェクトに変換できます。

> library(GenomicRanges)
> ccds = read.csv("CCDS.current.tsv", sep="\t")
> ccds = makeGRangesFromDataFrame(ccds, strand.field="cds_strand", keep.extra.columns=T)
> ccds
GRanges object with 287116 ranges and 8 metadata columns:
           seqnames               ranges strand   | nc_accession      gene   gene_id     ccds_id ccds_status  cds_from    cds_to
              <Rle>            <IRanges>  <Rle>   |     <factor>  <factor> <integer>    <factor>    <factor> <integer> <integer>
       [1]        1     [801942, 802433]      -   |  NC_000001.8 LINC00115     79854     CCDS1.1   Withdrawn    801942    802433
       [2]        1     [861321, 861392]      +   | NC_000001.10    SAMD11    148398     CCDS2.2      Public    861321    879532
       [3]        1     [865534, 865715]      +   | NC_000001.10    SAMD11    148398     CCDS2.2      Public    861321    879532
       [4]        1     [866418, 866468]      +   | NC_000001.10    SAMD11    148398     CCDS2.2      Public    861321    879532
       [5]        1     [871151, 871275]      +   | NC_000001.10    SAMD11    148398     CCDS2.2      Public    861321    879532
       ...      ...                  ...    ... ...          ...       ...       ...         ...         ...       ...       ...
  [287112]        Y [16168463, 16168738]      +   |  NC_000024.9     VCY1B    353513 CCDS56618.1      Public  16168169  16168738
  [287113]        Y [16835028, 16835148]      +   |  NC_000024.9    NLGN4Y     22829 CCDS56619.1      Public  16835028  16953141
  [287114]        Y [16936067, 16936252]      +   |  NC_000024.9    NLGN4Y     22829 CCDS56619.1      Public  16835028  16953141
  [287115]        Y [16941609, 16942398]      +   |  NC_000024.9    NLGN4Y     22829 CCDS56619.1      Public  16835028  16953141
  [287116]        Y [16952292, 16953141]      +   |  NC_000024.9    NLGN4Y     22829 CCDS56619.1      Public  16835028  16953141
           match_type
             <factor>
       [1]  Identical
       [2]  Identical
       [3]  Identical
       [4]  Identical
       [5]  Identical
       ...        ...
  [287112]  Identical
  [287113]  Identical
  [287114]  Identical
  [287115]  Identical
  [287116]  Identical
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

遺伝子毎のCDS領域を得たければ、Gene IDのアノテーションがccdsにありますのでsplitでGRangesListオブジェクトにできます。

> ccdsbygene = split(ccds.public, ccds.public$gene_id)
> ccdsbygene
GRangesList object of length 18383:
$1 
GRanges object with 8 ranges and 8 metadata columns:
      seqnames               ranges strand | nc_accession     gene   gene_id     ccds_id ccds_status  cds_from    cds_to
         <Rle>            <IRanges>  <Rle> |     <factor> <factor> <integer>    <factor>    <factor> <integer> <integer>
  [1]       19 [58858387, 58858394]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [2]       19 [58858718, 58859005]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [3]       19 [58861735, 58862016]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [4]       19 [58862756, 58863052]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [5]       19 [58863648, 58863920]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [6]       19 [58864293, 58864562]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [7]       19 [58864657, 58864692]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
  [8]       19 [58864769, 58864802]      - |  NC_000019.9     A1BG         1 CCDS12976.1      Public  58858387  58864802
      match_type
        <factor>
  [1]  Identical
  [2]  Identical
  [3]  Identical
  [4]  Identical
  [5]  Identical
  [6]  Identical
  [7]  Identical
  [8]  Identical

...
<18382 more elements>
-------
seqinfo: 24 sequences from an unspecified genome; no seqlengths

りんごがでている

何か役に立つことを書きます

Consensus CDSの領域をGRangesオブジェクトにする