angsd Abbababa2: negative blockstart/end

[joanocha]

I am currently using Abbababa2 in angsd and I found what seems to be a bug in the output.Angsd file. This is an example of how a good chr looks like when run for 5M blocks (angsd default):

CHR BLOCKstart BLOCKend Numer Denom numSites
chr1 0 4999999 -133.344553 2123.291220 64251.000000
chr1 5000000 9999999 -78.996484 1165.778082 33456.000000
chr1 10000000 14999999 -28.215509 1489.065549 40988.000000

This is how a weird chr looks like (I have a few looking like this):
chr4 -478885888 -473885889 -291.655995 10272.699393 348202.000000

The blocks shouldn't be negative. I also noticed that the chrs that show this pattern change as I change the number of populations/comparisons.

Here is the command I run:
angsd -doAbbababa2 1 -bam list_bams.txt -sizeFile sizefile.txt -doCounts 1 -out output.Angsd -anc refgenome.fasta -sites quality_control_snps.pos -rf quality_control_snps_scaffoldslist.rf
-useLast 0 -minMapQ 30 -minQ 20 -remove_bads 1 -uniqueOnly 1 -only_proper_pairs 1 -p 1

Instead of using a -rf in the format of:
chr1:6-24
chr1:26-71
chr1:73-140

I feed a position file and a scaffold list:
-sites file looks like:
chr1 6
chr1 7
chr1 8
chr1 9
chr1 10
(...)

-rf file is a simple list of scaffolds:
chr1
chr2
(...)

So far it seemed to work well for most of the chr. What do you think could be creating this issue of negative blocks for some of the chromosomes? Please correct me if this is just me not running Angsd command properly.
Thank you!

[SamueleSoraggi]

Hi joanocha,
sorry for such a late answer, your issue got down in the list and I overlooked it!

This is very weird. But I noticed that the difference between 478885888 and 473885889 is exactly the size of 1 block of 5000000 bases. I check what happened in the abbababa2 script and hopefully find the error, if it is due to that script.

However, could you try to write

chr1:
chr2:
...

using also the semicolumns in the -rf file?

[joanocha]

Hi Samuele! Thanks for your reply! Better late than never :) The difference is always the block size. I did another experience running with a 200 kb block for example and I got the same result for chromosome 4. This time I run with a -rf file with the type chr4:posStart-posEnd. Which led me to exclude the -rf as the issue. best Joana chr4 -392185856 -391985857 chr4 -392185856 -391985857 chr4 -392185856 -391985857

…

On Mon, Feb 3, 2020 at 8:11 AM Samuele Soraggi ***@***.***> wrote: Hi joanocha, sorry for such a late answer, your issue got down in the list and I overlooked it! This is very weird. But I noticed that the difference between 478885888 and 473885889 is exactly the size of 1 block of 5000000 bases. I check what happened in the abbababa2 script and hopefully find the error, if it is due to that script. However, could you try to write chr1: chr2: using also the semicolumns in the -rf file? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#273?email_source=notifications&email_token=AKH5573DE7MQZA46ZLCCWELRBA64RA5CNFSM4JNTRNQKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEKUNMSI#issuecomment-581490249>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKH5576KPE6XPNXPXUD6ZJ3RBA64RANCNFSM4JNTRNQA> .

[SamueleSoraggi]

I wonder if you might have an ordering problem of the positions into your data. But it seems unlikely this is the problem.
It is just weird because I have never seen this happening. Now I am looking also in the abbababa2 script.

[joanocha]

Would it help if I send you the region file?

To give you all relevant info:
If I do 5M or 2Kb blocks with the same -rf I always get the chr4 error for a name list of 6 pops.
If I just use 4 pops in the comparison the error of the negative window appears somewhere else.
My angsd version is 0.922

[SamueleSoraggi]

I think your region file is fine. But what did you write in the -sizeFile file and in the list of bam files?

[joanocha]

Here is the example for 4 populations where last pop is the outgroup.

sizefile.txt
30
25
10
5

list_bams.txt:
/path/to/sample1.bam
/path/to/sample2.bam
(...)

the list is ordered so that first 30 sample.bam are from pop1, followed by 25 pop2 etc and last 5 bams are the outgroup.

I just followed the workflow of the wiki.

I also did abbababa1 with the same bam list and I don't get bugs. So the bam list seems good. It's also the bam list I use for all angsd analysis which seem to be ok.

[ANGSD]

Has this issue been resolved?