%%bash
## create directories for data assembly
mkdir -p analysis_raxml
mkdir -p analysis_bucky/locs ## Bucky analysis was not included in manuscript
mkdir -p figures/
Phylogenetic trees are inferred by maximum likelihood using raxml with bootstrap support inferred from 200 random starting trees. Input files are in phylip format containing concatenated RADseq loci output by PyRAD.
raxcmd = "raxmlHPC-PTHREADS-AVX -f a" +\
" -m GTRGAMMA"+\
" -N 200 -x 12345 -p 12345 -T 20 -o NI,HE "+\
" -w /home/deren/Documents/Oaks/Virentes/analysis_raxml/ "
phy = "-s analysis_pyrad/outfiles/virentes_c85d6m4p5.phy -n virentes_m4p5"
stderr = ! $raxcmd $phy
phy = "-s analysis_pyrad/outfiles/virentes_c85d6m20p5.phy -n virentes_m20p5"
stderr = ! $raxcmd $phy
phy = "-s analysis_pyrad/outfiles/MGV_c85d6m4p5.phy -n MGV_m4p5"
stderr = ! $raxcmd $phy
phy = "-s analysis_pyrad/outfiles/MGV_c85d6m16p5.phy -n MGV_m16p5"
stderr = ! $raxcmd $phy
phy = "-s analysis_pyrad/outfiles/OS_c85d6m4p5.phy -n OS_m4p5"
stderr = ! $raxcmd $phy
phy = "-s analysis_pyrad/outfiles/OS_c85d6m13p5.phy -n OS_m13p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/FB_c85d6m4p5.phy -n FB_m4p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/FB_c85d6m12p5.phy -n FB_m12p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/OSMGV_c85d6m4p5.phy -n OSMGV_m4p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/OSMGV_c85d6m20p5.phy -n OSMGV_m20p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/FBMGV_c85d6m4p5.phy -n FBMGV_m4p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/FBMGV_c85d6m20p5.phy -n FBMGV_m20p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/FBOS_c85d6m4p5.phy -n FBOS_m4p5"
stderr = ! $raxcmd $phy
phy = " -s analysis_pyrad/outfiles/FBOS_c85d6m16p5.phy -n FBOS_m16p5"
stderr = ! $raxcmd $phy
import itertools
import ete2
import numpy as np
import pandas as pd
import glob
from collections import OrderedDict, Counter
import matplotlib.pyplot as plt
%matplotlib inline
## read in the two full-taxa data sets
infile4 = "analysis_pyrad/outfiles/virentes_c85d6m4p5.loci"
infile20 = "analysis_pyrad/outfiles/virentes_c85d6m20p5.loci"
loci4 = open(infile4).read().split("|\n")
loci20 = open(infile20).read().split("|\n")
## how many loci
nloci4 = len(loci4)
nloci20 = len(loci20)
## read in a tree file for the order of tips
tree = ete2.Tree("analysis_raxml/RAxML_bipartitions.virentes_m4p5", format=2)
## sort names by their order in the tree
tree.ladderize()
names = [i.strip() for i in tree.iter_leaf_names()]
## create list of locus x sample
def create_array(nsamp, nloci, loci):
lxs = np.zeros((nsamp, nloci))
for loc in xrange(len(loci)):
for seq in loci[loc].split("\n"):
if ">" in seq:
lxs[names.index(seq.split()[0].replace(">","")),loc] += 1
return lxs
## create arrays for the two data sets
lxs1 = create_array(len(names), nloci4, loci4)
lxs2 = create_array(len(names), nloci20, loci20)
## empty matrix
share = np.zeros((len(names),len(names)))
For each column (locus) in lxs iterate over all pairs of samples and if both have data at that locus add 1 to their shared data in the matrix [share]. This fills the upper triangle of the matrix. I fill the lower triangle with shared data from the other data set.
## fill matrix share above diagonal (Allmin20 data set)
for col in xrange(lxs2.shape[1]):
for samp1,samp2 in itertools.combinations(names,2):
if (lxs2[names.index(samp1),col] and lxs2[names.index(samp2),col]):
share[names.index(samp1),names.index(samp2)] += 1
## fill matrix share below diagonal (Allmin4 data set)
for col in xrange(lxs1.shape[1]):
for samp2,samp1 in itertools.combinations(names,2):
if (lxs1[names.index(samp1),col] and lxs1[names.index(samp2),col]):
share[names.index(samp1),names.index(samp2)] += 1
f,ax = plt.subplots(1,1, figsize=(15,10))
cax = ax.imshow(share, interpolation='none', cmap='gist_yarg')
ax.grid(True, color='white', lw=0.9)
ax.get_xaxis().set_ticks([i-0.5 for i in range(35)])
ax.get_xaxis().set_ticklabels([])
ax.get_yaxis().set_ticks([i-0.5 for i in range(35)])
ax.get_yaxis().set_ticklabels([])
cb = plt.colorbar(cax);
for t in cb.ax.get_yticklabels():
t.set_fontsize(20)
f.savefig("figures/heatmap.svg")
f.savefig("figures/heatmap.png", dpi=300)
## get number of raw reads for each sample
with open("analysis_pyrad/stats/s2.rawedit.txt") as rawstat:
raws = rawstat.readlines()
dat = [(line.split("\t")[0],int(line.split("\t")[-1].strip())) for \
line in raws if "_" in line]
dat2 = OrderedDict()
for t in names:
for j in dat:
if t == j[0].split("_")[0]:
if t not in dat2:
dat2[t] = j[1]
else:
dat2[t] += j[1]
dat = pd.DataFrame({"raws":pd.Series(dat2)})
## make plot of raw reads
f, ax = plt.subplots(1,1, figsize=(10,2))
ax.bar(range(len(dat)),(dat.raws), color="#262626")
ax.grid('off')
f.savefig("figures/bars.svg")
from skbio import DistanceMatrix
from skbio.stats.distance import mantel
## func to measure jaccard's distance
def jaccard(a, b):
c = a.intersection(b)
return float(len(c)) / (len(a) + len(b) - len(c))
## write matrices to file for loading into R for mantel test
Pairs4 = OrderedDict()
for samp1,samp2 in itertools.product(names,repeat=2):
if samp1 == samp2:
Pairs4[str(samp1+"_"+samp2)] = [0.,0.,0.]
else:
jaccar = jaccard(set([i for i in range(lxs1.shape[1]) if \
lxs1[:,i][names.index(samp1)] == 1.]),
set([i for i in range(lxs1.shape[1]) if \
lxs1[:,i][names.index(samp2)] == 1.]))
meanraw = np.mean([dat.raws[samp1],dat.raws[samp2]])
phydist = tree.get_distance(samp1,samp2)
Pairs4[str(samp1+"_"+samp2)] = [jaccar,meanraw,phydist]
sharedist = DistanceMatrix(np.array([i[0] for i in Pairs4.values()]).reshape(34,34))
rawdist = DistanceMatrix(np.array([i[1] for i in Pairs4.values()]).reshape(34,34))
phydist = DistanceMatrix(np.array([round(i[2],7) for i in Pairs4.values()]).reshape(34,34))
coeff, p_value, n = mantel(sharedist, rawdist,
method='spearman',
permutations=9999,
alternative="two-sided")
print "shared loci ~ mean raw data"
print "coeff =",coeff
print "p-value =",p_value
shared loci ~ mean raw data coeff = 0.372246206797 p-value = 0.0106
coeff, p_value, n = mantel(sharedist, phydist,
method='spearman',
permutations=9999,
alternative="two-sided")
print "shared loci ~ phylo dist"
print "coeff =",coeff
print "p-value =", p_value
shared loci ~ phylo dist coeff = -0.144901236203 p-value = 0.2401
f,ax = plt.subplots(1,2, figsize=(8,4), sharey=True)
## remove values on the diagonal
r = rawdist.data[rawdist.data.nonzero()]
s = sharedist.data[sharedist.data.nonzero()]
p = phydist.data[phydist.data.nonzero()]
ax[0].scatter(r/1e6,s, color="#262626", alpha=0.3, s=50)
ax[1].scatter(p,s, color="#262626", alpha=0.3, s=50)
ax[0].set_ylabel("Jaccard's distance of shared loci", fontsize=14)
ax[0].set_xlabel("Pair-wise mean raw reads (x10$^6$)", fontsize=14)
ax[1].set_xlabel("Pair-wise phylogenetic distance", fontsize=14)
f.savefig("figures/rawVphylo_min4.png")
#f.savefig("figures/rawVphylo_min4.svg")
## write matrices to file for loading into R for mantel test
Pairs20 = OrderedDict()
for samp1,samp2 in itertools.product(names,repeat=2):
if samp1 == samp2:
Pairs20[str(samp1+"_"+samp2)] = [0.,0.,0.]
else:
jaccar = jaccard(set([i for i in range(lxs2.shape[1]) if \
lxs2[:,i][names.index(samp1)] == 1.]),
set([i for i in range(lxs2.shape[1]) if \
lxs2[:,i][names.index(samp2)] == 1.]))
meanraw = np.mean([dat.raws[samp1],dat.raws[samp2]])
phydist = tree.get_distance(samp1,samp2)
Pairs20[str(samp1+"_"+samp2)] = [jaccar,meanraw,phydist]
sharedist20 = DistanceMatrix(np.array([i[0] for i in Pairs20.values()]).reshape(34,34))
rawdist20 = DistanceMatrix(np.array([i[1] for i in Pairs20.values()]).reshape(34,34))
phydist20 = DistanceMatrix(np.array([round(i[2],7) for i in Pairs20.values()]).reshape(34,34))
coeff, p_value, n = mantel(sharedist, rawdist,
method='spearman',
permutations=9999,
alternative="two-sided")
print "shared loci ~ mean raw data"
print "coeff =",coeff
print "p-value =",p_value
shared loci ~ mean raw data coeff = 0.478713410308 p-value = 0.0028
coeff, p_value, n = mantel(sharedist, phydist,
method='spearman',
permutations=9999,
alternative="two-sided")
print "shared loci ~ phylo data"
print "coeff =",coeff
print "p-value =",p_value
shared loci ~ phylo data coeff = 0.0869666708908 p-value = 0.5273
f,ax = plt.subplots(1,2, figsize=(8,4), sharey=True)
## remove values on the diagonal
r = rawdist20.data[rawdist.data.nonzero()]
s = sharedist20.data[sharedist.data.nonzero()]
p = phydist20.data[phydist.data.nonzero()]
ax[0].scatter(r/1e6,s, color="#262626", alpha=0.3, s=50)
ax[1].scatter(p,s, color="#262626", alpha=0.3, s=50)
ax[0].set_ylabel("Jaccard's distance of shared loci", fontsize=14)
ax[0].set_xlabel("Pair-wise mean raw reads (x10$^6$)", fontsize=14)
ax[1].set_xlabel("Pair-wise phylogenetic distance", fontsize=14)
f.savefig("figures/rawVphylo_min20.png")
#f.savefig("figures/rawVphylo.svg")
## Dictionary to select one sample for each species
D = {"B":"BJSL25",
"Ft":"TXMD3",
"Fm":"MXED8",
"O":"CRL0030",
"S":"CUVN10",
"V":"FLBA140",
"G":"FLSF54",
"M":"FLSF47",
"out":"AR"}
def nexmake(D, N, outg):
outloc = open("analysis_bucky/locs/"+`N`+".nex", 'w')
## write header block
outloc.write("#NEXUS\nbegin data;\n"+\
"dimensions ntax="+`len(D)`+" nchar="+`len(D.values()[0])`+";\n"+\
"format datatype=dna interleave=yes gap=-;\n"+\
"matrix\n\n")
## write sequence data for locus loc
for i in D.items():
outloc.write( i[0][0:9]+" "*(10-len(i[0][0:9]))+"".join(i[1])+"\n" )
## write mrbayes block
outloc.write("\n\n ;\nend;\n"+\
"\n\nbegin mrbayes;\n"+\
"set autoclose=yes nowarn=yes;\n"+\
"lset nst=2 rates=gamma;\n"+\
"outgroup "+outg+";\n"+\
"mcmc ngen=4000000 samplefreq=4000;\n"+\
"sump burnin=1000000;\n"+\
"sumt burnin=1000000;\n"+\
"end;")
outloc.close()
## In latter versions of pyrad one could use the
## .alleles data file ('a' output from pyrad)
infile = open("/home/deren/Documents/Oaks/Virentes/analysis_pyrad/outfiles/virentes_c85d6m4p5.loci").read()
loci = iter(infile.strip().split("|\n")[:])
def unstruct(amb):
" returns bases from ambiguity code"
D = {"R":["G","A"],
"K":["G","T"],
"S":["G","C"],
"Y":["T","C"],
"W":["T","A"],
"M":["C","A"]}
if amb in D:
return D.get(amb)
else:
return [amb,amb]
def resolveambig(subseq):
N = []
for col in subseq:
N.append([unstruct(i)[np.random.binomial(1,0.5)] for i in col])
return np.array(N)
def PIS(seqsamp):
## check that locus is variable
pis = sum([Counter(col).most_common(2)[1][1] > 1 for col in seqsamp.T \
if len(Counter(col).most_common(2))>1])
if pis >= 2:
return pis
else:
return 0
N = 0
## create subsampled data set
while 1:
try: loc = loci.next().split("\n")[:-1]
except StopIteration: break
## if all tip samples have data in this locus
names = [i.split()[0].replace(">","") for i in loc]
## check that locus has required samples for each subtree
if all([i in names for i in D.values()]):
seqs = np.array([list(i.split()[1]) for i in loc])
seqsamp = seqs[[names.index(tax) for tax in D.values()]]
seqsamp = resolveambig(seqsamp)
pis = PIS(seqsamp)
if pis:
N += 1
## remove invariable columns given this subsampling
keep = []
for n,col in enumerate(seqsamp.T):
if all([i not in ["N","-"] for i in col]):
keep.append(n)
subseq = seqsamp.T[keep].T
## write to a nexus file
DD = dict(zip(D.keys(), [i.tostring() for i in subseq]))
nexmake(DD, N, "out")
print N, 'loci kept'
828 loci kept
## create a Parallel computing client
## can map a function across how ever many
## parallel engines are turned on in the notebook
## under the Clusters tab
from IPython.parallel import Client
PAR = Client()[:]
LOADBALANCE = Client().load_balanced_view()
## list of arguments to function single_mb
args = []
## create a new dir for log files
! mkdir -p "analysis_bucky"/logs/
## get all nex files from this dir
files = glob.glob("analysis_bucky/locs/*.nex")
## put nex file as input and other as output location into args
for ff in files:
num = ff.split("/")[-1]
args.append([ff, "analysis_bucky/logs/log."+str(num)])
print '------\n',len(args), '\ttotal mrbayes analyses'
------ 828 total mrbayes analyses
## decorator adds parallel capability if run in IPython notebook
@LOADBALANCE.parallel(block=True)
def single_mb(arg):
infile, outloc = arg
! mb "$infile" > "$outloc"
## run these jobs in parallel
quiet = single_mb.map(args)
## create dir for bucky input files
! mkdir -p analysis_bucky/ins/
## iterate over dirs
K = glob.glob("analysis_bucky/locs/*.nex")
nums = [int(i.split("/")[-1].split('.nex')[0]) for i in K]
## combine trees from the two replicate runs
for n in range(1,max(nums)+1):
cmd = """ -o analysis_bucky/ins/%d.in \
analysis_bucky/locs/%d.nex.run1.t \
analysis_bucky/locs/%d.nex.run2.t """ % (n,n,n)
stderr = ! mbsum -n 101 $cmd
! bucky -h
Bayesian Untangling of Concordance Knots (applied to yeast and other organisms) BUCKy version 1.4.3, 9 Jul 2014 Copyright (C) 2006-2014 by Bret Larget and Cecile Ane and Riley Larget This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. Usage: bucky [options] <input files> Options: Parameter | Usage | Default Value ------------------------------------------------------------------- alpha | -a number | 1 # of runs | -k integer | 2 # of MCMC updates | -n integer | 100000 # of chains | -c integer | 1 MCMCMC Rate | -r integer | 100 alpha multiplier | -m number | 10 subsample rate | -s integer | 1 output root file name | -o name | run1 input file list file | -i filename | random seed 1 | -s1 integer | 1234 random seed 2 | -s2 integer | 5678 CF cutoff for display | -cf number | 0.05 create sample file | --create-sample-file | false create joint file | --create-joint-file | false create single file | --create-single-file | false use independence prior | --use-independence-prior | false calculate pairs | --calculate-pairs | false taxon set | -p prune-file | common taxa skip genes with fewer taxa | -sg | false use update groups | --use-update-groups | true use update groups | --do-not-use-update-groups | Space optimization | --opt-space | false do not build population tree? | --no-population-tree | false grid-size for genomewide CF | --genomewide-grid-size | 1000 help | -h OR --help | version | --version | -------------------------------------------------------------------
! bucky -a 1 -k 4 -c 4 -n 1000000 -o run1 analysis_bucky/ins/*
! bucky -a 5 -k 4 -c 4 -n 1000000 -o run5 analysis_bucky/ins/*
! bucky -a 20 -k 4 -c 4 -n 1000000 -o run20 analysis_bucky/ins/*
! bucky --use-independence-prior -k 4 -n 1000000 -o runinf analysis_bucky/ins/*
! head -n 50 analysis_bucky/run1.concordance
translate 1 B, 2 Ft, 3 G, 4 M, 5 O, 6 S, 7 V, 8 Fm, 9 out; Population Tree: ((((1,8),2),((3,4),7)),(5,6),9); Primary Concordance Tree Topology: ((1,(2,8)),(((3,4),7),(5,6)),9); Population Tree, With Branch Lengths In Estimated Coalescent Units: ((((1:10.000,8:10.000):0.086,2:10.000):0.308,((3:10.000,4:10.000):0.086,7:10.000):0.245):0.023,(5:10.000,6:10.000):0.395,9:10.000); Primary Concordance Tree with Sample Concordance Factors: ((1:1.000,(2:1.000,8:1.000):0.206):0.219,(((3:1.000,4:1.000):0.199,7:1.000):0.136,(5:1.000,6:1.000):0.371):0.103,9:1.000); Four-way partitions in the Population Tree: sample-wide CF, coalescent units and Ties(if present) {1,2,5,6,8,9; 7|3; 4} 0.388, 0.086, {1; 8|2; 3,4,5,6,7,9} 0.388, 0.086, {1,2,3,4,7,8; 9|5; 6} 0.551, 0.395, {1,2,8; 5,6,9|3,4; 7} 0.478, 0.245, {1,8; 2|3,4,7; 5,6,9} 0.510, 0.308, {1,2,8; 3,4,7|5,6; 9} 0.349, 0.023, Splits in the Primary Concordance Tree: sample-wide and genome-wide mean CF (95% credibility), SD of mean sample-wide CF across runs {1,2,3,4,7,8,9|5,6} 0.371(0.326,0.429) 0.370(0.313,0.439) 0.027 {1,2,8|3,4,5,6,7,9} 0.219(0.174,0.262) 0.219(0.166,0.273) 0.012 {1,3,4,5,6,7,9|2,8} 0.206(0.147,0.280) 0.206(0.141,0.288) 0.038 {1,2,5,6,7,8,9|3,4} 0.199(0.156,0.240) 0.199(0.149,0.249) 0.016 {1,2,5,6,8,9|3,4,7} 0.136(0.105,0.167) 0.136(0.098,0.176) 0.010 {1,2,8,9|3,4,5,6,7} 0.103(0.054,0.164) 0.103(0.049,0.166) 0.026 Splits NOT in the Primary Concordance Tree but with estimated CF > 0.050: {1,8|2,3,4,5,6,7,9} 0.206(0.151,0.273) 0.206(0.145,0.281) 0.032 {1,2,3,5,6,8,9|4,7} 0.153(0.092,0.191) 0.153(0.088,0.201) 0.009 {1,2,4,5,6,8,9|3,7} 0.148(0.099,0.185) 0.148(0.094,0.194) 0.012 {1,2|3,4,5,6,7,8,9} 0.137(0.081,0.186) 0.137(0.077,0.194) 0.029 {1,2,3,4,7,8|5,6,9} 0.102(0.072,0.133) 0.102(0.067,0.141) 0.015 {1,9|2,3,4,5,6,7,8} 0.088(0.047,0.116) 0.088(0.046,0.123) 0.009 {1,2,5,6,8|3,4,7,9} 0.084(0.052,0.129) 0.084(0.047,0.134) 0.009 {1,2,3,4,5,6,8|7,9} 0.083(0.035,0.115) 0.083(0.033,0.123) 0.011 {1,2,3,4,5,7,8|6,9} 0.074(0.057,0.098) 0.074(0.050,0.104) 0.004 {1,2,3,4,6,7,8|5,9} 0.072(0.043,0.109) 0.072(0.040,0.113) 0.014
! head -n 50 analysis_bucky/run5.concordance
! head -n 50 analysis_bucky/run20.concordance
! head -n 50 analysis_bucky/runinf.concordance
%load_ext rmagic
## requires that the Python package Rpy2 is installed
The rmagic extension is already loaded. To reload it, use: %reload_ext rmagic
%%R
library(ape)
alpha1 = read.tree(text="((1,(2,8)),(((3,4),7),(5,6)),9);")
alpha5 = read.tree(text="((((1,8),2),((3,4),7)),(5,6),9);")
alpha20 = read.tree(text="(((1,8),2),(((3,4),7),(5,6)),9);")
alphainf = read.tree(text="(((((1,8),2),(5,6)),(3,4)),7,9);")
%%R -w 400 -h 400
par(mfrow=c(2,2),
mar=c(2,1,1,2))
plot(alpha1, edge.width=4, cex=2, adj=0.1)
plot(alpha5, edge.width=4, cex=2, adj=0.1)
plot(alpha20, edge.width=4, cex=2, adj=0.1)
plot(alphainf, edge.width=4, cex=2, adj=0.1)
%%R -w 600 -h 800
library(ape)
## save figure
#svg("figures/tree_virentes_m4p5.svg",height=8, width=6)
#tre = read.tree("analysis_raxml/RAxML_bipartitions.virentes_m4p5")
#plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
#nodelabels(tre$node.label, bg='grey', cex=0.8)
#dev.off()
## show figure
tre = read.tree("analysis_raxml/RAxML_bipartitions.virentes_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R
## save figure
svg("figures/tree_virentes_m20p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.virentes_m20p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
tre = read.tree("analysis_raxml/RAxML_bipartitions.virentes_m20p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -n
svg("figures/cladogram.svg", height=9, width=3)
plot(ladderize(tre), edge.width=3, cex=0.75, adj=0.1, use.edge.length=F)
dev.off()
png 2
%%R -w 300
svg("figures/ingrouptree.svg", height=9, width=3)
tre <- drop.tip(tre, c("HE","NI","CH","DU","DO","AR","EN"))
tre <- ladderize(tre)
tre <- rotate(tre,28)
tre <- rotate(tre,29)
tre <- rotate(tre,30)
plot(tre, edge.width=4, cex=0.75, adj=0.1)
dev.off()
png 2
%%R -w 500
## save figure
svg("figures/tree_FBMGV_m4p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.FBMGV_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_FBMGV_m20p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.FBMGV_m20p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_OSMGV_m4p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.OSMGV_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_OSMGV_m20p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.OSMGV_m20p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_FBOS_m4p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.FBOS_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_FBOS_m16p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.FBOS_m16p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_MGV_m4p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.MGV_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_MGV_m16p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.MGV_m16p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_OS_m4p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.OS_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_OS_m13p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.OS_m13p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_FB_m4p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.FB_m4p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
%%R -w 500
## save figure
svg("figures/tree_FB_m12p5.svg",height=8, width=6)
tre = read.tree("analysis_raxml/RAxML_bipartitions.FB_m12p5")
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)
dev.off()
## show figure
plot(ladderize(tre), edge.width=4, cex=0.8, adj=0.1)
nodelabels(tre$node.label, bg='grey', cex=0.8)