Question: Troubleshoot Running Coveragestats.Py In Tablet
0
gravatar for lzwright
5.7 years ago by
lzwright150
NYC
lzwright150 wrote:

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.

• 3.0k views
ADD COMMENTlink modified 13 months ago by Biostar ♦♦ 20 • written 5.7 years ago by lzwright150

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?

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by Peter5.8k

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?

ADD REPLYlink modified 5.7 years ago by Istvan Albert ♦♦ 80k • written 5.7 years ago by lzwright150
1

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.

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by Istvan Albert ♦♦ 80k

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....

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by lzwright150

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

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by lzwright150

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.

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by Istvan Albert ♦♦ 80k

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.

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by lzwright150

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.

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by Istvan Albert ♦♦ 80k

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)

ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by lzwright150

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.

ADD REPLYlink written 5.7 years ago by lzwright150

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

ADD REPLYlink written 5.7 years ago by Peter5.8k

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)
ADD REPLYlink modified 5.7 years ago • written 5.7 years ago by Istvan Albert ♦♦ 80k

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

ADD REPLYlink written 5.7 years ago by lzwright150
1
gravatar for lzwright
5.7 years ago by
lzwright150
NYC
lzwright150 wrote:

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.

ADD COMMENTlink written 5.7 years ago by lzwright150
0
gravatar for Istvan Albert
5.7 years ago by
Istvan Albert ♦♦ 80k
University Park, USA
Istvan Albert ♦♦ 80k wrote:

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

ADD COMMENTlink written 5.7 years ago by Istvan Albert ♦♦ 80k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1617 users visited in the last hour