Hi, I'm trying to run bam-readcount with region information (as a workaround to the 8000 read limitation). I use version 0.8 cloned from the repository.
The reference FASTA file (reference_1.fa) has this header:
It's indexed: bwa index reference_1.fa
I run the alignment:
bwa mem -t 8 reference_1.fa -p input.fastq | samtools sort -@8 -O BAM -o alignment_reference_1.bam -
Output is indexed:
samtools index alignment_reference_1.bam
At this point, I need to get base counts at each position in the reference. I created a separate file with region information (reference_1_regions.txt) that contains one line:
reference_1 1 35911 (tab-delimited)
The actual command:
bam-readcount -w 1 -d 1000000 -l reference_1_regions.txt -f reference_1.fa alignment_reference_1.bam > readcount_reference_1.txt
reference_1 not found in bam file. Region reference_1 1 35911 skipped.
I checked the bam file header with samtools view -H and here's the output:
@HD VN:1.5 SO:coordinate @SQ SN:reference_1:1-35911 LN:35911 @PG ID:bwa PN:bwa VN:0.7.17-r1188 CL:bwa mem -t 8 reference_1.fa -p input.fastq
If anyone has any suggestions about the reason for this error, please let me know. Is this a bug in bam-readcount or am I missing something?
Thank you in advance