アメリエフの技術ブログ

Amelieff Staff Blog

(Rで)VCFにアノテーション "VariantAnnotation" ②

Rに読み込んだVCFに、アノテーションを付加しましょう。

staffblog.amelieff.jp

VariantAnnotationの公式説明書はこちらです。

アノテーションの付加には、Bioconducterのアノテーションデータパッケージを利用します。
パッケージ一覧はこちらから見られます。 Bioconductor - 3.7 AnnotationData Packages

アノテーションパッケージの読み込み

例として、変異が遺伝子のどの領域に位置するか(例: coding, intron)の情報を付加します。
アノテーションパッケージ TxDb.Hsapiens.UCSC.hg38.knownGene を読み込みます。
このとき、VCFとパッケージのゲノムビルドを揃えることに注意してください。

> vcf <- readVcf("ERX100000.vcf", "hg38")  # VCFの読み込み

> biocLite("TxDb.Hsapiens.UCSC.hg38.knownGene")  #アノテーションパッケージ
> library(TxDb.Hsapiens.UCSC.hg38.knownGene)

アノテーション付加

公式の説明書では、読み込んだVCFの rowRanges(vcf) (paramRangeID, REF, ALT, QUAL, FILTER)を抽出し、それにアノテーションすることが推奨されています。
変異の位置情報を付与する場合は、 locateVariants() を使用します。

> rd <- rowRanges(vcf)   # rowRangesを抽出

# locateVariants(アノテーションされるデータ, アノテーションパッケージ, 付与する情報)
> loc_code <- locateVariants(rd, TxDb.Hsapiens.UCSC.hg38.knownGene, CodingVariants())  # CordingVariantの変異に付加
> loc_all <- locateVariants(rd, TxDb.Hsapiens.UCSC.hg38.knownGene, AllVariants())  # 全ての変異に付加

アノテーション結果の閲覧

LOCATION(変異の位置)列以降が付加されています。
例えば1行目の変異はintergenic 、遺伝子間であることがわかります。

> class(loc_all)  
[1] "GRanges"
attr(,"package")
[1] "GenomicRanges"

> head(loc_all,4)
GRanges object with 4 ranges and 9 metadata columns:
                 seqnames         ranges strand |   LOCATION  LOCSTART
                    <Rle>      <IRanges>  <Rle> |   <factor> <integer>
  chr1:14464_A/T     chr1 [14464, 14464]      * | intergenic      <NA>
  chr1:14522_G/A     chr1 [14522, 14522]      * | intergenic      <NA>
  chr1:14542_A/G     chr1 [14542, 14542]      * | intergenic      <NA>
  chr1:17005_A/G     chr1 [17005, 17005]      * | intergenic      <NA>
                    LOCEND   QUERYID        TXID         CDSID      GENEID
                 <integer> <integer> <character> <IntegerList> <character>
  chr1:14464_A/T      <NA>         1        <NA>                      <NA>
  chr1:14522_G/A      <NA>         2        <NA>                      <NA>
  chr1:14542_A/G      <NA>         3        <NA>                      <NA>
  chr1:17005_A/G      <NA>         4        <NA>                      <NA>
                                         PRECEDEID        FOLLOWID
                                   <CharacterList> <CharacterList>
  chr1:14464_A/T 100130417,100132287,100302278,...                
  chr1:14522_G/A 100130417,100132287,100302278,...                
  chr1:14542_A/G 100130417,100132287,100302278,...                
  chr1:17005_A/G 100130417,100132287,100302278,...                
  -------
  seqinfo: 25 sequences from hg38 genome

LOCATIONにはどのような分類があるか見てみます。

> levels(loc_all$LOCATION)
[1] "spliceSite" "intron"     "fiveUTR"    "threeUTR"   "coding"    
[6] "intergenic" "promoter"  

> summary(loc_all$LOCATION)
spliceSite     intron    fiveUTR   threeUTR     coding intergenic 
       187     168767       6091      34545      46521      17017 
  promoter 
     39771 

これまでで、データの列名や位置情報の分類名がわかりましたので、
条件を絞り込んで見ることもできます。

> head(subset(loc_all, LOCATION == "coding"))    #LOCATION(位置)がcordingである変異を表示
> head(subset(loc_all, !is.na(GENEID)))   # GENEIDがある(<NA>以外)変異を表示

変異の絞り込みは、VariantAnnotationの機能を利用して行うこともできます。今後ご紹介していきます。