fix bamdash to deal with BAM file with no mapped reads (#45)

* fix on bamdash process to deal with no mapped reads

* adjusting coveragePlot_out_ch channel

* Update sars-cov-2.params

setup.py

@@ -1,7 +1,7 @@
 from setuptools import setup
 
 setup(name='ViralFlow',
-      version='1.0.0',
+      version='1.0.1',
       description='''
       Workflows for viral genome Assembly at FioCruz/IAM
       ''',

test_files/denv.params

@@ -10,7 +10,7 @@ trimLen 0
 refGenomeCode NC_001474.2
 referenceGFF null
 referenceGenome null
-nextflowSimCalls null
+nextflowSimCalls 6
 fastp_threads 1
 bwa_threads 1
 mafft_threads 1

test_files/sars-cov-2.params

@@ -11,7 +11,7 @@ trimLen 0
 refGenomeCode null
 referenceGFF null
 referenceGenome null
-nextflowSimCalls null
+nextflowSimCalls 6
 fastp_threads 1
 bwa_threads 1
 mafft_threads 1

vfnext/containers/build_containers.py

@@ -58,13 +58,13 @@ def build_container(container, command):
     print("\nSome containers failed to build. Please check the error messages above.")
 
 if success:
-    print("\nExecuting additional steps:")
+    print("\nExecuting additional steps:\n")
 
-    print("  > Loading sars-cov2 nextclade dataset...")
+    print("  > Loading sars-cov2 nextclade dataset...\n")
     nextclade_command = "singularity exec -B nextclade_dataset/sars-cov-2:/tmp nextclade:2.4.sif nextclade dataset get --name 'sars-cov-2' --output-dir '/tmp'"
     try:
         subprocess.check_call(nextclade_command, shell=True)
-        print("    > Done <")
+        print("    > Done <\n")
     except subprocess.CalledProcessError as e:
         print("    > Failed <")
         print(f"Error: {e}")
@@ -85,10 +85,14 @@ def build_container(container, command):
     if not os.path.exists(unsquashfs_desired_location):
         print("\n\033[91mError:\n  > unsquashfs executable not found at expected location. You should create a symbolic link using the following command:\033[0m")
         unsquashfs_location = os.path.join(os.environ["HOME"], "miniconda3/envs/viralflow/bin/unsquashfs")
-        print(f"    > sudo ln -s {unsquashfs_location} /usr/local/bin/unsquashfs")
+        print(f"   >  sudo ln -s {unsquashfs_location} /usr/local/bin/unsquashfs\n")
+        print(f"  > If the first symbolic link does not solve the error you can alternatively create a symbolic link using the following command:")
+        print(f"   >  sudo ln -s /usr/bin/unsquashfs /usr/local/bin/unsquashfs\n")
+
         # Print a message indicating the unsquashfs location
-        print(f"  > unsquashfs executable is located at {unsquashfs_location}.")
+        print(f"  > unsquashfs executable is located at {unsquashfs_location} or at /usr/bin/unsquashfs\n")
+        print(f"  > After create the symbolic link for unsquashfs, please run 'viralflow -build_containers' command again to finish the additional steps necessary to run ViralFlow correctly.")
 
 if success:
     print("\nAll steps from '-build_containers' completed successfully. You can test ViralFlow using the following command:")
-    print("  > viralflow -run --params_file test_files/sars-cov-2.params")
+    print("   > viralflow -run --params_file test_files/sars-cov-2.params")

vfnext/main.nf

@@ -42,7 +42,7 @@ ANSI_RESET = "\033[0m"
 
 log.info """
   ===========================================
-  VFNEXT v1.0.0
+  VFNEXT v1.0.1
   parameters:
   -------------------------------------------
   --inDir            : ${params.inDir}
@@ -147,6 +147,10 @@ workflow {
   //Rendering the depth coverage plot
   coveragePlot(align2ref.out.regular_output,
               ref_gcode)
+  // Check if there are mapped reads
+  coveragePlot_out_ch = coveragePlot.out.result
+  coveragePlot_out_ch
+    | view(it -> log.warn("${it.text}"))
 
   if ((params.writeMappedReads == true)){
     // write mapped reads

vfnext/modules/generatePlots.nf

@@ -5,9 +5,10 @@ process coveragePlot {
 
       tuple val(sample_id), path(bam_files), path(bai_files)
       val(genome_code)
+
     output:
-        path("*coveragePlot*")
-
+        path("*coveragePlot*"), optional: true
+        path("coveragePlot_result.txt"), optional: true, emit: result, hidden: true
     script:
 
     bam = bam_files[0].toString()
@@ -16,17 +17,29 @@ process coveragePlot {
     html = "${sample_id}_coveragePlot.html"
     png = "${sample_id}_coveragePlot.png"
     svg = "${sample_id}_coveragePlot.svg"
-
+    
     """
-    bamdash -b ${bam} -r ${genomecode} -c ${depth} -e svg
-    mv ${genomecode}_plot.svg ${svg}
+    #!/usr/bin/env python
+    # ----- import libraries -------------------------------------------------
+    import pysam
+    import subprocess
+    
+    n_mapped_reads = pysam.AlignmentFile('${bam}','rb').count()
+    if n_mapped_reads > 0:
+      subprocess.run([f"bamdash -b ${bam} -r ${genomecode} -c ${depth} -e svg"], shell=True)
+      subprocess.run([f"mv ${genomecode}_plot.svg ${svg}"], shell=True)
 
-    bamdash -b ${bam} -r ${genomecode} -c ${depth} -e png
-    mv ${genomecode}_plot.png ${png}
+      subprocess.run([f"bamdash -b ${bam} -r ${genomecode} -c ${depth} -e png"], shell=True)
+      subprocess.run([f"mv ${genomecode}_plot.png ${png}"], shell=True)
 
-    bamdash -b ${bam} -r ${genomecode} -c ${depth}
-    mv ${genomecode}_plot.html ${html}
-    """
+      subprocess.run([f"bamdash -b ${bam} -r ${genomecode} -c ${depth}"], shell=True)
+      subprocess.run([f"mv ${genomecode}_plot.html ${html}"], shell=True)
+   
+    else:
+      result = "No mapped reads were found in the sorted BAM file for sample ${sample_id}. The coverage plot will not be generated for it."
+      with open('coveragePlot_result.txt', 'w') as f:
+        f.write(result)
+      """
 }
 
 process snpPlot {

vfnext/nextflow.config

@@ -26,7 +26,7 @@ params {
   referenceGenome=null
 
   // set available resources for nextflow
-  nextflowSimCalls = null
+  nextflowSimCalls = 6
 
   // set quality params
   mapping_quality = 30

wrapper/viralflow

@@ -8,7 +8,7 @@ from sys import exit
 
 __author__ = "Antonio Marinho da Silva Neto"
 __license__ = "GPL"
-__version__ = "1.0.0"
+__version__ = "1.0.1"
 __maintainer__ = "Antonio Marinho da Silva Neto"
 __email__ = "[email protected]"
 
@@ -124,7 +124,7 @@ if args.run == True:
     if args.params_file == None:
         print("ERROR: no params_file was provided")
         exit(1)
-    print("ViralFlow v1.0.0")
+    print("ViralFlow v1.0.1")
     wrapper.run_vfnext(VF_ROOT_PATH, args.params_file)
     exit(0)