18 Lecture 08 (0811)

class 01

• rentrez 패키지를 사용한 NCBID 데이터 다운로드

library(rentrez)

entrez_dbs()
entrez_db_summary("nuccore")

covid_paper <- entrez_search(db="pubmed", term="covid-19", retmax=40)
class(covid_paper)
methods(class="esearch")
names(covid_paper)
covid_paper$ids

katipo_search <- entrez_search(db="popset", term="Latrodectus katipo[Organism]")
class(katipo_search)
names(katipo_search)
katipo_search$ids

entsum <- entrez_summary(db="popset", id = katipo_search$ids)
entsum[[2]]$title
entsum[[2]]$article

extract_from_esummary(entsum, "title")

• 텍스트 형태로 다운로드 받으므로 파일에 쓰고 readDNAString 형태로 다시 읽어옴

library(Biostrings)

coi <- entrez_fetch(db="popset", id=katipo_search$ids, rettype="fasta")
write(coi, file = "coi.fasta")
coiseq <- readDNAStringSet(filepath = "coi.fasta")

• 뎅기바이러스 예제

dv <- c("NC_001477", "NC_001474", "NC_001475", "NC_002640")
dvseq <- entrez_fetch(db="nuccore", id=dv, rettype = "fasta")
dvseq

write(dvseq, file="dvseq.fasta") 
dvseq <- readDNAStringSet("dvseq.fasta")  
dvseq

strwrap(as.character(dvseq[[1]]), width=10)

• Covid-19 예제

sresult <- entrez_search(db="popset", term="Covid-19", retmax=40)
class(sresult)
names(sresult)
sresult$ids

mysummary <- entrez_summary("popset", sresult$ids)
mysummary

#lapply
#sapply

mysummary
extract_title <- function(x){
  return(x$title)
}
mysummary[[1]]$title
extract_title(x=mysummary[[1]])


sapply(mysummary, extract_title)
tmp <- extract_from_esummary(mysummary, "title")
tmp
class(tmp)
mydata <- data.frame(title=tmp)
mydata


## sequence download

tmpstr <- entrez_fetch("popset", id=sresult$ids, rettype = "fasta")
write(tmpstr, "covid19.fasta")

covid19seq <- readDNAStringSet("covid19.fasta")
covid19seq
names(covid19seq)

class 02

• popset 데이터베이스에서 선택한 각 아이템 (스터디별) 사용한 샘플 수 구하기
• sapply 및 사용자 정의 함수 extract_statistics_count 사용한 barplot 그리기

library(tidyverse)

mydatadf <- mydata %>% 
  rownames_to_column(var = "uid")

mysummary$`2279969783`$statistics$count[4]
mysummary$`2271163415`$statistics$count[4]
mysummary$`2252835834`$statistics$count[4]
## ..

sapply(mysummary, extract_statistics_count)


## build function
extract_statistics_count <- function(x){
  ## x == mysummary$`2279969783`
  z <- x$statistics$count[4]
  return(z)
}

cnt <- sapply(mysummary, extract_statistics_count)
cntdf <- data.frame(cnt) %>% 
  rownames_to_column("uid")

# mydata %>% 
#   rownames_to_column(var = "uid") %>% 
#   left_join(y=cnt, by="uid")

mydatadf
cntdf

mydata2 <- left_join(y=mydatadf, x=cntdf, by="uid") 


## barplot

ggplot(mydata2, aes(x=uid, y=cnt)) +
  geom_bar(stat = "identity") +
  theme(
    axis.text.x = element_text(angle=90)
  ) +
  ylab("Count")

• 코로나 바이러스 ID들에 대해서 정렬을 수행하고 writePairwiseAlignments를 이용한 가시화

covid <- data.frame(
species = c(rep("Human", 7), c("Civet", "Civet"), rep("Bat", 3), "Pangolin"),
coronavirus = c("SARS-CoV-2", "SARS-CoV-2", "SARS-CoV-1", "SARS-CoV-1", "SARS-CoV-1", "H-CoV-OC43", "MERS-CoV", "SARS-CoV", "SARS-CoV", "SL-CoV", "SL-CoV", "SL-CoV", "SL-CoV"),
isolate = c("Wuhan Hu-1", "USA-WA-1", "Urbani", "Tor2", "GD03T10013", "UK/London",  "EMC-2012", "SZ3", "Civet007", "ZXC21", "WIV16", "RaTG13", "MP789"),
year = c("2020", "2020", "2002", "2002", "2003", "2011", "2011", "2003", "2004", "2015", "2013", "2013", "2020"),
gbacc = c("NC_045512.2", "MN985325.1", "AY278741.1", "AY274119.3", "AY525636.1", "KU131570.1", "NC_019843.3", "AY304486.1", "AY572034.1", "MG772934.1", "KT444582.1", "MN996532.1", "MT084071.1"))
write.csv(covid, file = "covid_table.csv", quote = F, row.names=F)


covid19info <- read.csv("covid_table.csv")

str1 <- entrez_fetch("nuccore", covid19info$gbacc, rettype = "fasta")
str2 <- entrez_fetch("nuccore", covid19info$gbacc, rettype = "gb")
write(str1, "covid19seqs.fasta")
write(str2, "covid19seqs.gb")

library(Biostrings)

covid19seqs <- readDNAStringSet("covid19seqs.fasta")
covid19seqs

aln <- pairwiseAlignment(covid19seqs[1], covid19seqs[2])
aln

writePairwiseAlignments(aln, block.width=50)
writePairwiseAlignments(aln, file = "covidaln.txt", block.width=50)

summary(aln)

consensusMatrix(aln)[,1:10]
consensusString(aln)

• (참고) 문자열 나누기 코드, sapply와 사용자 정의 함수 extract_first_element 사용법 익히기

covid19seqs
names(covid19seqs)

strsplit("AT GC", split=" ")

strsplit(names(covid19seqs)[1], split=" ")
strsplit(names(covid19seqs)[1:2], split=" ")

tmpstr <- strsplit(names(covid19seqs), split=" ")
tmpstr[1:5]
sapply(tmpstr, extract_first_element)

extract_first_element <- function(x){
  return(x[1])
}
myacc <- sapply(tmpstr, extract_first_element)
mytitle <- sapply(tmpstr, function(x){
  z <- paste(x[-1], collapse=" ")
  return(z)
  })

mydf <- data.frame(myacc, mytitle)
names(covid19seqs) <- myacc
covid19seqs

aln <- pairwiseAlignment(covid19seqs[1], covid19seqs[2])
writePairwiseAlignments(aln, block.width=50)

class 03

• DECIPHER 패키지 활용한 다중 서열 정렬

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

BiocManager::install("DECIPHER")

library(DECIPHER)

• DECIPHER 패키지 서열 관리

dbConn <- dbConnect(SQLite(), ":memory:")
Seqs2DB(covid19seqs, "XStringSet", dbConn, "covid19")
BrowseDB(dbConn)

coi <- readDNAStringSet("coi.fasta")
Seqs2DB(coi, "XStringSet", dbConn, "coi")
BrowseDB(dbConn)

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

covid19seq2 <- SearchDB(dbConn, identifier = "covid19")

dbDisconnect(dbConn)
BrowseDB(dbConn)

• DECIPHER 패키지 AlignSeqs 활용한 서열정렬 및 BrowseSeqs 이용한 가시화

covid19sub <- subseq(covid19seq2, 1, 2000)
BrowseSeqs(covid19seq2, colWidth = 80, patterns=DNAStringSet(c("ACTG", "CSC")))


aln2 <- AlignSeqs(covid19sub)

BrowseSeqs(aln2, colWidth = 80)
BrowseSeqs(aln2, colWidth = 80, patterns=DNAStringSet(c("ACTG", "CSC")))
BrowseSeqs(aln2, colWidth = 80, patterns="-", colors="black")

• RESTRICTION_ENZYMES 활용 예제

data(RESTRICTION_ENZYMES)
RESTRICTION_ENZYMES

rsite <- RESTRICTION_ENZYMES["BsaI"]
d <- DigestDNA(rsite, covid19seq2[1])
unlist(d)

pos <- DigestDNA(rsite, covid19seq2[1], type="positions")
unlist(pos)

BrowseSeqs(covid19seq2[1], colWidth = 100, patterns=rsite)
library(stringr)
str_extract_all(rsite, "[[A-Z]]", simplify = T)


rsite2 <- paste(str_extract_all(rsite, "[[A-Z]]", simplify = T), collapse="")
rsite3 <- as.character(reverseComplement(DNAString(rsite2)))
BrowseSeqs(covid19seq2[1], colWidth = 100, patterns=c(rsite2, rsite3))

    아래 관련 알아보기 
    - DECIPHER AA 적용 여부
    - DECIPHER --> ggtree

• 웹에서 BLAST 실행 후 결과 데이터 분석

myhit <- read.csv("F9BA23TY01R-Alignment-HitTable.csv", header = F)
myhit
myseq <- readDNAStringSet("seqdump.txt")
myseq

myaln <- AlignSeqs(myseq)
#AAStringSet(myaln)
BrowseSeqs(myaln)

myconm <- consensusMatrix(myaln)
mycons <- consensusString(aln)
dim(myconm)
myconm[1:10, 1:10]

class 04

• ggtree 활용

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

BiocManager::install("ggtree")
BiocManager::install("tidytree")
BiocManager::install("genbankr")

library(ggtree)
library(tidytree)


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

ggplot(tree) +
  geom_tree() +
  theme_tree()
  

ggtree(tree, branch.length="none")

ggtree(tree, layout="circular") +
  geom_tiplab(size=3, aes(angle=angle))

ggtree(tree, layout="circular", branch.length="none") +
  geom_tiplab(size=3, aes(angle=angle))

ggtree(tree) +
  theme_tree2()


ggtree(tree) +
  theme_tree2() +
  geom_tiplab() +
  geom_hilight(node=23, fill="steelblue", alpha=.4) 

as_tibble(tree) ## tidytree


ggtree(tree) +
  theme_tree2() +
  geom_tiplab() +
  geom_hilight(node=23, fill="steelblue", alpha=.4) 


dat <- data.frame(node=c(23, 35), type=c("AA", "BB"))

ggtree(tree) +
  theme_tree2() +
  geom_tiplab() +
  geom_hilight(dat, aes(node=node, fill=type), alpha=.4) +
  scale_fill_manual(values=c("steelblue", "darkgreen"))

• phylo 클래스 데이터를 tibble 형태로 변환 후 데이터 참고

as_tibble(tree) %>% 
  filter(label==c("t2", "t4"))

• genbank 타입 파일에서 정보 읽어오기 (다음시간 계속)

library(rentrez)
require(genbankr)

mygb <- entrez_fetch("nuccore", id="NC_045512.2", rettype = "gb")
write(mygb, "co.gb")

mygbinfo <- readGenBank("co.gb")
class(mygbinfo)
methods(class="GenBankRecord")
chrseq <- getSeq(mygbinfo)
cdsrng <- cds(mygbinfo)
cdsrng
cdsseq <- getSeq(chrseq, cdsrng)
cdsseq

18.1 Lecture 09 (0901)

class 01

• msa
• NC_045512.2 NCBI 다운로드

library(Biostrings)
library(rentrez)

entrez_dbs()
tmps <- entrez_fetch(db="nuccore", id = "NC_045512.2", rettype = "fasta")
tmps
write.table(tmps, "covid19yu.fasta", quote = F, row.names = F, col.names = F)

download.file(url = "https://raw.githubusercontent.com/greendaygh/kribbr2022/main/covid_table.csv", destfile = "covid_table2.csv")

covid19 <- read.csv("covid_table2.csv")
selid <- covid19$gbacc[1:4]

tmps <- entrez_fetch(db="nuccore", id = selid, rettype = "fasta")
write.table(tmps, "covid19.fasta", quote = F, row.names = F, col.names = F)

• msa 설치
• Bioconductor landing page (구글에서 bioconductor msa)
• installation 코드 복사 후 콘솔에서 붙여넣기

## install 
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("msa")


library(msa)


mycovidseq <- readDNAStringSet("covid19.fasta") 
class(mycovidseq)
mycovidseq

• 일부 서열만 가지고 실습
• (팁) 블럭지정 후 ctrl + enter 블럭만 실행

#subseq
#substr


mycovidsubseq <- subseq(mycovidseq, 1, 2000)

aln <- msa(mycovidsubseq, method="Muscle")
aln
DNAStringSet(aln)

?msa

print(aln, show="complete")
#msaPrettyPrint(aln, output = "pdf")

alphabetFrequency(aln)

• IRanges 설치Bioconductor landing page (구글에서 bioconductor IRanges 검색 ) installation 코드 복사 후 콘솔에서 붙여넣기

library(IRanges)
nchar(consensusString(aln))

mym <- IRanges(start=c(1, 73), end=c(36, nchar(consensusString(aln))))
colmask(aln) <- mym
print(aln, show="complete")

alphabetFrequency(aln)
unmasked(aln)

• DECIPHER

library(DECIPHER)

aln <- AlignSeqs(mycovidsubseq)
BrowseSeqs(aln)
BrowseSeqs(aln, colWidth = 100)

class 02

• ggtree (bioconductor), treeio (Bioconductor), ape (CRAN) 설치
• 1. MSA - 2) Distance - 3) clustering - 4) Tree

aln <- msa(mycovidsubseq, method="Muscle")
class(aln)

?stringDist

mydist <- stringDist(DNAStringSet(aln))
class(mydist)

• label 축소
• treeio error 발생시, ggtree 패키지 unload (detach(“package:ggtree”, unload = TRUE)) 후에 treeio, ggtree 순서로 로딩

library(treeio)
library(ggtree)


tmps <- strsplit(names(mycovidsubseq), split=" ")
mynames <- sapply(tmps, function(x){x[1]})
names(mycovidsubseq) <- mynames
mycovidsubseq

aln <- msa(mycovidsubseq)
mydist <- stringDist(DNAStringSet(aln))
myclust <- hclust(mydist)
plot(myclust)

myphylo <- as.phylo(myclust)

ggtree(myphylo)  +
  geom_tiplab(color="firebrick")

ggtree(myphylo, layout="circular")  +
  geom_tiplab(color="firebrick")

• DECIPHER 버전업 되면서 바뀐부분 체크

aln <- AlignSeqs(mycovidsubseq)
#mydist <- DistanceMatrix(aln)
# myclust <- IdClusters(aln, cutoff=10, method="NJ", type="dendrogram")
# class(mydist)

• as.tibbe

library(tidyverse)

as.tibble(myphylo)

• 여러 그림 그리기

mytree <- rtree(30)
class(mytree)

p <- ggtree(mytree) +
  geom_tiplab()

p

a <- runif(30, 0, 1)
b <- 1 - a
df <- data.frame(mytree$tip.label, a, b)
df2 <- pivot_longer(df, -mytree.tip.label)

p2 <- p + geom_facet(panel="bar", 
                     data=df2, 
                     geom=geom_bar, 
                     mapping = aes(x = value, fill = as.factor(name)),
                     orientation = 'y', 
                     width=0.8,
                     stat = 'identity') +
  xlim_tree(9)

facet_widths(p2, widths = c(1, 2))

• Genome 분석
  
• NC_045512.2 NCBI 다운로드

library(Biostrings)
library(rentrez)


tmps <- entrez_fetch(db="nuccore", id = "NC_045512.2", rettype = "fasta")
tmps
write.table(tmps, "covid19yu.fasta", quote = F, row.names = F, col.names = F)

covid19seq <- readDNAStringSet("covid19yu.fasta")

• read genbank

require(genbankr)

tmps <- entrez_fetch(db="nuccore", id = "NC_045512.2", rettype = "gb")
tmps
write.table(tmps, "covid19yu.gb", quote = F, row.names = F, col.names = F)

covidgb <- readGenBank("covid19yu.gb")
covidgb

cds(covidgb)
class(covidgb)
methods(class="GenBankRecord")

covidseq <- getSeq(covidgb)
class(covidgb)

cdsseq <- getSeq(covidseq, cds(covidgb))
class(covidseq)

class 03

library(IRanges)

ir <- IRanges(start = 1, end = 10)
ir

ir <- IRanges(start = 1:10, width = 10:1)
ir
class(ir)

x <- "GATTGCCCCCCTAG"
y <- unlist(strsplit(x, split=""))
z <- Rle(y)
write.table(z, "tmp.txt")

cds(covidgb)

library(ggbio)

ir <- ranges(cds(covidgb))
ir

start(ir)
end(ir)
width(ir)

autoplot(ir) +
  theme_bw()

autoplot(ir, aes(fill=width)) +
  theme_bw()




ir <- IRanges(c(1, 10, 100), end=c(20, 50, 110))
autoplot(ir)


disjointBins(ir)
autoplot(disjoin(ir))
autoplot(IRanges::reduce(ir))

• ggbio 설치

IRanges(start = c(1,1,1,1), width = c(100, 15, 30, 4))

library(GenomicRanges)
cds(covidgb)

gr <- GRanges(seqnames = "chr1", 
              ranges = IRanges(1:10, 21:30),
              strand = "+")
gr


gr <- GRanges(seqnames = "chr1", 
              ranges = IRanges(1:10, 21:30),
              strand = "+", 
              names = paste0("g", 1:10), 
              score = 1:10)
gr

seqnames(gr)
ranges(gr)
strand(gr)
start(gr)
end(gr)

granges(gr) 
mcols(gr)
seqlengths(gr) <- 50

autoplot(gr)


gr2 <- GRanges(seqnames = "chr2", 
              ranges = IRanges(1:10, 21:30),
              strand = "+", 
              names = paste0("g", 11:20), 
              score = 1:10)
gr3 <- c(gr, gr2)
seqnames(gr3)
seqlengths(gr3)

split(gr3, c(rep("a", 10), rep("b", 10)))

• disjoin
• reduce

class 04

• 대장균 지놈 plot

tmps <- entrez_fetch("nuccore", id="NC_000913.3", rettype="gbwithparts")
write(tmps, "ecoli-mg1655.gb")

ecoligb <- readGenBank("ecoli-mg1655.gb")
ecoli_cds <- cds(ecoligb)
autoplot(ecoli_cds)
seqlengths(ecoli_cds)
ggbio() +
  circle(ecoli_cds, geom="ideo", fill = "gray70") +
  circle(ecoli_cds, geom="scale", size=3) +
  circle(ecoli_cds, geom = "text", aes(label = locus_tag), vjust = 0, size = 3) +
  theme(
    axis.text.x = element_text(angle=90)
  )


gr1 <- granges(ecoli_cds)
gr2 <- granges(ecoli_cds)
mcols(gr2)$test <- rnorm(length(ecoli_cds))

ggplot() + 
  layout_circle(ecoli_cds, geom = "ideo", fill = "gray70", radius = 9, trackWidth = 1) +
  layout_circle(ecoli_cds, geom = "scale", size = 3, trackWidth = 1, scale.n=20) +
  layout_circle(gr1, geom = "rect", color = "steelblue",  radius = 5) +
  layout_circle(gr2, geom = "bar", aes(y=test), size = 3, trackWidth = 1, scale.n=20, radius = 4) 
  

gr1disj <- disjoin(gr1)

ggplot() + 
  layout_circle(ecoli_cds, geom = "ideo", fill = "gray70", radius = 9, trackWidth = 1) +
  layout_circle(ecoli_cds, geom = "scale", size = 3, trackWidth = 1, scale.n=20) +
  layout_circle(gr1disj, geom = "rect", color = "steelblue",  radius = 5) +
  layout_circle(gr2, geom = "bar", aes(y=test), size = 3, trackWidth = 1, scale.n=20, radius = 4) 

library(plyranges)

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

BiocManager::install("plyranges")

covidseq <- getSeq(covidgb)
cdsrng <- cds(covidgb)
cdsseq <- getSeq(covidseq, cdsrng)

gcratio <- rowSums(letterFrequency(cdsseq, c("G", "C"))) / nchar(cdsseq)

cdsrng %>% 
  dplyr::select(gene, product) %>% 
  mutate(gcratio)
  

• gggene install –> CRAN
• Rstudio > Packages > install

library(gggenes)

mydf <- data.frame(seqname = seqnames(cdsrng),
           start = start(cdsrng), 
           end = end(cdsrng), 
           strand = case_when(
             as.vector(strand(cdsrng))=="+"~ TRUE,
             as.vector(strand(cdsrng))=="-"~ FALSE
           ),
           gene = mcols(cdsrng)$gene)

ggplot(mydf, aes(xmin = start, 
                 xmax = end, 
                 y = seqname,
                 label = gene, 
                 fill = gene,
                 forward = strand)) +
  geom_gene_arrow(arrowhead_height = grid::unit(12, "mm"),
    arrowhead_width = grid::unit(6, "mm"),
    arrow_body_height = grid::unit(12, "mm")) +
  geom_gene_label() +
  theme(legend.position = "none")
      

18.2 Lecture 10 (0915)

class 01

명령창에서 fasterq-dump 실행시 에러가 발생할 경우 아래 실행

vdb-config --interactive

Enable Remote Access 가 (X)로 선택되어 있으면 save 후 exit. tab 키로만 움직일 수 있음. Space 키가 선택 (ok).

• Cmd 창에서 dir 입력하면 파일 리스트가 보임

prefetch --

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

BiocManager::install("GEOquery")

• GDS

library(GEOquery)
.libPaths()

gds <- getGEO(filename = system.file("extdata/GDS507.soft.gz",package="GEOquery"))
gds  
class(gds)
methods(class="GDS")
Table(gds)
Columns(gds)

• GSM

gsm <- getGEO(filename=system.file("extdata/GSM11805.txt.gz",package="GEOquery"))
methods(class=class(gsm))
head(Meta(gsm))
Table(gsm)
Columns(gsm)

• GPL

gpl <- getGEO(filename=system.file("extdata/GPL97.annot.gz",package="GEOquery"))
class(gpl)
methods(class="GPL")
Table(gpl)

• GES class

gse <- getGEO(filename=system.file("extdata/GSE781_family.soft.gz",package="GEOquery"))
methods(class=class(gse))
Meta(gse)
head(GSMList(gse))
gsm <- GSMList(gse)[[1]]
Meta(gsm)
Table(gsm)
Columns(gsm)

class 02

• ExpressionSet

myexp <- Table(gds)
mysample <- Columns(gds)
mygene <- Table(gpl)

myeset <- list()
myeset$exp <- myexp
myeset$sample <- mysample
myeset$gene <- mygene

myeset

• Download GSE

gse2553 <- getGEO('GSE2553', GSEMatrix=TRUE)
gse2553
class(gse2553)
class(gse2553[[1]])
myeset <- gse2553[[1]]

exprs(myeset)
fData(myeset)
pData(myeset)

gds <- getGEO(filename=system.file("extdata/GDS507.soft.gz",package="GEOquery"))
class(gds)
methods(class="GDS")
GPL(gds)
eset <- GDS2eSet(gds, do.log2=TRUE)
eset
fData(eset)

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

BiocManager::install("SummarizedExperiment")

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

BiocManager::install("airway")

library(airway)

data(airway)
airway

myexp <- assay(airway)
dim(myexp)

myfeature <- rowRanges(airway)
length(myfeature)
myfeature[[1]]


mysample <- colData(airway)
mysample

• SummarizedExpriment 생성

nrows <- 200
ncols <- 6 
counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows)
dim(counts)

rowRanges <- GRanges(rep(c("chr1", "chr2"), c(50, 150)),
                     IRanges(floor(runif(200, 1e5, 1e6)), width=100),
                     strand=sample(c("+", "-"), 200, TRUE),
                     feature_id=sprintf("ID%03d", 1:200))

colData <- DataFrame(Treatment=rep(c("ChIP", "Input"), 3),
                     row.names=LETTERS[1:6])

se <- SummarizedExperiment(assays = list(counts = counts),
                           rowRanges = rowRanges,
                           colData = colData
                           )
se

class 03

• 표준정규분포 모수:

library(tidyverse)

x <- seq(-4, 4, 0.01)
y <- dnorm(x, 0, 1)
?dnorm
dat <- data.frame(x, y)
ggplot(dat, aes(x=x, y=y)) +
  geom_point()

y2 <- dnorm(x, 0, 0.5)
dat <- data.frame(x, y, y2)

dat %>% 
  pivot_longer(-x) %>% 
  ggplot(aes(x=x, y=value, color=name)) +
  geom_point()

nsample <- 16
x <- rnorm(nsample, 0, 1)
y <- rep(1, nsample)
hist(x)
xbar <- mean(x)
dat <- data.frame(x, y)
ggplot(dat, aes(x, y)) +
  geom_point()+
  geom_point(aes(x=mean(x), y=1), color="blue", size=5, shape=15)

nsample <- 16
x <- rnorm(nsample*2, 0, 1)
y <- c(rep(1, nsample), rep(0.9, nsample))
g <- factor(c(rep(1, nsample), rep(2, nsample)))
dat <- data.frame(x, y, g)

ggplot(dat, aes(x, y)) +
  geom_point() +
  geom_point(aes(x=mean(x[1:nsample]), y=1), colour="blue", size=5, shape=15) +
  geom_point(aes(x=mean(x[(nsample+1):length(x)]), y=0.9), colour="blue", size=5, shape=15) +
  scale_y_continuous(limits=c(0, 1.2))

nsample <- 16
nrep <- 10

x <- rnorm(nsample*nrep, 0, 1)
tmpy <- seq(0.1, 1, length.out=nrep)
y <- rep(tmpy, each=nsample)
## ?rep
g <- factor(y)

dat <- data.frame(x, y, g)
head(dat)
## sample means
dat_mean <- dat %>% 
  group_by(g) %>% 
  summarise(mean=mean(x))
head(dat_mean)

ggplot() + 
  geom_point(data=dat, aes(x, y)) +
  geom_point(data=dat_mean, 
             aes(x=mean, y=as.numeric(as.character(g))), 
             colour="blue", 
             size=5, 
             shape=15) +
  theme_bw()
  

• 1000개 샘플링 해서 평균 값을 모집단의 평균값이다 라고 하는 것 보다 100개 샘플링 후 평균을 10번 반복 후 (샘플의) 평균의 분포를 가지고 모집단의 평균을 말하는 것이 더 정확하다

head(dat)
head(dat_mean)

x <- seq(-4, 4, by=0.01)
# distribution of the samples 
y <- dnorm(x, mean(dat$x), sd(dat$x))
# distribution of the sample means
y2 <- dnorm(x, mean(dat_mean$mean), sd(dat_mean$mean))
dat2 <- data.frame(x, y, y2)
dat2_long <- dat2 %>% 
  pivot_longer(cols=c(y, y2))
head(dat2_long)

ggplot(dat2_long, aes(x=x, y=value, color=name)) +
  geom_line(size=1.2) +
  labs(y="Density") +
  theme_bw() +
  scale_color_manual(name="Type", 
                     labels=c("Samples", "Sample Means"),
                     values=c("red", "blue"))

nsample <- 100
x <- rnorm(nsample, 100, 16)
hist(x)

z <- (90 - 100)/(10/sqrt(10))
z

x <- seq(-4, 4, 0.01)
y <- dnorm(x, 0, 1)
dat <- data.frame(x, y)
ggplot(dat, aes(x=x, y=y)) +
  geom_point() +
  geom_point(aes(x=z, y=0), size=5, color="blue")

class 04

• 모집단의 분포를 가정을 하고 샘플링을 통해서 통계량을 구한 후 해당 통계량이 가정한 모집단 분포의 어디에 위치하는지를 보고 유의성 판단

x1 <- seq(-6, 6, by=0.01)
y1 <- dnorm(x1, 0, 1)
z <- data.frame(x1, y1)
pval <- round(1-pnorm(2,0,1), 4)

ggplot(z) +
  geom_line(aes(x1, y1), color="purple") +
  geom_vline(xintercept = 2) +
  geom_hline(yintercept = 0) +
  geom_area(data=filter(z, x1 > 2), 
            aes(x1, y1), 
            fill="#80008055") +
  annotate("text", x=3, y=0.1, label=pval) 

x1 <- seq(-6, 6, by=0.01)
y1 <- dnorm(x1, 0, 1)
x2 <- seq(-1, 11, by=0.01)
y2 <- dnorm(x2, 3, 1)
z <- data.frame(x1, y1, x2, y2)
#pval <- round(1-pnorm(2,0,1), 4)
ggplot(z) +
  geom_line(aes(x1, y1), color="blue") +
  geom_line(aes(x2, y2), color="red") +
  geom_vline(xintercept = 2) +
  geom_hline(yintercept = 0) +
  geom_area(data=filter(z, x1 > 2), 
            aes(x1, y1), 
            fill="#0000ff55") +
  geom_area(data=filter(z, x2 < 2), 
            aes(x2, y2), 
            fill="#ff000055") +
  annotate("text", x=3, y=0.1, label=pval) 

• example

mpg <- c(11.4, 13.1, 14.7, 14.7, 15, 15.5, 15.6, 15.9, 16, 16.8)

tstat <- (mean(mpg)-17)/(sd(mpg)/sqrt(length(mpg)))
tstat

x <- seq(-5, 5, length=100)
y <- dt(x, df=length(mpg)-1)
dat <- data.frame(x, y)
ggplot(dat, aes(x, y)) +
  geom_line() +
  geom_vline(xintercept = tstat)

pt(tstat, df=9)

• DEG (lecture note 06, class 03) 분석 예제
• 아래 예제는 코드를 보충해서 본문으로 옮겼습니다

gds <- getGEO(filename = system.file("extdata/GDS507.soft.gz",package="GEOquery"))
gds
eset <- GDS2eSet(gds, do.log2=TRUE)
eset

myexp <- data.frame(exprs(eset))
myfeature <- fData(eset)
mypheno <- pData(eset)

table(mypheno$disease.state)

class(myexp)

## transformation
mydat1 <- myexp %>% 
  rownames_to_column() %>% 
  pivot_longer(cols = -rowname) %>% 
  pivot_wider(names_from = rowname, values_from = value)  

mydat2 <- mypheno %>% 
  dplyr::select(sample, disease.state) %>% 
  left_join(mydat1, by = c("sample" = "name"))


mymean <- mydat2 %>% 
  group_by(disease.state) %>% 
  summarise(across(where(is.numeric), mean))



mymean2 <- mymean %>% 
  pivot_longer(-disease.state) %>% 
  pivot_wider(names_from=disease.state, values_from = value)

ggplot(mymean2, aes(x=normal, y=RCC)) +
  geom_point()