marbl / VALET

A pipeline for detecting mis-assemblies in metagenomic assemblies.
MIT License
14 stars 4 forks source link

generate summary table fails #2

Closed jbadomics closed 7 years ago

jbadomics commented 8 years ago

I am validating a metagenomic assembly with a single FASTA file of reads and keep getting this error. The contigs are numbered sequentially in the final assembly, e.g.:

>951_contig1
>951_contig2
...
>951_contig10682

Assembly here

The pipeline keeps failing when trying to generate the summary table. Full command and STDOUT are copied below. I am running (on a cluster):

Python 2.7.9 bowtie2 2.2.4 samtools 1.2 numpy 1.9.3 bedtools 2.24

$ python ~/Software/cmhill_VALET/src/py/valet.py -a 951.final.renamed.headers.fasta -r ../02-pear/951pear.assembled.fasta
###########################################################################
PROCESSING ASSEMBLY: asm_0 (951.final.renamed.headers.fasta)
###########################################################################
COMMAND:     reapr facheck 951.final.renamed.headers.fasta output/asm_0/assembly_facheck
RESULTS:     output/asm_0/assembly_facheck
---------------------------------------------------------------------------
STEP:    FILTERING ASSEMBLY CONTIGS LESS THAN 1000 BPs
RESULTS:     output/asm_0/filtered_assembly.fasta
---------------------------------------------------------------------------
STEP:    ALIGNING READS
COMMAND:     bowtie2-build /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/output/asm_0/filtered_assembly.fasta /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/output/asm_0/indexes/temp_2ucgDE
COMMAND:     bowtie2 -a -x /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/output/asm_0/indexes/temp_2ucgDE -q -U ../02-pear/951pear.assembled.fasta --very-sensitive -a --reorder -p 8 --un /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/output/asm_0/unaligned_reads/unaligned.reads -S output/asm_0/sam/library.sam
---------------------------------------------------------------------------
STEP:    RUNNING SAMTOOLS
COMMAND:     samtools view -F 0x100 -bS output/asm_0/sam/library.sam
COMMAND:     samtools sort output/asm_0/bam/library.bam output/asm_0/bam/sorted_library
COMMAND:     samtools mpileup -C50 -A -f output/asm_0/filtered_assembly.fasta output/asm_0/bam/sorted_library.bam
RESULTS:     output/asm_0/coverage/mpileup_output.out
COMMAND:     samtools index output/asm_0/bam/sorted_library.bam
---------------------------------------------------------------------------
STEP:    CALCULATING CONTIG COVERAGE
RESULTS:     output/asm_0/coverage/temp.cvg
---------------------------------------------------------------------------
STEP:    PARTITIONING COVERAGE FILE
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/split_pileup.py -p output/asm_0/coverage/mpileup_output.out -c 8
---------------------------------------------------------------------------
STEP:    DEPTH OF COVERAGE
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/depth_of_coverage.py -m output/asm_0/coverage/mpileup_output.out -w 501 -o output/asm_0/coverage/errors_cov.bed -g -e -c 8
COMMAND:     bedtools sort -i output/asm_0/coverage/errors_cov.bed
Error: The requested file (output/asm_0/coverage/errors_cov.bed) could not be opened. Error message: (No such file or directory). Exiting!
RESULTS:     output/asm_0/coverage.bed
---------------------------------------------------------------------------
STEP:    BREAKPOINT
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/breakpoint_splitter.py -u /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/output/asm_0/unaligned_reads/ -o output/asm_0/breakpoint/split_reads/
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/breakpoint_finder.py -a output/asm_0/filtered_assembly.fasta -r output/asm_0/breakpoint/split_reads/ -b 50 -o output/asm_0/breakpoint/ -c output/asm_0/coverage/temp.cvg -p 8
COMMAND:     bedtools sort -i output/asm_0/breakpoint/interesting_bins.bed
Error: The requested file (output/asm_0/breakpoint/interesting_bins.bed) could not be opened. Error message: (No such file or directory). Exiting!
RESULTS:     output/asm_0/breakpoint/../breakpoints.bed
---------------------------------------------------------------------------
STEP:    SUMMARY
RESULTS:     output/asm_0/summary.bed
RESULTS:     output/asm_0/suspicious.bed
Traceback (most recent call last):
  File "/home/bonddr/badalame/Software/cmhill_VALET/src/py/valet.py", line 1169, in <module>
    main()
  File "/home/bonddr/badalame/Software/cmhill_VALET/src/py/valet.py", line 239, in main
    contig_lengths, contig_abundances, final_misassemblies)
  File "/home/bonddr/badalame/Software/cmhill_VALET/src/py/valet.py", line 1008, in generate_summary_table
    table_file.write(contig + '\t' + str(filtered_contig_lengths[contig]) + '\t' + str(contig_abundances[contig]) + '\t' + \
KeyError: '951_contig6066'
cmhill-zz commented 8 years ago

Thanks for the detailed output. In the meantime, would you also have a link to a sample (or the whole) of the reads?

jbadomics commented 8 years ago

You should be able to pull them from here: https://www.dropbox.com/s/3zp6wzv8shd96j8/951pear.assembled.fasta.gz?dl=0

cmhill-zz commented 8 years ago

The issue was that the input reads were erroneously treated as FASTQ. I've fixed the issue in the latest commits. Add a -f to specify that the input reads are in FASTA format. Thanks for bringing this to my attention!

python src/py/valet.py -a data/951.final.renamed.headers.fasta.txt -r data/951pear.assembled.fasta -f --threads 16
jbadomics commented 8 years ago

Awesome, thanks! The pipeline completes successfully now, except I can't tell for sure if if the output I'm getting is correct - the log is saying a particular .bed file could not be found, and the overall coverage.bed file is empty (as are summary.bed and suspicious.bed). So the final summary.tsv file is showing no assembly errors (would be great if this is in fact true), but I'm not convinced this is the case since some of the files are empty. Any ideas or suggestions? Results here (4.5 GB total) if it would help.

$ python /home/bonddr/badalame/Software/cmhill_VALET/src/py/valet.py -a 951.final.renamed.headers.fasta -r ../02-pear/951pear.assembled.fasta -f -p 24 -o VALET_out -g 5 -z 1000 -b 50
###########################################################################
PROCESSING ASSEMBLY: asm_0 (951.final.renamed.headers.fasta)
###########################################################################
---------------------------------------------------------------------------
STEP:    FILTERING ASSEMBLY CONTIGS LESS THAN 1000 BPs
RESULTS:     VALET_out/asm_0/filtered_assembly.fasta
---------------------------------------------------------------------------
STEP:    ALIGNING READS
COMMAND:     bowtie2-build /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/VALET_out/asm_0/filtered_assembly.fasta /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/VALET_out/asm_0/indexes/temp_b5OHaI
COMMAND:     bowtie2 -a -x /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/VALET_out/asm_0/indexes/temp_b5OHaI -f -U ../02-pear/951pear.assembled.fasta --very-sensitive -a --reorder -p 24 --un /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/VALET_out/asm_0/unaligned_reads/unaligned.reads -S VALET_out/asm_0/sam/library.sam
---------------------------------------------------------------------------
STEP:    RUNNING SAMTOOLS
COMMAND:     samtools view -F 0x100 -bS VALET_out/asm_0/sam/library.sam
COMMAND:     samtools sort VALET_out/asm_0/bam/library.bam VALET_out/asm_0/bam/sorted_library
COMMAND:     samtools mpileup -C50 -A -f VALET_out/asm_0/filtered_assembly.fasta VALET_out/asm_0/bam/sorted_library.bam
RESULTS:     VALET_out/asm_0/coverage/mpileup_output.out
COMMAND:     samtools index VALET_out/asm_0/bam/sorted_library.bam
---------------------------------------------------------------------------
STEP:    CALCULATING CONTIG COVERAGE
RESULTS:     VALET_out/asm_0/coverage/temp.cvg
---------------------------------------------------------------------------
STEP:    PARTITIONING COVERAGE FILE
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/split_pileup.py -p VALET_out/asm_0/coverage/mpileup_output.out -c 24
---------------------------------------------------------------------------
STEP:    DEPTH OF COVERAGE
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/depth_of_coverage.py -m VALET_out/asm_0/coverage/mpileup_output.out -w 501 -o VALET_out/asm_0/coverage/errors_cov.bed -g -e -c 24
COMMAND:     bedtools sort -i VALET_out/asm_0/coverage/errors_cov.bed
Error: The requested file (VALET_out/asm_0/coverage/errors_cov.bed) could not be opened. Error message: (No such file or directory). Exiting!
RESULTS:     VALET_out/asm_0/coverage.bed
---------------------------------------------------------------------------
STEP:    BREAKPOINT
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/breakpoint_splitter.py -u /panfs/roc/groups/1/bonddr/badalame/CoDL/DDH951/14-VALET/VALET_out/asm_0/unaligned_reads/ -o VALET_out/asm_0/breakpoint/split_reads/
COMMAND:     /home/bonddr/badalame/Software/cmhill_VALET/src/py/breakpoint_finder.py -a VALET_out/asm_0/filtered_assembly.fasta -r VALET_out/asm_0/breakpoint/split_reads/ -b 50 -o VALET_out/asm_0/breakpoint/ -c VALET_out/asm_0/coverage/temp.cvg -p 24 -f
COMMAND:     bedtools sort -i VALET_out/asm_0/breakpoint/interesting_bins.bed
RESULTS:     VALET_out/asm_0/breakpoint/../breakpoints.bed
---------------------------------------------------------------------------
STEP:    SUMMARY
RESULTS:     VALET_out/asm_0/summary.bed
RESULTS:     VALET_out/asm_0/suspicious.bed
RESULTS:     VALET_out/asm_0/summary.tsv
###########################################################################
GENERATING ASSEMBLY COMPARISON PLOTS
###########################################################################
COMMAND:     Rscript /home/bonddr/badalame/Software/cmhill_VALET/src/R/compare_assemblies.R VALET_out/asm_0/summary.tsv asm_0 VALET_out/comparison_plots
RESULTS:     VALET_out/comparison_plots.pdf
nate-d-olson commented 8 years ago

There was an error with the coverage analysis (see log file output below), but it was not reported during the run. The error should have been reported, I will work on fixing the bug. The error is associated with numpy, try re-installing numpy then rerunning VALET and let me know if it works.

Here is the text from the coverage/error.log file

Traceback (most recent call last):
  File "/home/bonddr/badalame/Software/cmhill_VALET/src/py/depth_of_coverage.py", line 7, in <module>
    import numpy as np
  File "/home/bonddr/badalame/local/lib/python2.7/site-packages/numpy-1.9.2-py2.7-linux-x86_64.egg/numpy/__init__.py", line 170, in <module>
    from . import add_newdocs
  File "/home/bonddr/badalame/local/lib/python2.7/site-packages/numpy-1.9.2-py2.7-linux-x86_64.egg/numpy/add_newdocs.py", line 13, in <module>
    from numpy.lib import add_newdoc
  File "/home/bonddr/badalame/local/lib/python2.7/site-packages/numpy-1.9.2-py2.7-linux-x86_64.egg/numpy/lib/__init__.py", line 8, in <module>
    from .type_check import *
  File "/home/bonddr/badalame/local/lib/python2.7/site-packages/numpy-1.9.2-py2.7-linux-x86_64.egg/numpy/lib/type_check.py", line 11, in <module>
    import numpy.core.numeric as _nx
  File "/home/bonddr/badalame/local/lib/python2.7/site-packages/numpy-1.9.2-py2.7-linux-x86_64.egg/numpy/core/__init__.py", line 6, in <module>
    from . import multiarray
ImportError: /home/bonddr/badalame/local/lib/python2.7/site-packages/numpy-1.9.2-py2.7-linux-x86_64.egg/numpy/core/multiarray.so: undefined symbol: PyUnicodeUCS2_AsASCIIString
nate-d-olson commented 7 years ago

I'm closing this issue. Let us know if you have any other problems with VALET.