Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[loadChr] Error loading fasta info from chr #617

Open
dlema82 opened this issue Mar 14, 2024 · 0 comments
Open

[loadChr] Error loading fasta info from chr #617

dlema82 opened this issue Mar 14, 2024 · 0 comments

Comments

@dlema82
Copy link

dlema82 commented Mar 14, 2024

Hello, I am trying to get some Fst stats,
I've been all over google and github and sourceforge and biostars and etc. I cannot find a solution. Please help me find a workaround, or explain what I am doing wrong. What am I missing here????

I keep getting this error:

[loadChr] Error loading fasta info from chr:'KP202265.1_Panthera_pardus'
-> Problem with length of fastafile vs length of chr in BAM header
-> Chromosome name: 'KP202265.1_Panthera_pardus' length from BAM header:24850 length from fai file:0
Trying to access fasta efter end of chromsome+200:KP202265.1_Panthera_pardus/KP202265.1_Panthera_pardus pos=250 ref_len=0
Trying to access fasta efter end of chromsome+200:KP202265.1_Panthera_pardus/KP202265.1_Panthera_pardus pos=433 ref_len=0

I have tried using the -checkBamHeaders flag, but I get the error regardless. I've also tried samtools reheader to try to make them compatible, but it doesn't work either, I still get the error.

This is for haploid mitochondrial data. Here is my ref fai:
KP202265.1_Panthera_pardus 24850 28 70 71

And the header for the ref fasta:

KP202265.1_Panthera_pardus

And here is the head of my bam files called with:
samtools view 2469_mtgenome.mapped.srt.dedup.KP.bam | head -n 2
A00312:75:HLLFVDSXX:4:1215:13810:31923 163 KP202265.1_Panthera_pardus 160 4M1I145M = 170 330 CGATGGGTTAATGACTAATCAGCCCATGATCACACATAACTGTGGTGTCATGCATTTGGTATCTTTAATTTTTAGGGGGTCGAACTTGCTATGACTCAGCTATGACCTAAAGGTCCTGACTCAGTCAAATATAATGTAGCTGGACTTATT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFF,FFFFFFFFFF:FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF:FFFFFFFFFFFFFFFFFFFFFFF:FFF MC:Z:96M8D43M3D11M MD:Z:9G2G30A17T87 PG:Z:MarkDuplicates NM:i:5 AS:i:126 XS:i:0
A00312:75:HLLFVDSXX:4:2351:6253:11318 145 KP202265.1_Panthera_pardus 160 4M1I145M = 16482 16334 CGATGGGTTAATGACTAATCAGCCCATGATCACACATAACTGTGGTGTCATGCATTTGGTATCTTTAATTTTTAGGGGGTCGAACTTGCTATGACTCAGCTATGACCTAAAGGTCCTGACTCAGTCAAATATAATGTAGCTGGACTTATT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF MC:Z:150M MD:Z:9G2G30A17T87 PG:Z:MarkDuplicates NM:i:5 AS:i:126 XS:i:0

Here is my prompt:
ANGSD="/nas3/dlema/angsd/angsd"
REAL="/nas3/dlema/angsd/misc/realSFS"
ANC="/nas3/dlema/NC028302.1_Pleo_mtgenome.fasta"
REF="/nas3/dlema/Ppardus_mtgenome_li_etal_plateau.fa"
FILTERS="-minMapQ 30 -minQ 20"
OPT=" -dosaf 1 -gl 1" ##I have tried with and without -checkBamHeaders and I get the error regardless

$ANGSD -b africa.filelist -anc $ANC -out total_africa $FILTERS $OPT -ref $REF
$ANGSD -b east.filelist -anc $ANC -out east $FILTERS $OPT -ref $REF
$ANGSD -b west.filelist -anc $ANC -out west $FILTERS $OPT -ref $REF
$ANGSD -b west1.filelist -anc $ANC -out west1 $FILTERS $OPT -ref $REF
$ANGSD -b northeast.filelist -anc $ANC -out northeast $FILTERS $OPT -ref $REF
$ANGSD -b south.filelist -anc $ANC -out south $FILTERS $OPT -ref $REF
$ANGSD -b east1.filelist -anc $ANC -out east1 $FILTERS $OPT -ref $REF
$ANGSD -b east2.filelist -anc $ANC -out east2 $FILTERS $OPT -ref $REF

$REAL total_africa.saf.idx > total_africa.sfs
$REAL east.saf.idx > east.sfs
$REAL west.saf.idx > west.sfs
$REAL west1.saf.idx > west1.sfs
$REAL northeast.saf.idx > northeast.sfs
$REAL south.saf.idx > south.sfs
$REAL east1.saf.idx > east1.sfs
$REAL east2.saf.idx > east2.sfs

$REAL east.saf.idx west.saf.idx > east_west.ml
$REAL east.saf.idx west1.saf.idx > east_west1.ml
$REAL east.saf.idx northeast.saf.idx > east_northeast.ml
$REAL east.saf.idx south.saf.idx > east_south.ml
$REAL east1.saf.idx east2.saf.idx > east1_east2.ml
$REAL west.saf.idx west1.saf.idx > west_west1.ml
$REAL west.saf.idx northeast.saf.idx > west_northeast.ml
$REAL west.saf.idx south.saf.idx > west_south.ml
$REAL west1.saf.idx northeast.saf.idx > west1_northeast.ml
$REAL west1.saf.idx south.saf.idx > west1_south.ml
$REAL northeast.saf.idx south.saf.idx > northeast_south.ml

$REAL fst index east.saf.idx west.saf.idx -sfs east_west.ml -fstout east.west
$REAL fst index east.saf.idx west1.saf.idx -sfs east_west1.ml -fstout east.west1
$REAL fst index east.saf.idx northeast.saf.idx -sfs east_northeast.ml -fstout east.northeast
$REAL fst index east.saf.idx south.saf.idx -sfs east_south.ml -fstout east.south
$REAL fst index east1.saf.idx east2.saf.idx -sfs east1_east2.ml -fstout east1.east2
$REAL fst index west.saf.idx west1.saf.idx -sfs west_west1.ml -fstout west.west1
$REAL fst index west.saf.idx northeast.saf.idx -sfs west_northeast.ml -fstout west.northeast
$REAL fst index west.saf.idx south.saf.idx -sfs west_south.ml -fstout west.south
$REAL fst index west1.saf.idx northeast.saf.idx -sfs west1_northeast.ml -fstout west1.northeast
$REAL fst index west1.saf.idx south.saf.idx -sfs west1_south.ml -fstout west1.south
$REAL fst index northeast.saf.idx south.saf.idx -sfs northeast_south.ml -fstout northeast.south

Thank you!!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant