R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Molekyylisystematiikkaa R:llä

Molekyylisystematiikka on tieteen ala, jossa tutkitaan lajien sukulaisuussuhteita geneettisten tuntomerkkien, kuten DNA-sekvenssien avulla. Alaa on perinteisesti leimannut työkalujen sekalaisuus, millä tarkoitan sitä, että tulosten saaminen on edellyttänyt kohtuullisen määrän eri ohjelmistojen opettelua ja niiden ominaisuuksien yhdistelyä.

Ohjelmistoja, jotka ovat jokseenkin yleisessä käytössä on kymmeniä, joten aiemmin toimin ihan tavanomaisesti esimerkiksi seuraavasti: 1) hain sekvenssit Genbank-tietokannasta selaimella, 2) rinnastin sekvenssit ClustalX, Muscle- ja/tai mafft-ohjelmalla, 3) muutin sekvenssirinnastuksen formattin eri ohjelmille sopivaksi seqin-ohjelmalla, 4) etsin kullekin rinnastukselle parsimoniapuut PAUP*-ohjelmalla, ja maximum likelihood-puut RAxML-ohjelmalla, 5) vertailun puita topologia-etäisyyden perusteella Phylip-ohjelmistolla ja shimodaira-hasegawan testillä PAML-ohjelmistossa, ja 6) lopuksi visualisoin puut TreeView- tai Mega-ohjelmistolla. Useimpia mainituista ohjelmistoista on edelleen tarpeen käyttää, mutta R:llä ne voi näppärästi sitoa yhteen.

Seuraavassa on esimerkki R-skriptistä, joka hoitaa edellä mainitut työvaiheet R:ssä tai kutsuu tarvittaessa muita ohjelmistoja komentoriviltä. Se olettaa, että ohjelmistot ovat työpöydällä polussa ”exepath”, ja aineistot työpöydällä polussa ”datapath”, mutta näitä voi muuttaa skriptin alussa olevilla polkuasetuksilla. Skripti hakee sekvenssit Genbank:sta automaattisesti tunnistenumeroilla, ajaa niille analyysin neljälle eri menetelmällä ja lopuksi vertaa eri menetelmien tuloksia toisiinsa.

# Automatic phylogenetic analysis
# JTT 2015-04-25
 
#################################
#
# Initialization
#
#################################
 
# Load libraries
library(ape)
library(Biostrings)
library(dendextend)
library(phyloch)
library(phangorn)
 
# Set paths
exeroot<-"C:\\Users\\Jarno Tuimala\\Desktop\\progs\\"
datapath<-"C:\\Users\\Jarno Tuimala\\Desktop\\"
 
 
#################################
#
# Download sequences
#
#################################
 
## Get sequences
# Download sequences from Genbank
acc<-paste0("X", 67626:67638)
gb<-read.GenBank(acc, species.names=T)
names(gb)<-attr(gb, "species")
#names(gb)<-gsub("_","",names(gb))
 
# Write sequences to disk
write.dna(gb, paste0(datapath, "ratite.fas"), format="fasta")
 
 
#################################
#
# Align sequences
#
#################################
 
## Align sequences
# Alignment
muscleexe<-dir(paste0(exeroot, "muscle"))
musclepath<-paste0(exeroot, "muscle\\", muscleexe)
 
datafile<-dir(datapath, pattern=".fas")
datafilepath<-paste0(datapath, datafile)
 
aligncall<-paste0('"', musclepath, '"', ' -in ', '"', datafilepath, '"', ' -out ', '"', paste0(datapath, gsub(".fas", ".fasta", datafile)), '"')
system(aligncall)      # system, system2, shell
 
# Examine alignment
seqs<-read.dna(file=paste0(datapath, gsub(".fas", ".fasta", datafile)), format="fasta")
#as.alignment(seqs)
#image(seqs)
aln = readDNAMultipleAlignment(filepath=paste0(datapath, gsub(".fas", ".fasta", datafile)), format="fasta")
#alphabetFrequency(aln)
detail(aln)
 
 
#################################
#
# Distance method
#
#################################
 
## Distance analysis
d<-dist.dna(seqs, "logdet")
dist.tree<-nj(d)
dist.tree<-root(dist.tree, "Eudromia_elegans")
 
 
#################################
#
# Model selection
#
#################################
 
mt<-modelTest(phyDat(seqs), dist.tree)
mt$deltaAIC<-mt$AIC-min(mt$AIC)
mt<-mt[order(mt$deltaAIC),]
 
#     Model df    logLik      AIC      BIC   deltaAIC
#23   GTR+G 32 -1449.390 2962.780 3090.904   0.000000
#24 GTR+G+I 33 -1449.361 2964.723 3096.851   1.942628
#22   GTR+I 32 -1451.277 2966.554 3094.678   3.773623
#15   HKY+G 28 -1460.479 2976.959 3089.068  14.178898
#16 HKY+G+I 29 -1460.350 2978.700 3094.813  15.920142
#19   SYM+G 29 -1460.722 2979.443 3095.556  16.663162
#14   HKY+I 28 -1462.568 2981.136 3093.245  18.355926
#20 SYM+G+I 30 -1460.720 2981.440 3101.557  18.660098
#18   SYM+I 29 -1462.403 2982.807 3098.919  20.026670
#11   K80+G 25 -1482.762 3015.525 3115.622  52.744703
#12 K80+G+I 26 -1482.740 3017.479 3121.580  54.699310
#10   K80+I 25 -1483.928 3017.856 3117.953  55.075774
#21     GTR 31 -1521.333 3104.666 3228.786 141.885582
#17     SYM 28 -1528.307 3112.615 3224.724 149.835028
#13     HKY 27 -1540.132 3134.264 3242.369 171.484523
#7    F81+G 27 -1550.204 3154.408 3262.513 191.628428
#9      K80 24 -1553.468 3154.935 3251.028 192.155118
#6    F81+I 27 -1550.795 3155.590 3263.695 192.810514
#8  F81+G+I 28 -1550.204 3156.409 3268.517 193.628646
#3     JC+G 24 -1561.871 3171.742 3267.835 208.962269
#2     JC+I 24 -1562.335 3172.670 3268.763 209.890002
#4   JC+G+I 25 -1561.871 3173.742 3273.839 210.962367
#5      F81 26 -1615.387 3282.774 3386.875 319.994275
#1       JC 23 -1624.635 3295.270 3387.359 332.489701
 
 
 
#################################
#
# Maximum likelihood method
#
#################################
 
## Maximum likelihood analysis
raxmlexe<-dir(paste0(exeroot, "raxml"))
raxmlpath<-paste0(exeroot, "raxml\\", raxmlexe)
 
datafile<-dir(datapath, pattern=".fasta")
datafilepath<-paste0(datapath, datafile)
 
rnd<-sample(1:100, 1)
setwd(tempdir())
raxmlcall<-paste0('"', raxmlpath, '"', ' -s ', '"', datafilepath, '"', ' -n demo ', ' -m GTRGAMMAI ', ' -o Eudromia_elegans ', ' -p ', rnd)
system(raxmlcall)      # system, system2, shell
raxml.tree<-read.tree(dir(tempdir(), full.names=T)[grep("bestTree", dir(tempdir()))])
 
 
#################################
#
# Maximum parsimony method
#
#################################
 
## Maximum parsimony
write.nexus.data(as.matrix(aln), file=paste0(datapath, "ratite.nex"))
write.dna(as.matrix(aln), file=paste0(datapath, "ratite.phy"), format="interleaved", nbcol=-1, colsep="")
phy<-readLines(paste0(datapath, "ratite.phy"))
 
write("xread", file=paste0(datapath, "ratite.tnt"))
write(paste(rev(strsplit(phy[1], split=" ")[[1]]), collapse=" "), file=paste0(datapath, "ratite.tnt"), append=T)
write(phy[-1], file=paste0(datapath, "ratite.tnt"), append=T)
write(";", file=paste0(datapath, "ratite.tnt"), append=T)
 
tntexe<-dir(paste0(exeroot, "tnt"), pattern=".exe")
tntpath<-paste0(exeroot, "tnt\\", tntexe)
 
datafile<-dir(datapath, pattern=".tnt")
datafilepath<-paste0(datapath, datafile)
 
file.copy(datafilepath, paste0(exeroot, "tnt\\", datafile), overwrite=TRUE)
 
write("mxram 200 ;", file=paste0(exeroot, "tnt\\", "script"))
write("nstates DNA ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("nstates NOGAPS ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("procedure ratite.tnt ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("taxname =;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("log ratitelog.txt ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("hold 1000 ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("mult=replic 10;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("bbreak=tbr ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("tsave *tbr-trees.out; save; tsave /;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("export - tbr-trees.tre;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("nelsen * ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("export * ratiteout.nex ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("resample replications 200 ;", file=paste0(exeroot, "tnt\\", "script"), append=T)
write("quit", file=paste0(exeroot, "tnt\\", "script"), append=T)
 
setwd(paste0(exeroot, "tnt"))
tntcall<-c("tnt procedure script")
system(tntcall)
 
pars.tree<-readLines(paste0(exeroot, "tnt\\tbr-trees.tre"))
pars.tree<-pars.tree[grep("\\(", pars.tree)]
pars.tree<-read.tree(text=pars.tree)
 
 
#################################
#
# Bayesian method
#
#################################
 
# Bayesian analysis
bayesexe<-dir(paste0(exeroot, "mrbayes"), pattern=".exe")
bayespath<-paste0(exeroot, "mrbayes\\", bayesexe)
 
write.nexus.data(as.matrix(aln), file=paste0(exeroot, "mrbayes\\ratite.nex"))
tmp<-readLines(paste0(exeroot, "mrbayes\\ratite.nex"))
tmp<-gsub("DATATYPEDNA", "DATATYPE=DNA", tmp)
write(tmp, file=paste0(exeroot, "mrbayes\\ratite.nex"))
 
write("begin mrbayes;", file=paste0(exeroot, "mrbayes\\script.txt"))
write("log start filename=log.txt;", file=paste0(exeroot, "mrbayes\\script.txt"))
write("exe ratite.nex;", file=paste0(exeroot, "mrbayes\\script.txt"), append=TRUE)
write("lset nst=6 rates=invgamma;", file=paste0(exeroot, "mrbayes\\script.txt"), append=TRUE)
write("mcmc ngen=20000 samplefreq=100 printfreq=100 diagnfreq=1000;", file=paste0(exeroot, "mrbayes\\script.txt"), append=TRUE)
write("sump;", file=paste0(exeroot, "mrbayes\\script.txt"), append=TRUE)
write("sumt burnin=5000;", file=paste0(exeroot, "mrbayes\\script.txt"), append=TRUE)
write("end;", file=paste0(exeroot, "mrbayes\\script.txt"), append=TRUE)
 
setwd(paste0(exeroot, "mrbayes\\))
shell("mrbayes_x64.exe script.txt")
bayes.tree<-read.nexus("ratite.nex.trprobs")[[1]]
 
 
#################################
#
# Bootstrapping
#
#################################
 
f <- function(x) nj(dist.dna(x, "logdet"))
tr <- boot.phylo(dist.tree, seqs, f, 1000, trees=TRUE)
plot(dist.tree)
nodelabels(round(tr$BP / 1000, 2), bg="white")
 
 
#################################
#
# Tree comparsions
#
#################################
 
# Topological distances
m<-matrix(ncol=4, nrow=4, data=NA)
diag(m)<-0
rownames(m)<-c("pars.tree","dist.tree","raxml.tree","bayes.tree")
colnames(m)<-c("pars.tree","dist.tree","raxml.tree","bayes.tree")
m[2,1]<-dist.topo(pars.tree, dist.tree)
m[3,1]<-dist.topo(pars.tree, raxml.tree)
m[4,1]<-dist.topo(pars.tree, bayes.tree)
m[3,2]<-dist.topo(dist.tree, raxml.tree)
m[4,2]<-dist.topo(dist.tree, bayes.tree)
m[4,3]<-dist.topo(raxml.tree, bayes.tree)
 
#           pars.tree dist.tree raxml.tree bayes.tree
#pars.tree          0        NA         NA         NA
#dist.tree          2         0         NA         NA
#raxml.tree         2         0          0         NA
#bayes.tree         2         0          0          0
 
 
# Shimodaira-Hasegawa test
SH.test(
optim.pml(pml(dist.tree, phyDat(seqs)), optGamma=TRUE, optInv=TRUE, model="GTR", optBf=TRUE, optNni=FALSE),
optim.pml(pml(acctran(optim.parsimony(pars.tree, phyDat(seqs)), phyDat(seqs)), phyDat(seqs)), optGamma=TRUE, optInv=TRUE, model="GTR", optBf=TRUE, optNni=FALSE),
optim.pml(pml(raxml.tree, phyDat(seqs)), optGamma=TRUE, optInv=TRUE, model="GTR", optBf=TRUE, optNni=FALSE),
optim.pml(pml(acctran(optim.parsimony(bayes.tree, phyDat(seqs)), phyDat(seqs)), phyDat(seqs)), optGamma=TRUE, optInv=TRUE, model="GTR", optBf=TRUE, optNni=FALSE)
)
 
#     Trees      ln L    Diff ln L p-value
#[1,]     1 -1451.267 0.000000e+00  0.7743
#[2,]     2 -1451.696 4.289465e-01  0.3690
#[3,]     3 -1451.267 1.312255e-07  0.7122
#[4,]     4 -1451.696 4.289486e-01  0.3684
 
 
#################################
#
# Visualizations
#
#################################
 
# Visualization
X11()
par(mfrow=c(2,2))
plot(ladderize(dist.tree), main="NJ + logdet")
plot(ladderize(raxml.tree), main="ML + GTRGAMMAI")
plot(ladderize(pars.tree), main="MP")
plot(ladderize(bayes.tree), main="Bayes")
 
X11()
plot(consensus(dist.tree, raxml.tree, pars.tree), main="Strict consensus")
 
X11()
cophyloplot(ladderize(raxml.tree), ladderize(pars.tree), as.matrix(data.frame(raxml.tree$tip.label, raxml.tree$tip.label)), space=180, cex=0.75, length.line=20, lty=1, col="grey75")
 
#X11()
#cophyloplot(ladderize(raxml.tree), ladderize(dist.tree), as.matrix(data.frame(raxml.tree$tip.label, raxml.tree$tip.label)), space=180, cex=0.75, length.line=20, lty=1, col="grey75")

Lopputuloksena on esimerkiksi seuraavanlainen kuva, jossa on verrattu RAxML:n ja TNT:n tuottamia maximum likelihood- ja parsimoniapuita toisiinsa. Harmaat viivat yhdistävät lajit eri puissa, jolloin topologioiden erot on helpompi hahmottaa.

EDIT 2015-04-30: Lisätty uusi skripti phylo6 ja example_output.

EDIT 2015-05-02: Mallinvalintaan soveltuva online-palvelu löytyy osoitteesta https://jtuimala.shinyapps.io/modeltestR2/. Se ottaa syötteenä FastA-muotoisen sekvenssirinnastuksen ja palauttaa taulukon, jossa on testattujen mallien listaus mm. AIC-arvoineen. Parhaiten sopiva malli on taulukossa ylimpänä. Testaus teshdään teknisesti logdet-etäisyyksien pohjalta luodulla neigbor joining -puulla.


Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *