CmhaDSO 範囲データの操作

ゲノム範囲

ゲノム参照配列を座標系とし、その範囲として各種特徴配列の位置を指定することができる. 此時の範囲をゲノム範囲といい、次の3つの情報からなる.

配列名

染色体名ともいう. scaffoldやcontigなどの配列を含めた広範な名称として配列名と称する.

範囲

範囲は特定の染色体/scaffold/contig配列上において単一の部分配列を表す整数区間である. 各範囲は開始位置と終了位置からなり、範囲の表現体系には次の2種類が存在する.

0ベースの座標系でかつ半開区間

0ベース半開区間は、コンピュータ内部においてメモリアドレスを指定するのに好適である。 この座標系では最初の塩基位置は0、範囲は[Start, end)と表現され、塩基数(範囲幅)は end - start で算出される。 制限酵素の切断位置など幅が0となる特徴配列を表現することもできる。

1ベースの座標系でかつ閉区間

1ベース閉区間では、最初の塩基位置は1、範囲は[Start, end]と表現され、塩基数(範囲幅)は end - start + 1 で算出される。

0ベース半開区間と1ベース閉区間の対比の明示・変換

0ベース半開区間はBED形式、対して1ベース開区間はPosition形式とも呼称され、それぞれchr1 100 120, chr1:101-120のように表現の仕方を変える場合がある。 前者を後者に変換する場合はstartに1を加算し、後者を前者に変換する場合はstartに1を減算する。 詳しくは ここ を参照するとよい。

ファイル形式など 範囲の表現体系
GTF/GFF 1ベース閉区間
SAM 1ベース閉区間
BAM 0ベース半開区間
VCF 1ベース閉区間
BCF 0ベース半開区間
BED 0ベース半開区間
Wiggle 1ベース閉区間
UCSC Genome Browser (WEB) 1ベース閉区間
UCSC Table Browser (Data) 0ベース半開区間
GenomicRanges 1ベース閉区間
BLAST 1ベース閉区間
GenBank/EMBL Feature Table 1ベース閉区間
R 文字列処理 1ベース閉区間
Python 文字列処理 0ベース半開区間
鎖

染色体DNAは二重鎖であるため、特徴配列は正鎖(順鎖、プラス鎖)、相補鎖(逆鎖、マイナス鎖)のいずれにも存在しうる. 染色体上の特徴配列の多くは鎖固有である. *は任意の鎖を表す.

BLAST結果などを例外としてほぼすべての範囲形式で、範囲の座標は参照配列の正鎖上で与えられる.

ゲノム範囲は参照ゲノムの特定のアセンブリバージョンと関連しており、どのバージョンに対するものかを指定する必要がある. 各種アセンブリバージョンの座標系間で様々なデータ形式の変換を行うツールにはCrossMapやLiftOver等がある.

GenomicRangesによる範囲データの操作

IRangesの基本操作

IRangesパッケージは整数区間の範囲情報(1ベース閉区間)を操作するための汎用手段を提供する. 範囲情報はIRangesクラスのオブジェクトに格納し操作する.

IRangesパッケージをロードする. 

> library(IRanges)

IRangesクラスのオブジェクトを作成する. 

> s <- c(4, 7, 2, 20) 
> e <- c(13, 7, 5, 23)  
> x <- IRanges(start = s, end = e)   # IRanges()コンストラクタ
> name(x) <- letters[1:4]
> x

各種情報にアクセスする. 

> class(x)                              # IRangesクラス
> start(x)                              # アクセサ関数
> end(x)                                # アクセサ関数
> width(x)                              # アクセサ関数

範囲操作を行う. 

> x +  nL         # 各範囲拡大
> x -  nL         # 各範囲縮小
> x *  nL         # 各範囲zoom-in
> x * -nL         # 各範囲zoom-out

> range(x)                # 範囲全体（両端）
> x[cdt]                  # subset 
> c(x, y)                 # concat

> restrict(x, s, e)           # [start, end]の範囲に制限

> flank(x, width=n)           # 幅nで上流領域
> flank(x, width=n, start=F)  # 幅nで下流領域

> reduce(x)                   # 重複のない範囲に集約
> gap(x)                      # 範囲間のギャップ
> gap(x, start=m, end=n)      # 両端を含めた範囲間のギャップ

> union(x, y)        # 和集合
> intersect(x, y)    # 積集合
> setdiff(x, y)      # 差集合

重複範囲

findOverlaps()関数を用いて、2つのIRangesオブジェクトの集合の重なりを見つける. queryHits列とsubjectHits列からなるHitsクラスのオブジェクトを返す. 重複はクエリと対象との間の写像を表す. 

> h <- findOverlaps(q, s) 
  # デフォルトはtype="any", 一部でも重なっていれば重複あり
  # デフォルトはselect="all", 全ての重複範囲を返す.
> h <- findOverlaps(q, s, type="within") # qの中でsに内包されているものが重複あり
> h <- findOverlaps(q, s, select="first")

固定されたsの範囲に対して多数のqから重複を見つけたい場合は区間木(NCListクラス)を利用することで計算コストを低減する. 

> s <- NCList(s)  
> h <- findOverlaps(q, s)  # 同様に処理可

アクセサ関数など 

> names(q)[queryHits(h)] 
> names(s)[subjectHits(h)]

> setNames(countQueryHits(h), names(q))     # countOverlaps(q, s)が同様の結果を返す.
> setNames(countSubjectHits(h), names(s))

> q[unique(queryHits(h))]                   # subsetByOverlaps(q, s)が同様の結果を返す. # 重複のあったqの範囲エントリのみ返す.

近傍範囲の検出と距離

近傍範囲エントリの検出 

> nearest(q, s)  # 最近傍(重複可)のsのエントリ番号(ベクトル)を返す  
> precede(q, s)  # qが先行(重複不可)する
> follow(q, s)   # qが後続(重複不可)する

距離の算出 

> distanceToNearest(q, s)  # Hitsオブジェクトにdistance列を追加したものを返す.
> distance(q, s)           # qとsの対になっているエントリ間の距離.

ランレングス符号化とview

カバレッジ(範囲重複の深度)データなどはしばしば同じ値の連続(run)を示す. runはランレングス符号化により圧縮できる. これはLengthsとValuesが組になったものである. Rle()関数はベクトルを入力とし、Rleオブジェクトを返す. 

> rle <- Rle(x)  # xはベクトル.
> runLength(rle)    # アクセサ関数.
> runValue(rle)     # アクセサ関数.
> as.vector(rle)    # 元に戻す.

coverage()関数はIRangesオブジェクトを入力とし、そのカバレッジをRleオブジェクトで返す. Rleオブジェクトに対してもサブセット化や統計量要約関数を適応できる. 

> rle <- coverage(x)  # xはIRanges.
> rle[rle > 3]        # subsetの例. 
> rle[y]                 # subsetの例. yはIRanges. yと重なる部分のみ返す.
> mean(rle[y])           # 要約関数適応例.

配列のサブセットを使用して新たな範囲を定義する. カバレッジピークの定義など. slice()関数はRleオブジェクトを入力とし、Viewsオブジェクトを返す. Viewsはランレングス符号化ベクトル(=配列ベクトル)と範囲を結合する. 

> v <- slice(rle, lower=n)  # n以上の値(カバレッジ等)を有するものを出力.
> ranges(v)                    # Viewsの範囲の部分のみ抽出.

Viewsは配列ベクトル(=Values)を特定の範囲(=Keys)ごとに集計要約することが容易である. 

> viewApply(v, mean)

ビンごとに配列を分割する場合. 

> views(rle, ir)   # rleをirで分割しViewsオブジェクトを作る.irはビンごとに分割された範囲.

ちなみにirは以下のように作るとよい. 

> bwidth <- nL                                # ビン幅を設定.
> end <- bwidth * floor(length(rle)/bwidth)   # ビン終端を設定.
> ir <- IRanges(start=seq(1, end, bwidth), width=bwidth)

GRangesの基本操作

GenomicRangesパッケージはゲノム範囲を操作するための手段を提供する. ゲノム範囲情報はGRangesクラスのオブジェクトに格納し操作する. GRangesはIRangesに配列名と鎖の情報を加えたものであり、任意でメタデータも追加できる.

GenomicRangesパッケージをロードする. 

> library(GenomicRanges)

GRangesオブジェクトはGRanges()コンストラクタにより作成できる. 追加の名前付き引数を指定することで任意のメタデータ列を追加できる. 

> gr <- GRanges(seqname = c("chr1", "chr1", "chr2", "chr3"),
                      ranges  = IRanges(start=5:8, width=10),
                      strand  = c("+", "-", "-", "+"), 
                      gc      = round(runif(4), 3)
              )
> seqlengths(gr) <- c(chr1=152, chr2=432, chr3=903)

すべてのメタデータはDataFrameに格納される. これはdata.frameと同様の振舞いをするが対応しているフィールド型が多い. カバレッジやギャップを計算するには配列長の情報を指定する必要がある. seqlengths()関数を用いて配列長を設定している.

アクセサ関数. subset. 

> start(gr)
> end(gr)
> width(gr)
> seqnames(gr)
> strand(gr)
> ranges(gr)

> gr[start(gr) > n]    # subset

> mcol(gr)                                 # メタデータへのアクセサ関数
> mcol(gr)$gc                              # ショートカット構文gr$gcも使える.
> mean(mcols(gr[seqname(gr)=="chr1"])$gc)  # chr1の全ての範囲のGC含量の平均値.

GRangesList

GRangesオブジェクトのリストをGRangesListといい、GRangesList()コンストラクタで作成でき、unlist()で再結合できる. 特定要素のアクセスはRのリストの場合と同様、[]はGRangesListを返し、[[]]はGRangesオブジェクトを返す. 実務上では主にsplit()関数の適用結果としてGRangesListが得られる. 

> grl < split(gr, seqname(gr))  # 配列名ごとにリスト分割. split-apply-combine patternを使用するため.
> sapply(grl, length)              # 配列名ごとにエントリ数を出力. 
> reduce(grl)                      # reduce(), flank(), coverage(), findOverlaps()などは自動的にリストの要素レベルで動作するよう設計されている. 
> unsplit(grl, seqname(gr))        # unsplit()関数は分割時と同じ因子で再結合する.

参考文献:

範囲データを操作するためのBioconductorパッケージ

概要

アノテーションデータ

GenomicFeatures

GenomicFeaturesはTxDbオブジェクトを作成・操作するための関数を提供する.

利用するためにはまずGenomicFeaturesおよび目的の転写アノテーションをインストール後、ロードする. 

> BiocManager::install("GenomicFeatures")
> BiocManager::install("TxDb.Mmusculus.UCSC.mm10.ensGene")
> library(TxDb.Mmusculus.UCSC.mm10.ensGene)
> txdb <- TxDb.Mmusculus.UCSC.mm10.ensGene                # ロード後、パッケージ名と同名のtxDbにアクセスできる.

目的のtxDbオブジェクトを含むアノテーションがない場合は作成する. 

> makeTxDbFromBiomart("ensembl", spicies_name)
> makeTxDbFromUCSC()
> makeTxDbFromGFF()
> saveDb()
> loadDb()

次にTxDbから情報抽出するための操作を示す. genes(), transcripts(), exons(), cds(), promoter()はTxDbから各情報を抽出し、GRangesオブジェクトとして返す. transcriptsBy(), exonsBy(), cdsBy(), intronBy()などはGRangesListオブジェクトとして返す.??? 特定染色体から抽出を行うにはseqlevels()を用いる. 指定範囲と重複するものを抽出するにはtranscriptsByOverlaps(), exonsByOverlaps(), cdsByOverlaps()などを用いる. 

> genes(txdb)                        # 全ての遺伝子情報

> exonsBy(txdb, by="tx")             # 転写産物ごとにリスト分割したエキソン情報

> seqlevels(txdb) <- "chr1"        # chr1に制限.    
> exonsBy(txdb, by="tx")             # 抽出.
> txdb <- restoreSeqlevels(txdb)   # 制限解除.

> transcriptByOverlaps(txdb, gr)     # grと重複する全ての転写産物情報.
BSgenome

BSgenomeパッケージは配列データを格納するためのクラスとそのアクセサ関数を含む.

利用するためにはまずBSgenomeおよび目的のゲノムパッケージをインストール後、ロードする. 

> BiocManager::install("BSgenome")
> BiocManager::install("BSgenome.Mmusculus.UCSC.mm10")
> library(BSgenome.Mmusculus.UCSC.mm10)
> mm10 <- BSgenome.Mmusculus.UCSC.mm10

アクセサ関数. 

> organism(mm10)            # Mus musculus
> providerVersion(mm10)     # mm10
> provider(mm10)            # UCSC
> seqinfo(mm10)             # 配列名や配列長などの情報
> mm10$chr1                 # chr1の配列
rtracklayer

ロードする. 

> library(rtracklayer)

import()関数は、様々な形式のファイルをその拡張子からデータ形式を検出し、GRangesオブジェクトとしてインポートする. import.bed(), import.gff(), import.wig()などの明示的に形式を指定した場合の関数もある. 

> mm_gtf <- import("in.gtf.gz")
> export(mm_gtf, con="out.bed", format="BED")

参考文献:

実践課題

プロモータ領域および配列の取得

BSgenomeオブジェクトbsgとゲノム範囲オブジェクトgrの染色体名の一致の確認と再マップ. 

> all(seqlevels(g) %in% seqlevels(gr))   # 一致確認. seqlevels()で配列名の確認・設定ができる.
> seqlevelsStyle(gr) <- "UCSC"         # seqlevelsStyle()により配列名をUCSCの命名体系に再マップ.

GTFファイルからタンパクコード遺伝子のプロモータ領域(上流n bp)を取得する. 

> gr_gtf <- import("in.gtf.gz")
> gr_pcg <- gr_gtf[gr_gtf$type=="gene" & gr_gtf$gtf_biotype=="protein_coding"] 
> flank(gr_pcg, width=n)                        # 既定で鎖を考慮する. ignore.strand=FALSE.
> promoters(gr_pcg, upstream=n, downstream=0)   # GenomicRangesのpromoters()関数で同様の処理.

ゲノム範囲オブジェクトgrとBSgenomeオブジェクトbsgを入力とし、grの各範囲の配列を取得する. 

> seqs <- getSeq(g, gr)

取得配列seqsをFASTAファイルにエクスポートする. 

> writeXStringSet(seqs, file="out.fasta", format="fasta")

遺伝子間およびイントロン領域の取得

gap()を用いたギャップ領域の抽出. Grangesオブジェクトgrに対し、gap(gr)とすると、全ての配列と全ての鎖の組合わせでギャップ範囲を出力する. 通常ギャップは鎖情報の如何を問わないため、*に置換する. 

> strand(gr) <- "*" 
> gaps(gr)[strand(gaps(gr))=="*"]

集合演算を用いたギャップ領域の抽出. 転写産物を含まない全ての遺伝子間領域を抽出する. 

> gr <- as(seqinfo(txdb), "GRanges")   # 染色体全体を表す範囲をGRamgesに強制変換
> collapsed_tx <- reduce(transcripts(txdb))
> strand(collapsed_tx) <- "*"
> setdiff(gr, collapsed_tx)

専用関数を用いたギャップ領域の抽出. 転写産物のイントロン領域を抽出する. 

> intronByTranscripts(txdb)

集合演算を用いたギャップ領域の抽出. 転写産物のイントロン領域を抽出する. エキソンと転写産物の両者を転写産物IDごとにリスト化し、setdiff()によりペアワイズ差集合を得る. 以下の例ではAmy1(ENSMUSG00000074264)のイントロン作成の場合に限定する. 

> amy1 <- transcriptsBy(txdb, "gene")$ENSMUSG00000074264 
> amy1_tx <- split(amy1, amy1$tx_id)  # Amy1転写産物の範囲をGRangesListに分割.リストの各要素名は転写産物ID.

> all_exons <- exonBy(txdb, "tx")  # GRangesListの各要素名は転写産物ID.

> amy1_exons <- all_exons[match(names(amy1_tx), names(all_exons))]  # amy1のエキソンのみ抽出.
> all(names(amy1_tx)==names(amy1_exons))   # すべて一致しているか確認.
> amy1_introns <- setdiff(amy1_tx, amy1_exons)   # ペアワイズ差集合.

GRangesにおける重複の検出・集計

染色体1の変異領域データgr_snpがあるものとする. 変異体長の長さ分布の要約統計量および最長のエントリのRS識別子を確認する. 

> summary(width(gr_snp))                    # 要約統計量
> names(gr_snp[which.max(width(gr_snp))])   # 最長エントリのRS識別子

次に変異がエキソン上に存在する比率を調べる. 幅が0の変異体(多くは参照ゲノムへの挿入に対応)は他の特徴配列と重複しないため、 重複演算を適用できるようにするために事前にサイズ変更をする必要がある. 

> collapsed_exons <- reduce(exons(txdb), ignore.strand=TRUE)
> chr1_collapsed_exons <- collapsed_exons[seqnames(collapsed_exons)=="chr1"]
> 
> gr_snp[width(gr_snp)==0] <- resize(gr_snp[width(gr_snp)==0], width=1) # 幅のリサイズ: 0 to 1
> 
> h <- findOverlaps(gr_snp, chr1_collapsed_exons)
> length(unique(queryHits(h))) / length(gr_snp)

エキソン上に存在する変異体の範囲データを出力する. 

> subsetByOverlaps(gr_snp, chr1_collapsed_exons, ignore.strand=TRUE)

各エキソン上に存在する変異体の数を数える. 

> chr1_collapsed_exons$N_snp <- countOverlaps(chr1_collapsed_exons, gr_snp, ignore.strand=TRUE)   # 引数順に注意

GRangesにおけるカバレッジ計算

以下の例ではマウスゲノム19番染色体上にランダムに150bpの偽リードを生成しカバレッジ計算をしている. 生成するリード数はLander-Watermanカバレッジ方程式(C=LN/G, C:カバレッジ、L:リード長、G:配列長、N:リード個数)から算出している. coverage(gr)は全ての染色体ごとに計算され、ランレングス符号化リストとして返される. 

> set.seed(0)
> G <- seqlengths(txdb)["chr19"]
> C <- 5 
> L <- 150 
> N <- floor(C*G/L)
> start <- sample(1:(G-L), N, replace=TRUE)
> gr <- GRanges("chr19", IRanges(start=start, width=L))
> 
> rle <- coverage(gr)
> 
> mean(rle)          # リードの平均カバレッジ
> sum(runLength(rle)[runValue(rle)==0])/G    # 19番染色体上でリードがカバーしていない領域の割合

Bedtoolsによる範囲データの操作

bedtools intersectを用いた重複抽出

bedtools intersect -a q.bed -b s.bed         # qの中から共通範囲を返す.
bedtools intersect -a q.bed -b s.bed -u      # 複数出力抑制.
bedtools intersect -a q.bed -b s.bed -v      # 重複しない範囲.
bedtools intersect -a q.bed -b s.bed -wo     # 重複塩基数も返す.
bedtools intersect -a q.bed -b s.bed -wa     # sに重なるqの範囲を返す.
bedtools intersect -a q.bed -b s.bed -wa -wb # 重複しているすべてのsとqの範囲を返す.
bedtools intersect -a q-sorted.bed -b s-sorted.bed --sorted # ソート済み入力データに対しメモリ消費の少ないアルゴリズム適用.

-sを用いた鎖情報.findOverlaps()との対比.

bedtools slop/bedtools flank

bedtools slopにより特徴範囲を拡張する. まずは終端を指定するための配列長ファイルg.txtを準備する. 

cat genome.fa | bioawk -c fastx '{print $name"\t"length($seq)}' > g.txt
 
bedtools slop -i q.bed -g g.txt -b m              # 両端をmbp拡張. 開始点や終了点を超過する場合はそれぞれ0又は配列長を返す.
bedtools slop -i q.bed -g g.txt -l m -r n  # 左端をmbp,右端をnbp拡張.

bedtools flankにより特徴範囲の近傍範囲を取得する. 

bioawk -cgff '{if($feature=="gene") print $0}' ref.gtf | grep 'gene biotype "protein_coding";' > pcg.gtf
bedtools flank -i pcg.gtf -g g.txt -l m -r n  # 各範囲の左mbp, 右nbpを抽出.

Promoter, Exon/Intron 

  awk 'BEGIN{FS=OFS="\t"} $3 == "transcript"' data.gtf & tx.gtf
  awk 'BEGIN{FS=OFS="\t"} $3 == "exon"' data.gtf       & exons.gtf 
  bedtools subtract -a tx.gtf -b exons.gtf             & introns.gtf
  bedtools getfasta -fi genome.fa -bed introns.gtf     & introns.fa

範囲データに対応する塩基配列を取得する

samtools faidx

samtools faidxにより索引を作ることでFASTAファイルから指定範囲の配列を効率的に取得することができる。 

samtools faidx xxx.fa                            # 索引xxx.fa.faiが生成される
samtools faidx xxx.fa chr:start:end ...   # 8:20000-25000 8:24000-32000のように範囲指定する。

範囲データの操作の際、終端座標を指定するために、各染色体長を知ることが必要となる場合が多くある。 この情報は参照ゲノムに対するFASTA索引ファイルの1,2列目を抽出することで得られる。 

samtools faidx genome.fa
cut -f 1,2 genome.fa.fai

bedtoolsを利用して塩基配列を取得することもできる。 

# 入力: 参照ゲノム塩基配列ファイル, BED(又はGTFがVCF) ファイル
# 出力: multi-fastaファイル
bedtools getfasta -fi genome.fa -bed region.bed