Chapter 11 Day3 강의 정리
11.1 Class 1 - NCBI 서열 download
library(rentrez)
entrez_dbs()
covid_paper <- entrez_search(db="pubmed", term="covid19")
class(covid_paper)
names(covid_paper)
covid_paper$ids
covid_link <- entrez_link(db="all", id=covid_paper$ids, dbfrom="pubmed")
names(covid_link$links)
covid_link$links$pubmed_pubmed실제 서열 데이터 다운로드는 entrez_fetch 함수 이용. fetch 함수를 이용한 다운로드는 text 파일을 그대로 받아오기 때문에 이를 fasta 등 적절한 포멧의 파일로 저장하고 다시 읽어들여야 함.
katipo_search <- entrez_search(db="popset", term="Latrodectus katipo[Organism]")
entrez_db_summary("popset")
katipo_search$ids
katipo_summs <- entrez_summary(db="popset", id=katipo_search$ids)
names(katipo_summs)
names(katipo_summs$`41350664`)
COI_ids <- katipo_search$ids[c(2,6)]
trnL_ids <- katipo_search$ids[5]
COI <- entrez_fetch(db="popset", id=COI_ids, rettype="fasta")
trnL <- entrez_fetch(db="popset", id=trnL_ids, rettype="fasta")
write(COI, "COI.fasta")
write(trnL, "trnl.fasta")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Biostrings")
library(Biostrings)
coi <- readDNAStringSet("COI.fasta")
trnl <- readDNAStringSet("trnl.fasta")
coi
class(coi)Covid19 accession 번호 탐색 후 관련 서열 다운로드
covid19 <- entrez_fetch(db="nuccore", id="NC_045512", rettype="fasta")
write(covid19, "covid19.fasta")
covid19seq <- readDNAStringSet("covid19.fasta")
as.character(covid19seq)
covid19 <- entrez_fetch(db="nuccore", id="NC_045512", rettype="gb")
covid19
write(covid19, "covid19.gb")Genbank 형태의 데이터 다운로드.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("genbankr")
library(genbankr)
covid19gb <- readGenBank(file = "covid19.gb")
class(covid19gb)
methods(class="GenBankRecord")
cds(covid19gb)11.2 class 2 - sequence alignment
특정 패턴이 존재하는지 검사
x <- DNAString("A")
methods(class="DNAString")
?matchPattern
coi
matchPattern("ATG", coi[[1]])
vmatchPattern("ATG", coi) 두 서열의 비교를 위한 pairwiseAlignment
coi[[1]]
aln <- pairwiseAlignment(pattern = coi[[1]], subject = coi[[2]])
class(aln)
methods(class="PairwiseAlignmentsSingleSubject")
?PairwiseAlignmentsSingleSubjectalignment된 서열을 효율적으로 가시화 해주는 툴 중에 하나 DECIOHER의 BrowseSeqs 함수
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DECIPHER")
library(DECIPHER)
BrowseSeqs(aln)
alignedPattern(aln)
alignedSubject(aln)
alnseqs <- c(alignedPattern(aln), alignedSubject(aln))
class(alnseqs)
BrowseSeqs(alnseqs)
BrowseSeqs(alnseqs, colWidth=100)
BrowseSeqs(alnseqs, colWidth=100, patterns = "TCCTGCCCGGGGCCT")DECIPHER는 AlignSeq함수를 이용해서 MSA을 수행할 수도 있음
library(DECIPHER)
dbConn <- dbConnect(SQLite(), ":memory:")
Seqs2DB(coi, "XStringSet", dbConn, "coi")
BrowseDB(dbConn)
Seqs2DB(trnl, "XStringSet", dbConn, "trnl")
BrowseDB(dbConn)
dna <- SearchDB(dbConn, identifier="coi")
dbDisconnect(dbConn)
alignedcoi <- AlignSeqs(coi)
BrowseSeqs(alignedcoi, colWidth = 100)
class(alignedcoi)
conseq <- ConsensusSequence(alignedcoi)MSA 결과 각 서열간의 거리가 계산될 수 있고 이러한 값들은 clustering이나 phylogenetic tree를 그리는데 사용될 수 있음.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ggtree")
library(ggtree)
alignedcoi <- AlignSeqs(coi)
class(alignedcoi)
dm <- DistanceMatrix(alignedcoi)
dim(dm)
tmp <- dm[1:4,1:4]
tree <- IdClusters(dm, cutoff=10, method="NJ", showPlot=TRUE, type="dendrogram")
class(tree)
methods(class="dendrogram")
tree2 <- as.hclust(tree)
class(tree2)
methods(class="hclust")
tree3 <- as.phylo(tree2)
class(tree3)만약 ggtree 실행이 잘 되지 않는 경우 다음 사이트를 참고해서 재설치 진행할 수 있음. issues, ggtree github
remotes::install_github("YuLab-SMU/ggtree")
ggtree(tr = tree3)
library(ggtree)
library(ape)
tr <- rtree(10)
class(tr)
require(ape)
tr <- rtree(10)
ggtree(tr)
ggtree(tr)11.3 class 3 - genomic data with ranges
지놈스케일의 데이터를 분석할 때 IRanges 패키지를 사용. 바이오 데이터 분석에 필수로 사용되는 패키지이며 Bioconductor software 중 4위에 올라있음.
library(IRanges)
ir <- IRanges(start = c(1,3,5), end = c(3,5,7))
ir
class(ir)
methods(class="IRanges")
start(ir)
disjointBins(ir)기존에는 IRange 객체의 가시화를 위해 직접 함수를 만들어 썼으나 이제 ggbio 패키지를 이용하면 됨.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ggbio")
library(ggbio)
ir <- IRanges(start = c(1,3,6), end = c(4,5,7))
ir
autoplot(ir, aes(fill=width)) +
theme_bw()유전체 정보는 IRrange를 활용한 위치 정보만으로는 부족하며 다양한 high-throughput 데이터로부터 얻어지는 정도를 효율적으로 저장하고 관리하기 위해서 GenomicRanges라는 패키지를 사용함.
library(GenomicRanges)
gr <- GRanges(
seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
score = 1:10,
GC = seq(1, 0, length=10),
product = rep("Aa", 10))
gr
rep("Aa", 10)GenomicRange 객체는 Rle와 IRange 객체를 조합해서 생성.
library(IRanges)
x <- "GATTGCCCCCCTAG"
y <- unlist(strsplit(x, split=""))
#unlist(strsplit(x, split=""))
y
yrle <- Rle(y)
yrle
runLength(yrle)
runValue(yrle)
nrun(yrle)
x <- Rle(values = c(1:3), lengths = c(1:3))
class(x)Meta data는 mcols 함수를 이용해서 관리
gr <- GRanges(
seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
score = 1:10,
GC = seq(1, 0, length=10))
gr
class(gr)
methods(class="GRanges")
seqnames(gr)
range(gr)
tmp <- mcols(gr)
tmp$GC tidyverse의 dplyr 패키지와 유사한 기능의 plyranges 패키지를 이용하면 genomicranges 데이터를 효율적으로 분석 가능.
11.3.1 Exercise
- 이번 시간 배운 일련의 패키지들을 이용해서 covid19 genbank 파일을 읽고 cds 서열을 추출한 후 autoplot을 이용해서 그래프를 그리시오.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("plyranges")
library(plyranges)
library(genbankr)
library(DECIPHER)
library(ggbio)
covid19 <- readGenBank("covid19.gb")
class(covid19)
covid19cds <- cds(covid19)
covidseq <- covid19@sequence
ranges(covid19cds)
?getSeq
covidcds <- getSeq(covidseq, covid19cds)
BrowseSeqs(covidcds, colWidth = 200)
names(covidcds) <- mcols(covid19cds)$product
BrowseSeqs(covidcds, colWidth = 100)
autoplot(covid19cds)- gene과product를 선택하고 gc ratio를 계산해서 metadata에 넣은 후 임의의 그룸으로 나누어 평균 gc 비율을 구하시오
library(plyranges)
gcr <- rowSums(letterFrequency(covidcds, letters = c("C", "G"), as.prob = T))
tmpd <- covid19cds %>%
select(gene, product) %>%
mutate(gc = gcr) %>%
filter(grepl(pattern = "ORF", gene)) %>%
mutate(g = sample(c("A", "B"), 9, replace=T))
gcmean <- tmpd %>%
group_by(g) %>%
summarise(m = mean(gc))png 함수를 사용하면 특정 plot을 그린 후 png 파일로 저장할 수 있음. dev.off까지 실행해 주어야 하며 블록으로 설정 후 실행해 주어야함 (한줄 한줄 실행하면 저장되지 않음).
png("covid19.png", width=10, height=5, units='in', res=300)
autoplot(tmpd, aes(fill=g)) +
theme_bw()
dev.off()11.4 class 4 - ORFfinder with docker
11.4.1 Exercise
- 리눅스 실행이 필요한 경우 활용할 수 있는 docker를 활용해 ORFfinder를 실행해보고 앞에서 배운 방법들을 이용해 covid19 지놈의 예측된 ORF를 Genomic range 형태로 만들고 가시화 하시오
파일브라우저로 작업 디렉토리에 (workding directory) 마우스 커서를 올리고 Shift + 마우스 오른쪽 버튼을 누르면 메뉴가 뜨고 그 중 Powershell 실행. 명령 프롬프트 상에서 아래와 같은 명령 실행. 참고로 docker는 사전에 설치되어 있어야 하며, unlhcc/orffinder 라는 아이디/프로그램파일명은 구글이나 docker hub 에서 orffinder로 검색해서 찾아볼 수 있음.
이미지 다운로드
docker pull unlhcc/orffinder다운로드 받은 이미지 리스트 확인
docker imagesORFfinder 실행
docker run --rm -v ${PWD}:/app -w /app unlhcc/orffinder:latest ORFfinder -in covid19.fasta -out covid19orf.fasta -outfmt 1docker run --rm -v ${PWD}:/app -w /app unlhcc/orffinder:latest 여기까지는 도커 컨테이너를 생성하는 부분이고 그 이후는 ORFfinder를 실행하는 명령어임. ${PWD} 는 현재 directory (input file인 covid19.fasta가 있는 디렉토리)를 나타냄. 위 명령어를 실행하면 covid19orf.fasta가 생성되며 다음 코드로 genomic range 데이터를 생성함.
covidorfs <- readDNAStringSet("covid19orf.fasta")
BrowseSeqs(covidorfs, colWidth = 100)
## average length
mydat <- data.frame(len=nchar(covidorfs))
ggplot(mydat, aes(x=len)) +
geom_histogram(bins=300) +
xlim(c(0, 2000))
## extract start end position
tmps <- names(covidorfs)
strsplit(tmps[1], split=":")
tmpl <- strsplit(tmps, split=":")
startpos <- lapply(tmpl, function(x){c(x[3])}) %>% unlist %>% as.numeric
endpos <- lapply(tmpl, function(x){c(x[4])}) %>% unlist %>% as.numeric
orfnames <- lapply(tmpl, function(x){c(x[2])}) %>% unlist
strnd <- rep("+", length(startpos))
strnd[startpos>endpos] <- "-"
## swap startpos and endpos
tmpp <- startpos
startpos[startpos>endpos] <- endpos[startpos>endpos]
endpos[tmpp>endpos] <- tmpp[tmpp>endpos]
## generate grange
covidorfir <- IRanges(start=startpos, end=endpos)
covidorfgr <- GRanges(seqnames = "covid19", covidorfir, strnd, orfnames)
autoplot(covidorfgr, aes(fill=strand))