Rに読み込んだVCFに、アノテーションを付加しましょう。
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の機能を利用して行うこともできます。今後ご紹介していきます。
