fasta配列操作
S.Konishi
2017-07-24
1 fasta ファイル操作
1.1 fastaファイル読み込み
- GRCm38.cdna(cDNA), GRCm38.pep(アミノ酸)
readDNAStringSet,readAAStringSet
# ライブラリロード ----
suppressWarnings(suppressMessages(library(Biostrings)))
# fasta読み込み gzファイル読み込み可 ----
fn <- "~/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz"
fa <- "~/db/protein/Mus_musculus.GRCm38.pep.all.fa.gz"
fna <- readDNAStringSet(fn) # cdna
faa <- readAAStringSet(fa) # amino acid
fna;faa## A DNAStringSet instance of length 107368
## width seq names
## [1] 14 GGGACTGGGGGGGC ENSMUST0000017886...
## [2] 12 GGGACAGGGGGC ENSMUST0000017853...
## [3] 9 ATGGCATAT ENSMUST0000019622...
## [4] 11 ATGGCATATCA ENSMUST0000017966...
## [5] 16 ATCGGAGGGATACGAG ENSMUST0000017756...
## ... ... ...
## [107364] 556 ATGAAGCAGCAGCTGCTG...TCCAGGCCCTAGAATGA ENSMUST0000020080...
## [107365] 664 ATGAGAAGAATGGCTCTT...ACTTGGGGACCAAGAAC ENSMUST0000019586...
## [107366] 545 ATGAAGCTGCTGCTGCTG...TCCAGGCCCTAGAATGA ENSMUST0000020150...
## [107367] 549 ATGCAGCTGCTGCTGCTG...TCCAGGCCCTAGAATGA ENSMUST0000020097...
## [107368] 545 ATGAAGCTGCTGCTGCTG...TCCAGGCCCTAGAATGA ENSMUST0000020268...## A AAStringSet instance of length 57751
## width seq names
## [1] 4 GTGG ENSMUSP0000014176...
## [2] 5 GTGGA ENSMUSP0000014131...
## [3] 5 IGGIR ENSMUSP0000014295...
## [4] 3 MAY ENSMUSP0000014254...
## [5] 4 LTGT ENSMUSP0000014215...
## ... ... ...
## [57747] 553 MAFPELLDRVGGLGRFQLF...KKVTHDTPDGSILMSTRL ENSMUSP0000014452...
## [57748] 274 XRFQLFQTVALVTPILWVT...IFLLQALIGIVDFPVKTG ENSMUSP0000014437...
## [57749] 223 MACLGLRRYKAQLQLPSRT...PTEPECEKQFQPYFIPIN ENSMUSP0000014153...
## [57750] 102 MACLGLRRYKAQLQLPSRT...CATTFTEKNTVGFLDLCH ENSMUSP0000014398...
## [57751] 200 MKPYFCHVFVFCFLIRLLT...FMAAVNTNKKSRLAGVTS ENSMUSP0000014122...1.2 配列数、配列長、fasta index
- 配列数
length, 配列長width - ファイル読み込み無しで配列長を調べる
fasta.seqlength(filepath) - fasta indexを作成
fasta.index
length(fna); width(fna)[1:50]## [1] 107368## [1] 14 12 9 11 16 11 16 10 12 17 10 17 29 10 17 29 10
## [18] 17 29 10 17 17 23 17 23 416 334 417 560 455 368 393 389 366
## [35] 334 385 340 340 368 370 398 374 566 335 381 290 347 344 367 344head(fasta.seqlengths(fn))## ENSMUST00000178862.1 cdna chromosome:GRCm38:6:41542163:41542176:1 gene:ENSMUSG00000094569.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trbd2 description:T cell receptor beta, D region 2 [Source:MGI Symbol;Acc:MGI:4439727]
## 14
## ENSMUST00000178537.1 cdna chromosome:GRCm38:6:41533201:41533212:1 gene:ENSMUSG00000095668.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trbd1 description:T cell receptor beta, D region 1 [Source:MGI Symbol;Acc:MGI:4439571]
## 12
## ENSMUST00000196221.1 cdna chromosome:GRCm38:14:54113468:54113476:1 gene:ENSMUSG00000096749.2 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trdd1 description:T cell receptor delta diversity 1 [Source:MGI Symbol;Acc:MGI:4439547]
## 9
## ENSMUST00000179664.1 cdna chromosome:GRCm38:14:54113468:54113478:1 gene:ENSMUSG00000096749.2 gene_biotype:TR_D_gene transcript_biotype:processed_transcript gene_symbol:Trdd1 description:T cell receptor delta diversity 1 [Source:MGI Symbol;Acc:MGI:4439547]
## 11
## ENSMUST00000177564.1 cdna chromosome:GRCm38:14:54122226:54122241:1 gene:ENSMUSG00000096176.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trdd2 description:T cell receptor delta diversity 2 [Source:MGI Symbol;Acc:MGI:4439546]
## 16
## ENSMUST00000179520.1 cdna chromosome:GRCm38:12:113430528:113430538:-1 gene:ENSMUSG00000094028.1 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene gene_symbol:Ighd4-1 description:immunoglobulin heavy diversity 4-1 [Source:MGI Symbol;Acc:MGI:4439801]
## 11# fasta.index()
faidx <- fasta.index(fn)
knitr::kable(head(faidx), format = "markdown", align = "r", caption="fasta index")| recno | fileno | offset | desc | seqlength | filepath |
|---|---|---|---|---|---|
| 1 | 1 | 0 | ENSMUST00000178862.1 cdna chromosome:GRCm38:6:41542163:41542176:1 gene:ENSMUSG00000094569.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trbd2 description:T cell receptor beta, D region 2 [Source:MGI Symbol;Acc:MGI:4439727] | 14 | /Users/shkonishi/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz |
| 2 | 1 | 259 | ENSMUST00000178537.1 cdna chromosome:GRCm38:6:41533201:41533212:1 gene:ENSMUSG00000095668.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trbd1 description:T cell receptor beta, D region 1 [Source:MGI Symbol;Acc:MGI:4439571] | 12 | /Users/shkonishi/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz |
| 3 | 1 | 516 | ENSMUST00000196221.1 cdna chromosome:GRCm38:14:54113468:54113476:1 gene:ENSMUSG00000096749.2 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trdd1 description:T cell receptor delta diversity 1 [Source:MGI Symbol;Acc:MGI:4439547] | 9 | /Users/shkonishi/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz |
| 4 | 1 | 772 | ENSMUST00000179664.1 cdna chromosome:GRCm38:14:54113468:54113478:1 gene:ENSMUSG00000096749.2 gene_biotype:TR_D_gene transcript_biotype:processed_transcript gene_symbol:Trdd1 description:T cell receptor delta diversity 1 [Source:MGI Symbol;Acc:MGI:4439547] | 11 | /Users/shkonishi/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz |
| 5 | 1 | 1041 | ENSMUST00000177564.1 cdna chromosome:GRCm38:14:54122226:54122241:1 gene:ENSMUSG00000096176.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trdd2 description:T cell receptor delta diversity 2 [Source:MGI Symbol;Acc:MGI:4439546] | 16 | /Users/shkonishi/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz |
| 6 | 1 | 1304 | ENSMUST00000179520.1 cdna chromosome:GRCm38:12:113430528:113430538:-1 gene:ENSMUSG00000094028.1 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene gene_symbol:Ighd4-1 description:immunoglobulin heavy diversity 4-1 [Source:MGI Symbol;Acc:MGI:4439801] | 11 | /Users/shkonishi/db/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz |
1.3 配列操作
toString配列を取り出すcomplement相補鎖,reverseComplement逆相補鎖,reverse逆鎖
toString(fna[[1]]) # 配列を取り出し as.character(fna[[1]])でもO.K.## [1] "GGGACTGGGGGGGC"toString(complement(fna[[1]])) # 相補鎖## [1] "CCCTGACCCCCCCG"toString(reverseComplement(fna[[1]])) # 逆相補鎖revcom## [1] "GCCCCCCCAGTCCC"toString(reverse(fna[[1]])) # 逆鎖## [1] "CGGGGGGGTCAGGG"1.4 塩基, アミノ酸ごとの出現頻度
alphabetFrequency,dinucleotideFrequency,trinucleotideFrequencyoligonucleotideFrequencyk-merカウント
alphabetFrequency(fna[[1]]) # 塩基## A C G T M R W S Y K V H D B N - + .
## 1 2 10 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0alphabetFrequency(faa[[1]]) # アミノ酸## A R N D C Q E G H I L K
## 0 0 0 0 0 0 0 3 0 0 0 0
## M F P S T W Y V U O B J
## 0 0 0 0 1 0 0 0 0 0 0 0
## Z X * - + . other
## 0 0 0 0 0 0 0dinucleotideFrequency(fna[[1]]) ## AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
## 0 1 0 0 0 0 0 1 1 1 8 0 0 0 1 0trinucleotideFrequency(fna[[1]])## AAA AAC AAG AAT ACA ACC ACG ACT AGA AGC AGG AGT ATA ATC ATG ATT CAA CAC
## 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## CAG CAT CCA CCC CCG CCT CGA CGC CGG CGT CTA CTC CTG CTT GAA GAC GAG GAT
## 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0
## GCA GCC GCG GCT GGA GGC GGG GGT GTA GTC GTG GTT TAA TAC TAG TAT TCA TCC
## 0 0 0 0 1 1 6 0 0 0 0 0 0 0 0 0 0 0
## TCG TCT TGA TGC TGG TGT TTA TTC TTG TTT
## 0 0 0 0 1 0 0 0 0 0oligonucleotideFrequency(fna[[1]], width=5)[1:30]## AAAAA AAAAC AAAAG AAAAT AAACA AAACC AAACG AAACT AAAGA AAAGC AAAGG AAAGT
## 0 0 0 0 0 0 0 0 0 0 0 0
## AAATA AAATC AAATG AAATT AACAA AACAC AACAG AACAT AACCA AACCC AACCG AACCT
## 0 0 0 0 0 0 0 0 0 0 0 0
## AACGA AACGC AACGG AACGT AACTA AACTC
## 0 0 0 0 0 01.5 ヘッダ部分取り出し編集
# ヘッダ部分
names(fna)[1:3] ## [1] "ENSMUST00000178862.1 cdna chromosome:GRCm38:6:41542163:41542176:1 gene:ENSMUSG00000094569.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trbd2 description:T cell receptor beta, D region 2 [Source:MGI Symbol;Acc:MGI:4439727]"
## [2] "ENSMUST00000178537.1 cdna chromosome:GRCm38:6:41533201:41533212:1 gene:ENSMUSG00000095668.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trbd1 description:T cell receptor beta, D region 1 [Source:MGI Symbol;Acc:MGI:4439571]"
## [3] "ENSMUST00000196221.1 cdna chromosome:GRCm38:14:54113468:54113476:1 gene:ENSMUSG00000096749.2 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:Trdd1 description:T cell receptor delta diversity 1 [Source:MGI Symbol;Acc:MGI:4439547]"# ENSEMBLE protein idを取り出し
ensp <- sapply(strsplit(names(fna), " "), '[', 1)
ensp[1:3]## [1] "ENSMUST00000178862.1" "ENSMUST00000178537.1" "ENSMUST00000196221.1"# gene_symbolを取り出し
sym <- sub("gene_symbol:","",
sapply(strsplit(names(fna), " "), function(x){x[grep("gene_symbol:",x)]}))
sym[1:3]## [1] "Trbd2" "Trbd1" "Trdd1"1.6 fastaファイル書き出し
writeXStringSet
# キーワード検索してfastaオブジェクト抽出
str <- "Bnip2" # GPCR
fstobj <- fna[grep(str, sym, ignore.case = T)]
## 書き出し
# writeXStringSet(fstobj, file="select.faa", format="fasta", width=50)2 マルチプルアラインメント
2.1 Example: Protein coding sequences
AlignSeqs,AlignTranslationBrowseSeqs
suppressMessages(suppressWarnings(library(DECIPHER)))
# example sequence 317seqs, 平均830bp
fas <- system.file("extdata", "50S_ribosomal_protein_L2.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
dna ## A DNAStringSet instance of length 317
## width seq names
## [1] 819 ATGGCTTTAAAAAATTTTAA...TATTGTAAAAAAAAGAAAA Rickettsia prowaz...
## [2] 822 ATGGGAATACGTAAACTCAA...CATTGAGAGAAGGAAAAAG Porphyromonas gin...
## [3] 822 ATGGGAATACGTAAACTCAA...CATTGAGAGAAGGAAAAAG Porphyromonas gin...
## [4] 822 ATGGGAATACGTAAACTCAA...CATTGAGAGAAGGAAAAAG Porphyromonas gin...
## [5] 819 ATGGCTATCGTTAAATGTAA...CGTACGTCGTCGTGGTAAA Pasteurella multo...
## ... ... ...
## [313] 819 ATGGCAATTGTTAAATGTAA...CGTACGTCGCCGTACTAAA Pectobacterium at...
## [314] 822 ATGCCTATTCAAAAATGCAA...TCGCGATCGTCGCGTCAAG Acinetobacter sp....
## [315] 864 ATGGGCATTCGCGTTTACCG...TCGCGGTGGTCGTCAGTCT Thermosynechococc...
## [316] 831 ATGGCACTGAAGACATTCAA...CCGCCACAAGCGGAAGAAG Bradyrhizobium ja...
## [317] 840 ATGGGCATTCGCAAATATCG...GACGGCTTCCGGGCGAGGT Gloeobacter viola...# マルチプルアラインメント
DNA <- AlignSeqs(dna, verbose = F) # align the sequences directly without translation
# AAに翻訳してマルチプルアラインメント
AA <- AlignTranslation(dna, type="AAStringSet", verbose=F) # align the translation
# htmlで出力 highlight NA (all seq), 0(differ seq from consensus)
BrowseSeqs(AA, highlight=0)3 fastaファイル(ゲノムシーケンス)操作
3.1 ファイル読み込み
- データ
GCA_000448345.1(Cgr1.0 17A/GY),GCA_000419365.1(C_griseus1.0),GCA_000223135.1( CriGri_1.0 )
# a. fastaファイル読み込み ----
f_cgr <- "~/db/genome/Cgr1.0/GCA_000448345.1_Cgr1.0_genomic.fna"
f_c_gr <- "~/db/genome/C_griseus1.0/GCA_000419365.1_C_griseus_v1.0_genomic.fna"
f_crigri <- "~/db/genome/CriGri_1.0/CriGri_1.0.fa"
system.time( f.cgr <- readDNAStringSet(f_cgr) ) ## user system elapsed
## 9.886 3.243 14.203system.time( f.c_gr <- readDNAStringSet(f_c_gr) ) ## user system elapsed
## 9.217 3.791 14.783system.time( f.crigri <- readDNAStringSet(f_crigri)) ## user system elapsed
## 10.765 4.389 21.0103.2 塩基ごとの出現頻度, Nの割合, GCの割合
# 塩基ごとの出現頻度 ----
b.cgr <- apply(alphabetFrequency(f.cgr), 2, sum) # 名前付きベクトル
b.c_gr <- apply(alphabetFrequency(f.c_gr), 2, sum)
b.crigri <- apply(alphabetFrequency(f.crigri), 2, sum)
nuc <- data.frame(cbind(Cgr1.0=b.cgr, C_griseus1.0=b.c_gr, CriGri_1.0=b.crigri))
knitr::kable(nuc, format = "markdown", align = "r", caption="C. griceus genome contents")| Cgr1.0 | C_griseus1.0 | CriGri_1.0 | |
|---|---|---|---|
| A | 612896951 | 674433945 | 679960020 |
| C | 430937235 | 476451287 | 479592549 |
| G | 431194278 | 476151522 | 479353877 |
| T | 613965601 | 674272879 | 679225796 |
| M | 11240 | 0 | 0 |
| R | 23120 | 0 | 0 |
| W | 10003 | 0 | 0 |
| S | 6824 | 0 | 0 |
| Y | 22657 | 0 | 0 |
| K | 11009 | 0 | 0 |
| V | 422 | 0 | 0 |
| H | 517 | 0 | 0 |
| D | 493 | 0 | 0 |
| B | 331 | 0 | 0 |
| N | 243693609 | 58820511 | 81654506 |
| - | 0 | 0 | 0 |
| + | 0 | 0 | 0 |
| . | 0 | 0 | 0 |
# Nの割合----
(nrate <- c(
sum(b.cgr[names(b.cgr) %in% c("N")])/sum(as.numeric(b.cgr)),
sum(b.c_gr[names(b.c_gr) %in% c("N")])/sum(as.numeric(b.c_gr)),
sum(b.crigri[names(b.crigri) %in% c("N")])/sum(as.numeric(b.crigri))
)
)## [1] 0.10446515 0.02492257 0.03402573# GC ----
(gcrate <- c(
sum(b.cgr[names(b.cgr) %in% c("C","G")])/
sum(b.cgr[names(b.cgr) %in% c("A","T","C","G")]),
sum(b.c_gr[names(b.c_gr) %in% c("C","G")])/
sum(b.cgr[names(b.c_gr) %in%c("A","T","C","G")]),
sum(b.crigri[names(b.crigri) %in% c("C","G")])/
sum(b.cgr[names(b.crigri) %in%c("A","T","C","G")])
)
)## [1] 0.4127018 0.4560103 0.45904703.3 N50
# N50 ----
MYNX <- function(fstobj, N){
sortlen <- sort(as.numeric(width(fstobj)), decreasing = T) # 配列を長い順に並べ
for(i in 1:length(sortlen)){
sumlen <- sum(sortlen[1:i]) # 配列を長い順に並べて上から順に足していく
if(sumlen >= sum(as.numeric(width(fstobj)))/(100/N)){
nx <- i
break
}
}
return(sortlen[nx])
}
MYNX(fstobj = f.cgr, N = 50); MYNX(fstobj = f.c_gr, N = 50); MYNX(fstobj = f.crigri, N = 50)## [1] 1236516## [1] 1558295## [1] 11472333.4 まとめ
- 全塩基数はas.numericにしないと整数オーバーフロー
dat <-
data.frame(
name =c("Cgr1.0", "C_griseus1.0","CriGri_1.0"),
assembly=c("GCA_000448345.1", "GCA_000419365.1", "GCA_000223135.1"),
# 全塩基数,
bp=c(sum(as.numeric(width(f.cgr))),
sum(as.numeric(width(f.c_gr))),
sum(as.numeric(width(f.crigri)))),
# 配列数,
scaffold=c(length(f.cgr),
length(f.c_gr),
length(f.crigri)),
# N80,N50,N20
N80=c(MYNX(fstobj = f.cgr, N = 80),
MYNX(fstobj = f.c_gr, N = 80),
MYNX(fstobj = f.crigri, N = 80)),
N50=c(MYNX(fstobj = f.cgr, N = 50),
MYNX(fstobj = f.c_gr, N = 50),
MYNX(fstobj = f.crigri, N = 50)),
N20=c(MYNX(fstobj = f.cgr, N = 20),
MYNX(fstobj = f.c_gr, N = 20),
MYNX(fstobj = f.crigri, N = 20)),
N = nrate,
GC =gcrate
)
knitr::kable(dat, format = "markdown", align = "r", caption="C.griceus genome")| name | assembly | bp | scaffold | N80 | N50 | N20 | N | GC |
|---|---|---|---|---|---|---|---|---|
| Cgr1.0 | GCA_000448345.1 | 2332774290 | 28749 | 430945 | 1236516 | 2787033 | 0.1044651 | 0.4127018 |
| C_griseus1.0 | GCA_000419365.1 | 2360130144 | 52710 | 689430 | 1558295 | 2955633 | 0.0249226 | 0.4560103 |
| CriGri_1.0 | GCA_000223135.1 | 2399786748 | 109152 | 287345 | 1147233 | 2671025 | 0.0340257 | 0.4590470 |
3.5 scaffold長のヒストグラム
# 各配列長(log10)
w.cgr <- log10(width(f.cgr))
w.c_gr <- log10(width(f.c_gr))
w.crigri <- log10(width(f.crigri))
## xlim
xlim.al <- range(w.cgr, w.c_gr, w.crigri)
## 透過色作成
adjcols <- adjustcolor(c("red", "blue", "green"), alpha.f = 0.2)
## ヒストグラム
par(mfrow=c(2,2))
hist(w.cgr, col=adjcols[1], breaks = 100, border=NA,
xlab = " log10 Length", main = "Cgr_1.0")
hist(w.c_gr, col=adjcols[2], breaks = 100, border=NA,
xlab = " log10 Length", main = "C_griseus1.0")
hist(w.crigri, col=adjcols[3], breaks = 100, border=NA,
xlab = "log10 Length",main = "CriGri_1.0")
hist(w.crigri, col=adjcols[3], breaks = 100, border=NA,
xlim=xlim.al, main="", xlab="log10 Length " )
hist(w.cgr, col=adjcols[1], breaks = 100, border=NA, add=T)
hist(w.c_gr, col=adjcols[2], breaks = 100, border=NA, add=T)
## ggplot
library(ggplot2)
source("~/pub/bin/r/myfunc.R")
dat <- cbind.fillna(list(w.cgr, w.c_gr, w.crigri))
names(dat) <- c("cgr", "c_gr", "crigri")
dat <- tidyr::gather(dat)
dat <- subset(dat, !is.na(dat$value))
ggplot(dat, aes(x=value, fill=key), colour=key) +
geom_histogram(binwidth=0.01, alpha=0.3, position="identity") +
theme_bw() +
facet_grid(key ~ .)
ggplot(dat, aes(x=value, fill=key ), colour=key) +
geom_histogram(binwidth=0.01, alpha=0.3, position="identity") +
theme_bw() 