Post date: Mar 27, 2017 5:28:45 PM
As a complement to the LDA approach, I tried detecting recent migration for L. melissa with BayAss. This fits a Bayesian population model similar in some respects to structure to infer contemporary migration rates (migrants within the last two generations) (see the manual and paper). The program is not meant to deal with a large number of SNPs, and in fact I had to alter the source code to increase it to 500. Here is what I did (from the /uufs/chpc.utah.edu/common/home/u6000989/projects/lmelissa_hostAdaptation/demog/ sub-directory)
1. Grab a random subset of 500 SNPs and round to convert to integers
R CMD BATCH subGens.R
2. Convert to BayAss input format (drop ABM)
perl makeBayAss.pl
This makes in_melissa_bayAss.txt
3. Run BayAss. Note that the command line options are not shown except in the source code. Here is my take on them:
-s int seed
-i int number of iterations; def = 1,000,000
-n int thinning; def = 100
-b int burnin; def = 100,000
-o string outfile
-m float deltaM; def = .1
-a float deltaA; def = .1
-f fload deltaF; def = .1
-v logic increases verbose; include
-u logic increases settings; do not include, makes interactive
-g logic increases genotypes; prints out genotypes and migrant ancestries for individuals
-t logic increases trace; include trace file
-p logic increases nolikelihood;
I ran three chains with different seeds, 500,000 iterations a 100,000 iteration burnin and a thinning interval of 200:
BA3 -s 10 -i 500000 -n 200 -b 100000 -o out_melissa_bayAss0.txt -v -g -t in_melissa_bayAss.txt
BA3 -s 137 -i 500000 -n 200 -b 100000 -o out_melissa_bayAss1.txt -v -g -t in_melissa_bayAss.txt
BA3 -s 742 -i 500000 -n 200 -b 100000 -o out_melissa_bayAss2.txt -v -g -t in_melissa_bayAss.txt
4. The output files are messy text files. I thus extracted the point estimates of the migration matrix (25 x 25 populations in alphabetical order). According the the manual the CIs should be given, but you just get the posterir variance (or maybe s.d.?, not clear) instead.
grep "m\[" out_melissa_bayAss0.txt | perl -p -i -e 's/m\[\d+\]\[\d+\]://g' | perl -p -i -e 's/\([0-9\.]+\)//g' > migMat0.txt
grep "m\[" out_melissa_bayAss1.txt | perl -p -i -e 's/m\[\d+\]\[\d+\]://g' | perl -p -i -e 's/\([0-9\.]+\)//g' > migMat1.txt
grep "m\[" out_melissa_bayAss2.txt | perl -p -i -e 's/m\[\d+\]\[\d+\]://g' | perl -p -i -e 's/\([0-9\.]+\)//g' > migMat2.txt
Here is one of the matrixes:
0.7368 0.0084 0.0146 0.0332 0.0188 0.0072 0.0189 0.0095 0.0087 0.0070 0.0082 0.0081 0.0081 0.0106 0.0091 0.0092 0.0075 0.0074 0.0091 0.0064 0.0061 0.0064 0.0238 0.0088 0.0081
0.0078 0.8012 0.0089 0.0067 0.0072 0.0063 0.0071 0.0076 0.0080 0.0100 0.0074 0.0074 0.0076 0.0059 0.0065 0.0069 0.0062 0.0078 0.0080 0.0284 0.0079 0.0085 0.0061 0.0075 0.0072
0.0119 0.0072 0.7158 0.0815 0.0090 0.0056 0.0085 0.0076 0.0076 0.0287 0.0082 0.0068 0.0091 0.0094 0.0065 0.0066 0.0086 0.0086 0.0065 0.0068 0.0082 0.0080 0.0070 0.0080 0.0086
0.0067 0.0055 0.0067 0.7907 0.0072 0.0063 0.0068 0.0055 0.0077 0.0444 0.0069 0.0097 0.0074 0.0075 0.0071 0.0068 0.0077 0.0075 0.0057 0.0057 0.0071 0.0069 0.0062 0.0084 0.0115
0.0086 0.0090 0.0131 0.0635 0.6975 0.0176 0.0101 0.0097 0.0104 0.0089 0.0091 0.0108 0.0101 0.0094 0.0108 0.0100 0.0094 0.0097 0.0095 0.0093 0.0099 0.0088 0.0098 0.0091 0.0160
0.0087 0.0067 0.0138 0.0776 0.0066 0.7242 0.0066 0.0071 0.0079 0.0124 0.0078 0.0112 0.0068 0.0082 0.0071 0.0075 0.0072 0.0077 0.0074 0.0082 0.0095 0.0062 0.0072 0.0075 0.0190
0.0100 0.0080 0.0098 0.0294 0.0075 0.0062 0.7340 0.0071 0.0064 0.0100 0.0090 0.0086 0.0076 0.0153 0.0078 0.0306 0.0079 0.0068 0.0067 0.0096 0.0087 0.0057 0.0252 0.0078 0.0143
0.0084 0.0071 0.0087 0.0353 0.0104 0.0084 0.0087 0.7796 0.0074 0.0074 0.0080 0.0083 0.0081 0.0073 0.0069 0.0106 0.0075 0.0141 0.0069 0.0083 0.0055 0.0072 0.0076 0.0064 0.0060
0.0074 0.0086 0.0084 0.0163 0.0084 0.0080 0.0087 0.0079 0.7900 0.0072 0.0077 0.0076 0.0085 0.0078 0.0090 0.0062 0.0075 0.0080 0.0092 0.0068 0.0073 0.0067 0.0074 0.0236 0.0056
0.0077 0.0058 0.0261 0.0872 0.0057 0.0076 0.0059 0.0056 0.0076 0.7301 0.0079 0.0074 0.0061 0.0068 0.0058 0.0081 0.0056 0.0064 0.0078 0.0078 0.0059 0.0075 0.0065 0.0080 0.0134
0.0076 0.0128 0.0069 0.0161 0.0102 0.0077 0.0334 0.0056 0.0082 0.0062 0.7274 0.0117 0.0080 0.0169 0.0076 0.0346 0.0065 0.0072 0.0073 0.0074 0.0077 0.0079 0.0193 0.0077 0.0077
0.0075 0.0069 0.0117 0.0503 0.0056 0.0186 0.0073 0.0075 0.0090 0.0095 0.0109 0.7572 0.0071 0.0062 0.0085 0.0070 0.0086 0.0086 0.0073 0.0062 0.0063 0.0095 0.0081 0.0076 0.0071
0.0071 0.0071 0.0074 0.0070 0.0066 0.0057 0.0075 0.0079 0.0091 0.0073 0.0082 0.0078 0.8167 0.0085 0.0084 0.0094 0.0070 0.0069 0.0081 0.0084 0.0076 0.0066 0.0080 0.0081 0.0076
0.0092 0.0077 0.0129 0.0482 0.0093 0.0079 0.0079 0.0070 0.0073 0.0105 0.0133 0.0057 0.0058 0.7404 0.0076 0.0147 0.0076 0.0076 0.0086 0.0073 0.0073 0.0077 0.0221 0.0069 0.0094
0.0066 0.0080 0.0087 0.0072 0.0079 0.0058 0.0072 0.0077 0.0068 0.0062 0.0066 0.0070 0.0077 0.0088 0.8033 0.0071 0.0086 0.0089 0.0080 0.0057 0.0083 0.0086 0.0096 0.0221 0.0075
0.0137 0.0071 0.0127 0.0306 0.0173 0.0149 0.0112 0.0093 0.0086 0.0101 0.0166 0.0077 0.0098 0.0163 0.0082 0.7332 0.0077 0.0079 0.0083 0.0081 0.0085 0.0079 0.0089 0.0076 0.0080
0.0073 0.0083 0.0096 0.0138 0.0071 0.0073 0.0075 0.0093 0.0080 0.0073 0.0088 0.0080 0.0059 0.0085 0.0166 0.0076 0.7807 0.0080 0.0183 0.0073 0.0078 0.0077 0.0071 0.0160 0.0064
0.0094 0.0068 0.0057 0.0138 0.0076 0.0067 0.0064 0.0288 0.0076 0.0073 0.0064 0.0073 0.0091 0.0078 0.0082 0.0088 0.0071 0.7670 0.0075 0.0078 0.0310 0.0075 0.0089 0.0071 0.0085
0.0065 0.0072 0.0068 0.0074 0.0072 0.0071 0.0070 0.0065 0.0154 0.0087 0.0068 0.0066 0.0078 0.0069 0.0062 0.0073 0.0093 0.0082 0.8077 0.0061 0.0074 0.0077 0.0071 0.0166 0.0086
0.0096 0.0107 0.0097 0.0086 0.0066 0.0101 0.0090 0.0085 0.0084 0.0094 0.0089 0.0092 0.0082 0.0086 0.0089 0.0092 0.0097 0.0081 0.0085 0.7864 0.0096 0.0096 0.0090 0.0079 0.0073
0.0087 0.0069 0.0053 0.0071 0.0082 0.0071 0.0071 0.0483 0.0074 0.0072 0.0077 0.0089 0.0075 0.0090 0.0085 0.0069 0.0078 0.0080 0.0084 0.0069 0.7752 0.0079 0.0075 0.0083 0.0080
0.0084 0.0089 0.0063 0.0073 0.0075 0.0096 0.0075 0.0063 0.0086 0.0075 0.0073 0.0075 0.0071 0.0073 0.0099 0.0053 0.0071 0.0067 0.0069 0.0077 0.0092 0.8095 0.0080 0.0160 0.0065
0.0080 0.0084 0.0065 0.0772 0.0179 0.0184 0.0064 0.0080 0.0055 0.0092 0.0075 0.0191 0.0068 0.0088 0.0077 0.0066 0.0063 0.0072 0.0079 0.0059 0.0082 0.0068 0.7200 0.0081 0.0074
0.0081 0.0060 0.0091 0.0065 0.0077 0.0078 0.0070 0.0071 0.0198 0.0062 0.0064 0.0076 0.0082 0.0087 0.0055 0.0080 0.0060 0.0079 0.0075 0.0071 0.0068 0.0068 0.0074 0.8139 0.0069
0.0090 0.0153 0.0109 0.0987 0.0083 0.0071 0.0066 0.0073 0.0066 0.0265 0.0066 0.0082 0.0072 0.0077 0.0080 0.0069 0.0059 0.0079 0.0086 0.0068 0.0078 0.0070 0.0077 0.0071 0.7002
5. I summarized this in R. A few points. Most off-diagonal migration rates are < 0.01 (that is those between populations). This ranges from 84-85% across chains. Indeed, few high migration rates are seen, here is the 95th quantile and max for each chain:
95% 100%
0.018905 0.081500
0.021755 0.050300
0.020620 0.04260
This all makes good sense. However, the proportion non-migrant (the diagonal) is lower than is probably sensible (summaries across chains):
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.6975 0.7301 0.7670 0.7615 0.7907 0.8167
0.6969 0.7306 0.7625 0.7630 0.7924 0.8156
0.7019 0.7414 0.7615 0.7649 0.7921 0.8154
I think this is just because there are so many populations and thus even with low migration for each population the sum of migrants across populations is unrealistically high. And we just can confidently exclude all migrants for most populations (thus this adds up). In other words, I don't think this is really working in this context.