Chapter 9 Working with DNA sequences II

9.1 Sequences from NCBI

전세계 연구자들이 서열 데이터를 분석하는데 가장 많이 이용하는 사이트 중 하나가 NCBI 이며 따라서 NCBI에서는 연구자들이 데이터베이스에 접근하기위한 편리한 방법을 제공하고 있고 그 중 하나가 Entrez 입니다.

R에서도 Entrez 기능을 도입한 package들이 제공되고 있으며 그 중 하나가 rentrez 입니다. https://www.ncbi.nlm.nih.gov/books/NBK25500/ 이 곳의 Downloading Full Records 를 참고하시면 좋습니다. Entrez는 대략적으로 다음 9개의 유틸리티를 제공합니다.

EInfo (database statistics)
ESearch (text searches)
EPost (UID uploads)
ESummary (document summary downloads)
EFetch (data record downloads)
ELink (Entrez links)
EGQuery (global query)
ESpell (spelling suggestions)
ECitMatch (batch citation searching in PubMed)

이 중 ESerach, EPost, ESummary, EFetch 등이 많이 사용하는 유틸이며 정보를 다운로드 받을 경우는 EFetch 를 주로 사용하게 됩니다. rentrez 는 위와 같은 NCBI Eutils API를 활용하여 R 환경에서 탐색이나 다운로드 등 NCBI 데이터베이스와 상호작용이 용이하도록 만들어 놓은 tool 입니다.

library(rentrez)

entrez_dbs()
entrez_db_summary("nuccore")


covid_paper <- entrez_search(db="pubmed", term="covid19")
covid_paper$ids

names(covid_paper)
covid_paper$ids


covid_link <- entrez_link(db="all", id=covid_paper$ids, dbfrom="pubmed")
names(covid_link)
names(covid_link$links)
head(covid_link$links$pubmed_pubmed)

특정 균주에 대한 정보를 찾은 후 두 개의 loci에 대한 서열 정보를 다운로드 하는 코드입니다. rettype (return type) 등 자세한 정보는 Eutils table 또는 NCBI Eutils 페이지를 참고하시기 바랍니다.

# popset database is a collection of related DNA sequences derived from population
katipo_search <- entrez_search(db="popset", term="Latrodectus katipo[Organism]")
katipo_search$ids

katipo_summs <- entrez_summary(db="popset", id=katipo_search$ids)
names(katipo_summs)
katipo_summs$`41350664`
class(katipo_summs)
methods(class="esummary_list")

titles <- extract_from_esummary(katipo_summs, "title")
unname(titles)

print(katipo_summs)
katipo_summs$`1790798044`$gi


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")

#library(Biostrings)
coi <- readDNAStringSet("COI.fasta")
trnl <- readDNAStringSet("trnl.fasta")

9.1.1 Exercise

  1. 뎅기바이러스 서열 4종에 대한 NCBI의 accession 번호가 다음과 같음 NC_001477, NC_001474, NC_001475, NC_002640 해당 DNA 서열을 fasta 형식으로 nuccore 데이터베이스에서 다운로드 하시오

  2. COVID-19 서열의 NCBI accession 번호를 찾고 nuccore 데이터베이스에서 fasta 포멧과 genbank 포멧의 정보를 다운로드 하고 파일에 저장하시오. 또한 이 파일들을 각각 Biostrings 패키지와 genbankr 패키지를 사용해서 읽어들이시오.

9.2 Pattern matching

Biostrings 패키지에는 하나의 subject 서열에 특정 pattern이 존재하는지 탐색하는 matchPattern함수를 제공합니다. 만약 여러개의 subject 서열에서 하나의 pattern을 찾을 경우에는 vmatchPattern함수를 사용하고 하나의 subject 서열에 여러개의 pattern을 찾는 경우에는 matchPDict 함수를 사용합니다.

library(Biostrings)

length(coi)
hits <- matchPattern("ATG", coi[[1]], min.mismatch=0, max.mismatch=0)
hits
class(hits)
methods(class="XStringViews")
ranges(hits)

hits <- vmatchPattern("ATG", coi, min.mismatch=0, max.mismatch=0)
stack(hits)

9.3 Align two sequences

Biostrings 패키지에는 다음과 같이 local, global alignment를 수행할 수 있는 함수를 제공하고 있습니다. 첫 번째 파라메터는 pattern이며 두 번째는 subject 로서 pattern은 query로서 해당 서열이 subject (target)에 있는지를 보는 것과 같습니다.

coi <- readDNAStringSet("COI.fasta")
coi

aln <- pairwiseAlignment(coi[[1]], coi[[2]])
alnseqs <- c(alignedPattern(aln), alignedSubject(aln))
class(aln)
class(alnseqs)

methods(class="PairwiseAlignmentsSingleSubject")
methods(class="DNAStringSet")

library(DECIPHER)
BrowseSeqs(alnseqs)
BrowseSeqs(alnseqs, colWidth=200)
BrowseSeqs(alnseqs, colWidth=200, patterns = "TCCTGCCCGGGGCCT")

DECIPHER 패키지는 서열 alignment나 primer design 등을 수행할 수 있는 패키지로 다음과 같이 별도 메모리에 서열을 저장하고 빠르게 alignment를 수행할 수 있어서 중소 규모의 서열에 대한 분석으로 유용하게 사용될 수 있습니다.

library(DECIPHER)
dbConn <- dbConnect(SQLite(), ":memory:")
Seqs2DB(coi, "XStringSet", dbConn, "coi")
BrowseDB(dbConn)

l <- IdLengths(dbConn)
Add2DB(l, dbConn)
BrowseDB(dbConn)

Seqs2DB(trnl, "XStringSet", dbConn, "trnl")
BrowseDB(dbConn)

## extract sequences
dna <- SearchDB(dbConn, identifier="coi")
BrowseSeqs(dna)

dbDisconnect(dbConn)

9.4 Multiple sequence alignment

Multiple sequence alignment(MSA) tool은 서열 데이터의 양과 계산량의 문제로 linux 기반 commandline 프로그램들이 많습니다. 대표적으로 CLUSTAL-Omega, MUSCLE. window 기반 환경에서는 docker 등을 활용해서 관련 분석을 수행할 수 있습니다. 본 강의에서는 DECIPHER 패키지를 활용합니다.

library(Biostrings)
library(DECIPHER)

coi <- readDNAStringSet("COI.fasta")
BrowseSeqs(coi)
alignedcoi <- AlignSeqs(coi)
BrowseSeqs(alignedcoi)
class(alignedcoi)

conseq <- ConsensusSequence(alignedcoi)
IUPAC_CODE_MAP

9.5 Phylogenetic trees with clustering

dm <- DistanceMatrix(alignedcoi)
class(dm)
dim(dm)
dm[1:2,1:2]

tree <- IdClusters(dm, cutoff=10, method="NJ", showPlot=TRUE, type="dendrogram")
class(tree)
methods(class="dendrogram")
plot(tree)

트리는 ggplot 형태의 ggtree, reference를 사용하면 쉽게 그릴 수 있습니다.

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ggtree")

library(ggtree)
tree <- rtree(n = 20)
class(tree)           
methods(class="phylo")

library("ape")

## convert to dendrogram -> hclust -> phylo 
cl <- as.hclust(tree)
class(cl)
methods(class="hclust")

py <- as.phylo(cl)
class(py)
ggtree(py)

ggtree(py, layout="circular")

ggtree(py, layout="circular") +
  geom_tiplab(size=1, aes(angle=angle))