Troubleshoot Running Coveragestats.Py In Tablet

Question

[lzwright]

I recently downloaded Tablet to assess read coverage of my de novo assembly of RNA-Seq of marine snail venom duct tissue. I have tried to run a script supplied by Tablet, coveragestats.py, which will give me: number of coverage values (=contig length), average depth, percent of bases with a depth greater than zero (percent coverage), and maximum depth.

I get the following error when I try to run the script

Traceback (most recent call last):
File "coveragestats.py", line 91, in <module>
stats = '\t'.join(templateStr) % tuple([func(coverages) for func in funcList])
File "coveragestats.py", line 71, in average
return float(sum(ints))/float(len(ints))
ZeroDivisionError: float division by zero

I am fairly certain this is happening because when I load my .bam and reference files to Tablet I get the following warning:

The index for this BAM file is reporting a total read count of zero. Although this may be correct, it is more likely that the index file was generated prior to read metrics being available. It is advisable to recreate this index using samtools 0.1.8 or higher.

I have verified my version of samtools as being the latest version and have reindexed my .bam file several times. Still the same problem (read count 0). However I do see the number of reads associated with my contigs in the Tablet header (and they are displayed), just not in the actual list of contigs in the left panel where this info should be associated to each contig.

Thanks in advance for any help.

[Peter]

Does 'samtools idxstats your_file.bam' confirm no mapped reads? This is the information Tablet is using to get the read count.

And what exactly is your version of samtools?

[lzwright]

re version problem is I am now home (I can answer this with more certainty tomorrow) but, I made sure I had the at least 0.1.18 by running sudo apt-get update and then sudo apt-get install samtools. then for good measure, since I saw there was a .19 version, I git cloned that on to my machine and indexed again. here at home I also have ubuntu running and samtools is most assuredly 0.1.18.

I did transport the files home with me (myfile.bam, myfile_sorted.bam and myfile_sorted.bam.bai) and when I run idxstats command I get this

~/trinity_bowtie_map$ samtools idxstats Trinity_bowtie.bam
[bam_index_load] fail to load BAM index.
[bam_idxstats] fail to load the index.

(then I reindexed one more time for good measure and same result) incidentally this is a bowtie2 mapped file, I have also done mapping with bwa-mem so I could see what happens with this .bam file as well

to make sure I am doing what you want, I should be running idxstats on .bam file and not sorted file?

[Istvan Albert]

There you have your answer that is not a properly sorted and indexed BAM file. Not even samtools can read it.

Go back to the beginning and redo the BAM file generation.

Also BAM files should always be sorted otherwise they can't be indexed.

[lzwright]

I have already tried rebuilding and yesterday I followed Tablet protocol exactly for generating BAM file

samtools faidx reference.fasta
samtools view -b -S -t reference.fasta.fai -o assembly.bam assembly.sam
samtools sort assembly.bam assembly_sorted.bam
samtools index assembly_sorted.bam

it all seemed to work fine....

[lzwright]

This was approach I had taken to building BAM file previously

samtools faidx dmel-all-chromosome-r5.37.fasta
samtools import dmel-all-chromosome-r5.37.fasta.fai RAL357_bwa.sam RAL357_bwa.bam
samtools sort RAL357_bwa.bam RAL357_bwa.sorted
samtools index RAL357_bwa.sorted.bam

both approaches produced the expected files

[Istvan Albert]

now you need to test that the process completed as expected, for example what do commands such as

samtools view -c -F 4 <bamfile>
samtools view <bamfile> contig:start-stop
samtools idxstats <bamfile>

produce? (substitute real values for contig/start/stop etc) those need to run without a hitch to assure that the bamfile is indeed correct.

[lzwright]

Once again telling me not indexed. weird. will try whole routine with my BWA-mem mapped file and see what happens?

-rw------- 1 liz liz 21G Aug 16 11:16 Trinity_bowtie.bam
-rw------- 1 liz liz 19G Aug 16 13:04 Trinity_bowtie_sorted.bam
-rw------- 1 liz liz 10M Aug 16 20:03 Trinity_bowtie_sorted.bam.bai

liz@liz-ThinkCentre-Edge-91Z:~/trinity_bowtie_map$ samtools view -c -F 4 Trinity_bowtie.bam
67308534

liz@liz-ThinkCentre-Edge-91Z:~/trinity_bowtie_map$ samtools view Trinity_bowtie.bam comp43206_c2_seq2:488-856
[bam_index_load] fail to load BAM index.
[main_samview] random alignment retrieval only works for indexed BAM files.

liz@liz-ThinkCentre-Edge-91Z:~/trinity_bowtie_map$ samtools idxstats Trinity_bowtie.bam
[bam_index_load] fail to load BAM index.
[bam_idxstats] fail to load the index.

[Istvan Albert]

the answer is simple - as the errors says you don't actually have an indexed BAM file.

now why and how to fix that is for you to sort out.

[lzwright]

great. still working on it, the mystery continues....

for the record I am using seqtk to extract 50K reads from my /1 and /2 fastq files (keeping pairing) and will map these to my assembly for sam to bam conversion, sort and index. will make file easier to pass around to willing victims (I mean candidates)

[lzwright]

An update: I was running idxstats on the wrong file (namely the .bam file and not the sorted.bam file). So as it turns out my BAM file is sorted and indexed correctly. I was able to take a different approach from Tablet (using edgeR and bedtools) http://ged.msu.edu/angus/tutorials-2013/rnaseq_bwa_counting.html in order to get read counts on my contigs. Am in discussion with Tablet developer to see why read counts will not display properly.

[Peter]

Duplicate post on SEQanswers, http://seqanswers.com/forums/showthread.php?t=32870

[Istvan Albert]

Run samtools and post the first two lines it prints, like so

$ samtools   

Program: samtools (Tools for alignments in the SAM format)
Version: 0.1.18 (r982:295)

[lzwright]

OK. will report in more detail tomorrow (in meantime see above)

Answer 1

[lzwright]

OK so this is not a problem with my bam file, or the sorting and indexing process. It is a Tablet issue. I downloaded Tablet for linux and this is where I was getting error message "read count is zero" when loading assembly and reference files. So I downloaded Tablet for mac, and loaded exactly the same files. It works perfectly, read counts are loaded. I am in touch with the developer to bring this to his attention and see if it can be fixed.

Answer 2

[Istvan Albert]

Just to mark this as answered - the comment seems to indicate there were problems with the BAM file to begin with.