Post date: Jun 25, 2018 5:57:23 PM
I re-estimate Ne, now based on the 2013-2017 data (and played some with subsets).
First I made infiles for varne using mkNeInfiles.pl
#!/usr/bin/perl
#
open(IN, "SampleSizes.txt") or die;
while(<IN>){
chomp;
@line = split(" ",$_);
$ps{$line[0]} = $line[1];
}
close(IN);
foreach $p (sort keys %ps){
$in = "../AlleleFreqs/p_"."$p".".txt";
$out = "in_"."$p".".txt";
open(IN, "cut -f 3 -d \" \" $in | ") or die "failed to read pipe\n";
open (OUT, "> $out") or die "failed to write\n";
while(<IN>){
chomp;
$n = 2 * $ps{$p};
print OUT "$_ $n\n";
}
close(IN);
close(OUT);
}
Then, I terated over pairs of popualtions with runNeEst.pl (takes the 2013 files as infiles).
#!/usr/bin/perl
#
foreach $in13 (@ARGV){
$in13 =~ m/in_([A-Z]+)/;
$pop = $1;
$in17 = $in13;
$in17 =~ s/13/17/ or die "failed sub $in13\n";
system "varne -a $in13 -b $in17 -l 30513 -t 4 -n 3000 -x 1000 > Ne_$pop"."_13-17.txt\n";
}
Point estimates (posterior means) are in SummaryNe13-17.txt.
I then made some summaries in R (commands.R).
## N from distance, number of years varies, from 2013-2016
bcr<-c(2382,1241.7)
bnp<-c(633.9,1273.2)
btb<-c(1838.7,1978.5,2763.3)
gnp<-c(1119.9,1024.5)
hnv<-c(5291.4)
mrf<-c(977.7)
psp<-366.6
rnv<-NA
ski<-c(1348.8,1242.2)
usl<-c(1708.2,2927.1)
ne<-scan("ne")
N<-c(mean(bcr),mean(bnp),mean(btb),mean(gnp),mean(hnv),mean(mrf),mean(psp),mean(rnv),mean(ski),mean(usl))
L<-30513
pfiles<-c("in_BCR-13.txt","in_BNP-13.txt","in_BTB-13.txt","in_GNP-13.txt","in_HNV-13.txt","in_MRF-13.txt","in_PSP-13.txt","in_RNV-13.txt","in_SKI-13.txt","in_USL-13.txt")
P<-matrix(NA,nrow=L,ncol=10)
for(i in 1:10){
x<-read.table(pfiles[i])
P[,i]<-x[,1]
}
het<-2*P * (1-P)
H<-apply(het,2,mean)