DNA配列をRで操作する
基本的なDNA配列の操作方法や、FASTA/FASTQ file を取り込む方法を解説します。また全ゲノム配列を読み込み操作する方法についても述べます。
準備
source("http://bioconductor.org/biocLite.R")
biocLite("Biostrings")
基本的なDNA, アミノ酸配列の操作
library("Biostrings")
x <- DNAString("actttGtag")
is(x)
## [1] "DNAString" "XString" "XRaw" "XVector" "Vector" "Annotated"
x
## 9-letter "DNAString" instance
## seq: ACTTTGTAG
XString クラスを継承した DNAString オブジェクトを作成しています。このオブジェケクトにしておくと、いろいろな配列操作が可能になります。RNA, アミノ酸は、それぞれ、RNAString and AAString を使います。
x[1:2]
## 2-letter "DNAString" instance
## seq: AC
aa <- translate(x)
is(aa)
## [1] "AAString" "XString" "XRaw" "XVector" "Vector" "Annotated"
aa
## 3-letter "AAString" instance
## seq: TL*
aa.revcomp <- translate(reverseComplement(x))
aa.revcomp
## 3-letter "AAString" instance
## seq: LQS
x.freq <- alphabetFrequency(x, baseOnly = TRUE)
is(x.freq)
## [1] "integer" "numeric" "vector"
## [4] "data.frameRowLabels" "EnumerationValue" "atomic"
## [7] "vectorORfactor"
x.freq
## A C G T other
## 2 1 2 4 0
DNA配列やコード表のデータ
もちろんDNAやコードのデータも用意されています。
DNA_BASES
## [1] "A" "C" "G" "T"
DNA_ALPHABET
## [1] "A" "C" "G" "T" "M" "R" "W" "S" "Y" "K" "V" "H" "D" "B" "N" "-" "+"
IUPAC_CODE_MAP
## A C G T M R W S Y K
## "A" "C" "G" "T" "AC" "AG" "AT" "CG" "CT" "GT"
## V H D B N
## "ACG" "ACT" "AGT" "CGT" "ACGT"
FASTA/FASTQ ファイルの読み込み
multi FASTA 形式のテキストファイルを読み込みます。デモデータは、
system.file("extdata", "someORF.fa", package="Biostrings")
にあります。
file <- system.file("extdata", "someORF.fa", package = "Biostrings")
seqs <- read.DNAStringSet(file, "fasta")
is(seqs)
## [1] "DNAStringSet" "XStringSet" "XRawList" "XVectorList"
## [5] "List" "Vector" "Annotated"
seqs
## A DNAStringSet instance of length 7
## width seq names
## [1] 5573 ACTTGTAAATATATCTTTTAT...ATCGACCTTATTGTTGATAT YAL001C TFC3 SGDI...
## [2] 5825 TTCCAAGGCCGATGAATTCGA...AAATTTTTTTCTATTCTCTT YAL002W VPS8 SGDI...
## [3] 2987 CTTCATGTCAGCCTGCACTTC...TACTCATGTAGCTGCCTCAT YAL003W EFB1 SGDI...
## [4] 3929 CACTCATATCGGGGGTCTTAC...CCCGAAACACGAAAAAGTAC YAL005C SSA1 SGDI...
## [5] 2648 AGAGAAAGAGTTTCACTTCTT...TAATTTATGTGTGAACATAG YAL007C ERP2 SGDI...
## [6] 2597 GTGTCCGGGCCTCGCAGGCGT...TTTTGGCAGAATGTACTTTT YAL008W FUN14 SGD...
## [7] 2780 CAAGATAATGTCAAAGTTAGT...AAGGAAGAAAAAAAAATCAC YAL009W SPO7 SGDI...
次に、デモ用の fastq file が system.file("extdata", "s_1_sequence.txt", package = "Biostrings") にありますので、これを取り込んでみます。
filepath <- system.file("extdata", "s_1_sequence.txt", package = "Biostrings")
fastq.geometry(filepath) # 統計情報を表示
## [1] 256 36
fastq <- read.DNAStringSet(filepath, format = "fastq")
全ゲノムDNA配列の操作
全ゲノム配列をネットからダウンロードします。時間がかかるので、日本ミラーを設定することをお勧めします。
参考: bioconductor のパッケージを高速にインストールする (bioconductor.jp の設定)
ここではマウスゲノム mm10 をダウンロードします。
source("http://bioconductor.org/biocLite.R")
biocLite("BSgenome.Mmusculus.UCSC.mm10")
ゲノムデータを読み込みます。
library("BSgenome.Mmusculus.UCSC.mm10")
Mmusculus
## Mouse genome
## |
## | organism: Mus musculus (Mouse)
## | provider: UCSC
## | provider version: mm10
## | release date: Dec. 2011
## | release name: Genome Reference Consortium GRCm38
## |
## | sequences (see '?seqnames'):
## | chr1 chr2 chr3
## | chr4 chr5 chr6
## | chr7 chr8 chr9
## | chr10 chr11 chr12
## | chr13 chr14 chr15
## | chr16 chr17 chr18
## | chr19 chrX chrY
## | chrM chr1_GL456210_random chr1_GL456211_random
## | chr1_GL456212_random chr1_GL456213_random chr1_GL456221_random
## | chr4_GL456216_random chr4_GL456350_random chr4_JH584292_random
## | chr4_JH584293_random chr4_JH584294_random chr4_JH584295_random
## | chr5_GL456354_random chr5_JH584296_random chr5_JH584297_random
## | chr5_JH584298_random chr5_JH584299_random chr7_GL456219_random
## | chrX_GL456233_random chrY_JH584300_random chrY_JH584301_random
## | chrY_JH584302_random chrY_JH584303_random chrUn_GL456239
## | chrUn_GL456359 chrUn_GL456360 chrUn_GL456366
## | chrUn_GL456367 chrUn_GL456368 chrUn_GL456370
## | chrUn_GL456372 chrUn_GL456378 chrUn_GL456379
## | chrUn_GL456381 chrUn_GL456382 chrUn_GL456383
## | chrUn_GL456385 chrUn_GL456387 chrUn_GL456389
## | chrUn_GL456390 chrUn_GL456392 chrUn_GL456393
## | chrUn_GL456394 chrUn_GL456396 chrUn_JH584304
## |
## | (use the '$' or '[[' operator to access a given sequence)
is(Mmusculus)
## [1] "BSgenome" "GenomeDescription"
seqnames(Mmusculus)
## [1] "chr1" "chr2" "chr3"
## [4] "chr4" "chr5" "chr6"
## [7] "chr7" "chr8" "chr9"
## [10] "chr10" "chr11" "chr12"
## [13] "chr13" "chr14" "chr15"
## [16] "chr16" "chr17" "chr18"
## [19] "chr19" "chrX" "chrY"
## [22] "chrM" "chr1_GL456210_random" "chr1_GL456211_random"
## [25] "chr1_GL456212_random" "chr1_GL456213_random" "chr1_GL456221_random"
## [28] "chr4_GL456216_random" "chr4_GL456350_random" "chr4_JH584292_random"
## [31] "chr4_JH584293_random" "chr4_JH584294_random" "chr4_JH584295_random"
## [34] "chr5_GL456354_random" "chr5_JH584296_random" "chr5_JH584297_random"
## [37] "chr5_JH584298_random" "chr5_JH584299_random" "chr7_GL456219_random"
## [40] "chrX_GL456233_random" "chrY_JH584300_random" "chrY_JH584301_random"
## [43] "chrY_JH584302_random" "chrY_JH584303_random" "chrUn_GL456239"
## [46] "chrUn_GL456359" "chrUn_GL456360" "chrUn_GL456366"
## [49] "chrUn_GL456367" "chrUn_GL456368" "chrUn_GL456370"
## [52] "chrUn_GL456372" "chrUn_GL456378" "chrUn_GL456379"
## [55] "chrUn_GL456381" "chrUn_GL456382" "chrUn_GL456383"
## [58] "chrUn_GL456385" "chrUn_GL456387" "chrUn_GL456389"
## [61] "chrUn_GL456390" "chrUn_GL456392" "chrUn_GL456393"
## [64] "chrUn_GL456394" "chrUn_GL456396" "chrUn_JH584304"
## chr19 の情報を取り出す
mm10.chr19 <- Mmusculus$chr19
is(mm10.chr19)
## [1] "MaskedDNAString" "MaskedXString"
4種類のゲノム配列マスクデータが含まれている。AGAPS は、いわゆるUCSCゲノムブラウザの Gap Track の部分が N でマスクされていることを示す。コンティグの間をNで埋めている。つまりアセンブリギャップ。AMB は、コンティグ内で読むことができなかった塩基を IUPAC ambiguity letter でマスクしている。RM は RepeatMasker でマスクしている。TRF は tandem repeat をマスクしたもの。
データを読み込んだときのデフォルトは、AGAPSとAMBだけが有効になっているゲノム配列が得られる。例えば、RM を有効にするには、
active(masks(mm10.chr19))["RM"] <- TRUE
## Warning: 置き換えるべき項目数が，置き換える数の倍数ではありませんでした
countPattern("AAGAACAT", mm10.chr19) # pattern をカウント
## [1] 2352
パターンマッチング
mm10.chr19.match <- matchPattern("AAGAACAT", mm10.chr19, max.mismatch = 1)
is(mm10.chr19.match)
## [1] "XStringViews" "Views" "List" "Vector"
## [5] "Annotated"
mm10.chr19.match
## Views on a 61431566-letter DNAString subject
## subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
## views:
## start end width
## [1] 3000189 3000196 8 [AAGAAAAT]
## [2] 3005005 3005012 8 [AAGAACCT]
## [3] 3005091 3005098 8 [AAGAACCT]
## [4] 3006251 3006258 8 [AAGATCAT]
## [5] 3009789 3009796 8 [AACAACAT]
## [6] 3009818 3009825 8 [AAGAATAT]
## [7] 3010972 3010979 8 [AAGAACAG]
## [8] 3011103 3011110 8 [CAGAACAT]
## [9] 3011138 3011145 8 [AAGAAGAT]
## ... ... ... ... ...
## [49007] 61323881 61323888 8 [AAGAAAAT]
## [49008] 61324346 61324353 8 [AAAAACAT]
## [49009] 61324511 61324518 8 [AAGAATAT]
## [49010] 61324707 61324714 8 [AAGAACAT]
## [49011] 61325474 61325481 8 [AAGATCAT]
## [49012] 61325642 61325649 8 [AAGAGCAT]
## [49013] 61327589 61327596 8 [AACAACAT]
## [49014] 61329266 61329273 8 [ATGAACAT]
## [49015] 61330773 61330780 8 [AAGGACAT]
配列を取り出す
mm10.chr19.dna <- as(mm10.chr19, "XStringViews")
is(mm10.chr19.dna)
## [1] "XStringViews" "Views" "List" "Vector"
## [5] "Annotated"
mm10.chr19.dna
## Views on a 61431566-letter DNAString subject
## subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
## views:
## start end width
## [1] 3000001 6685768 3685768 [GATCATACTCCTCATGCTGG...GCAAGCTTCCGGGTCAGAGA]
## [2] 6685869 6688105 2237 [GCACACACACATACACACAC...ACCTGTCTGTGTGGTACTGG]
## [3] 6813616 6829746 16131 [GCCTGCTCTCCTAGGAGTAC...GAGTGGGAGCCGCTGGCAAA]
## [4] 6829847 61331566 54501720 [GATCAGAGTGGGAGCCGCTG...AGGGTTAGGGTTAGGGTTAG]
start(mm10.chr19.dna[1]) # コンティグのスタートの座標
## [1] 3000001
mm10.chr19.dna[[1]][1:100] # 3000001 bp から 100 bp 取り出す
## 100-letter "DNAString" instance
## seq: GATCATACTCCTCATGCTGGACATTCTGGTTCCT...TCAACTTTCCCTCTTCCTATGCCAGTTATGCAT