samtools mpileup error

Question

[teresayahu]

I am using samtools function mpileup to get reads in sites. The command I used is:

samtools mpileup -l hg19.position -f hg19.fa -Q 20 -q 20 -I file.bam

in the command, hg19.position is position-based site file.

When running these is an error printed on screen: ...

[E::faidx_fetch_seq] The sequence "1" not found
[E::faidx_fetch_seq] The sequence "1" not found
[E::faidx_fetch_seq] The sequence "1" not found
[E::faidx_fetch_seq] The sequence "1" not found
[E::faidx_fetch_seq] The sequence "1" not found
[E::faidx_fetch_seq] The sequence "1" not found

...

May I ask the reason of this error and how could I fix it? Thanks! Ya

[Kevin Blighe]

Have you indexed the hg19.fa file using the same version of samtools as per your command above and using the samtools faidx command?

[BENYEBUGA]

• I'm getting the same error as the OP "[E::faidx_fetch_seq] The sequence "1" not found"

I tried the index command you suggested,

samtools faidx hg19.fa

After pressing Enter there's a few second pause, followed by a new line:

ben@DESKTOP-EMRLM4V:~$ samtools faidx hg19.fa

ben@DESKTOP-EMRLM4V:~$

... Does this mean the fa file is now indexed?

If so, why do I keep getting the same error, "[E::faidx_fetch_seq] The sequence "1" not found", when executing my initial command:

bcftools mpileup -Ou -f hg19.fa I19139.hg19.bam | bcftools call -vmO z -o Nubians.vcf.gz

Thanks for any help,

[BENYEBUGA]

For added context, here's a view of some lines from the BAM file using the samtools view H command:

@HD     VN:1.5  GO:none SO:coordinate
@SQ     SN:1    LN:249250621
@SQ     SN:2    LN:243199373
@SQ     SN:3    LN:198022430
@SQ     SN:4    LN:191154276
@SQ     SN:5    LN:180915260
@SQ     SN:6    LN:171115067
@SQ     SN:7    LN:159138663
@SQ     SN:8    LN:146364022
@SQ     SN:9    LN:141213431
@SQ     SN:10   LN:135534747
@SQ     SN:11   LN:135006516
@SQ     SN:12   LN:133851895
@SQ     SN:13   LN:115169878
@SQ     SN:14   LN:107349540
@SQ     SN:15   LN:102531392
@SQ     SN:16   LN:90354753
@SQ     SN:17   LN:81195210
@SQ     SN:18   LN:78077248
@SQ     SN:19   LN:59128983
@SQ     SN:20   LN:63025520
@SQ     SN:21   LN:48129895
@SQ     SN:22   LN:51304566
@SQ     SN:X    LN:155270560
@SQ     SN:Y    LN:59373566
@SQ     SN:MT   LN:16569

[Pierre Lindenbaum]

I would say you're bam or your reference using the 'chr1' notation, not just '1'

[teresayahu]

Thanks! both the bam and reference files are using the '1'. Is there a way that I can change bam notation from '1' to 'chr1'?

[GenoMax]

Please use ADD COMMENT/ADD REPLY when responding to existing posts to keep threads logically organized.

[Pierre Lindenbaum]

what about hg19.position ??