Denovo bug fixes

README.md

@@ -116,12 +116,11 @@ For example:
   - Summary and plots of reads, alignments and the transcript assembly and annotation
   - One file per run with results from each sample merged
   
-####Transcript isoform table
+Note transcript isoform table is currently only available for the reference-guilded assembly.
-Note: Currently only available for the reference-guilded assembly
-The report includes a searchable filterable table of transcript isoforms<br>
 Note: If using a large dataset with 10s of thousands of predicted isoforms, to speed up searching, the table can be limited
 by read coverage, which can be set with `transcript_table_cov_thresh = 50` (default 50x)
 
+#### Output files
 Each sample will also have it's own directory containing the following (dependent on assembly approach and input options)
 * {sampleid}_gffcompare
   * Directory containing all output from gffcompare
@@ -133,9 +132,6 @@ Each sample will also have it's own directory containing the following (dependen
   - A transcriptome derived from the query reads + reference annotation 
 * merged_transcriptome.fas
   -A transcriptome made from the combined query reads and reference annotation
-* final_polished_cds.fas (de novo only)
-  * CDS sequences derived from reads
-  
   
 
 ## Useful links

bin/merge_gff.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+
+"""
+Merge and fix gff files.
+
+Merge multiple gff files into single file.
+Rename gene and transcript ids to avoid attribute conflicts from
+independently-created files.
+"""
+import argparse
+from pathlib import Path
+import re
+
+from natsort import natsorted
+
+
+def main(gff_files: str, outfile: str):
+    """Entry point."""
+    regx_id = re.compile(r'gene_id "STRG\.(\d+)"')
+
+    start = 1
+    with open(outfile, 'w') as fh:
+        for gff in gff_files:
+
+            text = Path(gff).read_text()
+            if start != 1:
+                # Strip headers
+                text = [x for x in text.splitlines() if not x.startswith('#')]
+                text = '\n'.join(text)
+
+            ids = natsorted(set(re.findall(regx_id, text)))
+            new_gene_ids = list(range(start, start + len(ids)))
+            id_map = dict(zip(ids, new_gene_ids))
+
+            for old_id, new_id in id_map.items():
+                text = text.replace(
+                    f'gene_id "STRG.{old_id}"',
+                    f'gene_id "STRG.{new_id}"')
+                text = text.replace(
+                    f'transcript_id "STRG.{old_id}.',
+                    f'transcript_id "STRG.{new_id}.')
+            fh.write(text)
+            fh.write('\n')
+            start += len(ids)
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--gff_files", help="gff files to merge",
+                        required=True, nargs='+')
+    parser.add_argument("--out_file", help="where to save merged files",
+                        required=True)
+    args = parser.parse_args()
+    main(args.gff_files, args.out_file)

bin/report.py

@@ -6,7 +6,7 @@ from collections import Counter, defaultdict, OrderedDict
 import math
 from pathlib import Path
 
-from aplanat import bars, lines
+from aplanat import bars, hist, lines
 from aplanat.components import fastcat
 from aplanat.components import simple as scomponents
 from aplanat.report import WFReport
@@ -648,7 +648,8 @@ def transcript_table(report, df_tmaps, covr_threshold):
     ### Query transcript table
 
     Low coverage transcripts are removed to speed up the table viewing. <br>
-    This can be set with the parameter `args.min_isoform_cov` in the config.
+    This can be set with the parameter `transcript_table_cov_thresh` in the
+    config.
       ''')
     section.plot(cov_plt)
     # Filter on converge threshold
@@ -724,17 +725,9 @@ def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
 
                 transcript_lens.append(tr_len)
 
-        num_bins = max(isoforms_per_gene)
+        bar_isos = hist.histogram(
-        if num_bins > 12:
+            [isoforms_per_gene], colors=[Colors.cerulean],
-            num_bins = 12
+            title="isoforms per gene")
-
-        hist_, bins = np.histogram(isoforms_per_gene, bins=num_bins)
-
-        bar_isos = bars.simple_bar([str(round(x_, 1)) for x_ in bins], hist_,
-                                   title="Isoforms per gene",
-                                   colors=Colors.cerulean,
-                                   x_axis_label='Num. isoforms',
-                                   y_axis_label='Num. genes')
 
         bar_isos.xaxis.major_label_orientation = math.pi / 2.8
         plots.append(bar_isos)

denovo_assembly.nf

@@ -39,16 +39,13 @@ process dump_clusters {
 
 process build_backbones {
     /*
-    This step can fail in what seems to at the racon stage giving an 'empty overlap error set' message
+    Build coding
-    As a temporary fix, do racon_cmd || true to prevent it crashing the pipeline, and then move on to next cluster
-
-    This process needs some work
     */
     label "isoforms"
 	input:
 		tuple val(sample_id), path(cluster_fq)
 	output:
-		tuple val(sample_id), path("*final_polished_cds.fa"), emit: polished_cds
+		tuple val(sample_id), path("*final_polished_cds.fa"), emit: polished_cds, optional: true
 	script:
 	def cluster_fq_bl = new BlankSeparatedList(cluster_fq)
 	"""
@@ -74,46 +71,52 @@ process build_backbones {
             # polish 1
             samgz=\${cln}_aln.sam.gz
             racon_cons=\${cln}_racon.fa
-        racon_cons1=\${cln}_racon1.fa
+            tmpcons=\${cln}_tmcons.fa
 
             minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$spoa_cons \$sample | gzip - > \$samgz
-        racon -t ${params.threads} --no-trimming -u -w 2000 \$sample \$samgz \$spoa_cons > \$racon_cons || true
+            exitcode=0
-        if [ -f \$racon_cons ];
+            racon -t ${params.threads} --no-trimming -u -w 2000 \$sample \$samgz \$spoa_cons > \$racon_cons || exitcode=1
+            if [[ \$exitcode -eq 0 ]];
                 then
-            cat \$racon_cons | seqkit replace -p ".*" -r cluster_\${cln} > \$racon_cons1
+                    # Rename consensus sequence name with cluster id
+                    cat \$racon_cons | seqkit replace -p ".*" -r cluster_\${cln} > \$tmpcons
+                    mv \$tmpcons \$racon_cons
             else
+                echo "Polishing failed for \${cln}"
                 continue
             fi
 
             # polish 2
-        samgz1=\${cln}_aln_1.sam.gz
+            minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons \$sample  | gzip - > \$samgz
-        racon_cons2=\${cln}_racon2.fa
+            exitcode=0
-        minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons1 \$sample  | gzip - > \$samgz1
+            racon -t ${params.threads}  -u --no-trimming \$sample \$samgz \$racon_cons > \$tmpcons || exitcode=1
-        racon -t ${params.threads}  -u --no-trimming \$sample \$samgz1 \$racon_cons1 > \$racon_cons2 || true
+            if [[ \$exitcode -eq 0 ]];
-        if [ -f \$racon_cons2 ];
                 then
                     echo "success polish 2"
+                    mv \$tmpcons \$racon_cons
             else
+                echo "Polishing failed for \${cln}"
                 continue
             fi
 
              # polish 3
-        samgz2=\${cln}_aln_2.sam.gz
+            minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons \$sample | gzip - > \$samgz
-        minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons2 \$sample | gzip - > \$samgz2
+            exitcode=0
-        EXITCODE=0
+            racon -t ${params.threads}  -u \$sample \$samgz \$racon_cons > \$tmpcons || exitcode=1
-        racon -t ${params.threads}  -u \$sample \$samgz2 \$racon_cons2 >> ${sample_id}_\${UNID}_final_polished_cds.fa || EXITCODE=\$?
+            if [[ \$exitcode -eq 0 ]];
-        if [ \$EXITCODE -eq 0 ];
                 then
-            echo "finished"
+                    cat \$tmpcons >> ${sample_id}_\${UNID}_final_polished_cds.fa
+                    echo "polishing 3 success"
+            else
+                echo "Polishing failed for \${cln}"
             fi
 
         done
     echo "Finished backbones"
 	"""
-
-
 }
 
+
 process merge_cds {
     label "isoforms"
     input:
@@ -121,7 +124,6 @@ process merge_cds {
     output:
         tuple val(sample_id), path("${sample_id}_cds.fa"), emit: final_polished_cds
     script:
-    def merge_list
     """
     for FILE in *final_polished_cds.fa
     do
@@ -138,14 +140,13 @@ process cds_align {
         tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam
         tuple val(sample_id), path("${sample_id}_read_aln_stats.tsv"), emit: stats
     script:
-    def ab = "${sample_id}_reads_aln_sorted.bam"
     """
     minimap2 -t ${params.threads} \
         -ax splice ${params.minimap2_opts} $polished_cds $sorted_reads_dir/sorted_reads.fastq |\
         samtools view -b - |\
-		samtools sort -o $ab;
+    samtools sort -o "${sample_id}_reads_aln_sorted.bam";
-		samtools index $ab;
+    samtools index "${sample_id}_reads_aln_sorted.bam";
-		((seqkit bam -s -j ${params.threads} ${ab} 2>&1)  | tee ${sample_id}_read_aln_stats.tsv ) || true
+    ((seqkit bam -s -j ${params.threads}  "${sample_id}_reads_aln_sorted.bam" 2>&1)  | tee ${sample_id}_read_aln_stats.tsv ) || true
 	"""
 }
 
@@ -168,8 +169,6 @@ process make_batches {
 
     minimum_batch_size = 2000
     """
-    nr_bases=\$(seqkit stats -T $fastq|cut -f 5| sed '2q;d')
-
     b=0
     if [ ${params.batch_size} -lt \$b ];
     then
@@ -186,11 +185,12 @@ process make_batches {
     echo "Batch size:\$batch_size";
     echo "Num bases: \$nr_bases";
 
-    init_cls_options="--batch-size \$batch_size --kmer-size ${params.kmer_size} \
+    init_cls_options="--batch-size \$batch_size --kmer-size ${params.kmer_size} --window-size ${params.window_size} \
-    --window-size ${params.window_size} --min-shared ${params.min_shared} --min-qual ${params.min_qual} \
+        --min-shared ${params.min_shared} --min-qual ${params.min_qual} \
         --mapped-threshold ${params.mapped_threshold} --aligned-threshold ${params.aligned_threshold} \
         --min-fraction ${params.min_fraction} --min-prob-no-hits ${params.min_prob_no_hits} \
-                          -M ${params.batch_max_seq} -P ${params.consensus_period} -g ${params.consensus_minimum} -c ${params.consensus_maximum} -F ${params.min_left_cls} "
+        -M ${params.batch_max_seq} -P ${params.consensus_period} -g ${params.consensus_minimum} \
+        -c ${params.consensus_maximum} -F ${params.min_left_cls} "
     mkdir -p sorted; isONclust2 sort \$init_cls_options -v -o sorted $fastq;
     """
 
@@ -206,7 +206,6 @@ process clustering() {
     script:
     """
     run_isonclust2.py $sorted_batches
-
     """
 }
 
@@ -216,8 +215,8 @@ process cluster_quality() {
     label "isoforms"
 
     input:
-        tuple val(sample_id), path(reads_fl), path(final_clusters_dir)
         path reference
+        tuple val(sample_id), path(reads_fl), path(final_clusters_dir)
     output:
         tuple val(sample_id),
             path("${sample_id}_cluster_qc"), emit: cluster_qc_dir
@@ -229,7 +228,6 @@ process cluster_quality() {
     def bam = "${qc_dir}/ref_aln.bam"
 
     """
-    echo $qc_dir_raw
     mkdir $qc_dir
     mkdir $qc_dir_raw
     minimap2 -ax splice -t 2 $reference $reads_fl |\
@@ -248,37 +246,35 @@ workflow denovo_assembly {
        fastq_reads_fl
        reference
     main:
-
         make_batches(fastq_reads_fl)
 
-        clustering(make_batches.output.sorted_batches)
+        clustering(make_batches.out.sorted_batches)
 
-        dump_clusters(clustering.output.root_cluster
+        dump_clusters(
-                      .join(make_batches.output.sorted_reads_dir))
+            clustering.out.root_cluster
+            .join(make_batches.out.sorted_reads_dir))
 
-        build_backbones(dump_clusters.output.final_clusters
+        build_backbones(
+            dump_clusters.out.final_clusters
             .flatMap(map_sample_ids_cls)
-                                        .groupTuple(size: 10, remainder: true)
+            .groupTuple(size: 10, remainder: true))
-                        )
 
-        merge_cds(build_backbones.output.polished_cds
+        merge_cds(
+            build_backbones.out.polished_cds
             .flatMap(map_sample_ids_cls)
-                                            .groupTuple()
+            .groupTuple())
-                  )
 
-        cds_align(merge_cds.out.final_polished_cds.view()
+        cds_align(
-                  .join(make_batches.output.sorted_reads_dir))
+            merge_cds.out.final_polished_cds
+            .join(make_batches.out.sorted_reads_dir))
 
         if (!reference.name.startsWith('OPTIONAL_FILE')){
+            cluster_quality(reference, fastq_reads_fl
+            .join(dump_clusters.out.final_clusters_dir))
 
-            cluster_quality(fastq_reads_fl
+            cluster_quality.out.cluster_qc_dir.set { opt_qual_ch }
-                .join(dump_clusters.output.final_clusters_dir), reference)
 
-            cluster_quality.out.cluster_qc_dir
+            cluster_quality.out.cluster_qc_raw.set { opt_qual_raw_ch }
-                .set { opt_qual_ch }
-
-            cluster_quality.output.cluster_qc_raw
-                .set { opt_qual_raw_ch }
 
         } else{
             Channel.empty().set { opt_qual_ch }
@@ -286,9 +282,9 @@ workflow denovo_assembly {
         }
 
     emit:
-       bam = cds_align.output.bam
+       bam = cds_align.out.bam
        cds = merge_cds.out.final_polished_cds
-       stats = cds_align.output.stats
+       stats = cds_align.out.stats
        opt_qual_ch
        opt_qual_raw_ch
 }

environment.yaml

@@ -1,4 +1,4 @@
-name: epi2melabs-wf-isoforms-test
+name: epi2melabs-wf-isoforms
 channels:
     - epi2melabs
     - bioconda
@@ -26,3 +26,4 @@ dependencies:
     - isonclust2==2.3
     - parallel
     - scikit-learn==1.0.2
+    - natsort

main.nf

@@ -42,6 +42,7 @@ process getVersions {
     python -c "import pysam; print(f'pysam,{pysam.__version__}')" >> versions.txt
     python -c "import aplanat; print(f'aplanat,{aplanat.__version__}')" >> versions.txt
     python -c "import pandas; print(f'pandas,{pandas.__version__}')" >> versions.txt
+    python -c "import sklearn; print(f'scikit-learn,{sklearn.__version__}')" >> versions.txt
     fastcat --version | sed 's/^/fastcat,/' >> versions.txt
     minimap2 --version | sed 's/^/minimap2,/' >> versions.txt
     samtools --version | head -n 1 | sed 's/ /,/' >> versions.txt
@@ -158,16 +159,17 @@ process assemble_transcripts{
     script:
         def out_filename = bam.name.replaceFirst(~/\.[^\.]+$/, '') + "_${sample_id}.gff"
         def G_FLAG = ref_annotation.name.startsWith('OPTIONAL_FILE') ? '' : "-G ${ref_annotation}"
+        def prefix = bam.name.split('-')[0][5..-1]
     """
     stringtie --rf ${G_FLAG} -L -v -A gene_abund.tab -p ${params.threads} ${params.stringtie_opts} -o  ${out_filename} \
-    ${bam} 2>/dev/null
+        -l $prefix ${bam} 2>/dev/null
      """
 }
 
 
 process merge_gff_bundles{
     /*
-    Merge gff bundles into a single gff file.
+    Merge gff bundles into a single gff file per sample.
     */
     label 'isoforms'
 
@@ -231,7 +233,6 @@ process run_gffcompare{
 process get_transcriptome{
         /*
         Write out a transcriptome file based on the query gff annotations.
-        TODO: Do we need to touch merged_transcriptome or can we just pass it?
         */
         label 'isoforms'
 
@@ -244,13 +245,10 @@ process get_transcriptome{
         def transcriptome = "${sample_id}_transcriptome.fas"
         def merged_transcriptome = "${sample_id}_merged_transcriptome.fas"
         """
-
         gffread -g ${reference_seq} -w ${transcriptome} ${transcripts_gff}
         if  [ "\$(ls -A $gffcmp_dir)" ];
             then
-            echo "Yes"
+                gffread -F -g ${reference_seq} -w ${merged_transcriptome} $gffcmp_dir/str_merged.annotated.gtf
-            gffread -F -g ${reference_seq} -w ${merged_transcriptome} \
-            $gffcmp_dir/str_merged.annotated.gtf
         fi
         """
 }
@@ -274,7 +272,6 @@ process makeReport {
     script:
         // Convert the sample_id arrayList.
         sids = new BlankSeparatedList(sample_ids)
-
         def report_name = "wf-isoforms-report.html"
         def OPT_ALN = denovo ?  '' : "--alignment_stats ${aln_stats}"
         def OPT_DENOVO = denovo ? "--denovo" : ''
@@ -365,15 +362,16 @@ workflow pipeline {
         run_gffcompare(merge_gff_bundles.out.gff, ref_annotation)
 
         if (params.denovo){
-            // Use the perd-sample, de novo-assembled CDS
+            // Use the per-sample, de novo-assembled CDS
             seq_for_transcriptome_build = m.cds
         }else {
-            // If doing reference based assembly, there is only one reference
+            // For reference based assembly, there is only one reference
             // So map this reference to all sample_ids
             seq_for_transcriptome_build = sample_ids.flatten().combine(Channel.fromPath(params.ref_genome))
         }
 
-        makeReport(software_versions,
+        makeReport(
+            software_versions,
             workflow_params,
             params.denovo,
             summariseConcatReads.out.summary
@@ -381,18 +379,19 @@ workflow pipeline {
             .join(run_gffcompare.out.gffcmp_dir)
             .join(preprocess_reads.out.report)
             .join(merge_gff_bundles.out.gff)
-                    .toList().transpose().toList()
+            .toList().transpose().toList())
-                    )
 
         report = makeReport.out.report
 
-       get_transcriptome(merge_gff_bundles.out.gff
+        get_transcriptome(
+            merge_gff_bundles.out.gff
             .join(run_gffcompare.out.gffcmp_dir)
             .join(seq_for_transcriptome_build))
 
        if (use_ref_ann){
             results = preprocess_reads.out.report
-                        .concat(run_gffcompare.output.gffcmp_dir,
+                      .concat(
+                          run_gffcompare.output.gffcmp_dir,
                           m.stats,
                           get_transcriptome.out.flatMap(map_sample_ids_cls))
                       .map {it -> it[1]}

reference_assembly.nf

@@ -16,8 +16,8 @@ process map_reads{
        tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam
        tuple val(sample_id), path("${sample_id}_read_aln_stats.tsv"), emit: stats
     script:
-        def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context}, RightShift: ${params.poly_context},
+        def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context},
-           RegexEnd: "[Aa]{${params.max_poly_run},}",
+        RightShift: ${params.poly_context}, RegexEnd: "[Aa]{${params.max_poly_run},}",
         Stranded: True,Invert: True, Tsv: "internal_priming_fail.tsv"} """
     """
     seqkit fq2fa ${fastq_reads} -o "reads.fa";