Merge branch 'post-cw460' into 'dev'

Post cw460

See merge request epi2melabs/workflow-containers/wf-isoforms!43

bin/merge_gff.py

@@ -1,54 +0,0 @@
-#!/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

@@ -365,6 +365,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
     Features present in the query transcripts, but absent in the reference
     ''')
 
+    tabs = []
     for id_, dir_ in zip(sample_ids, gffcompare_outdirs):
         stats, _, miss, novel, total = \
             parse_gffcmp_stats(dir_ / 'str_merged.stats')
@@ -375,12 +376,13 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
         bar_missed = grouped_bar(miss, title="Missed")
         bar_novel = grouped_bar(novel, title="Novel")
 
-        grid = gridplot([bar_totals, bar_performance, bar_missed, bar_novel],
+        tabs.append(Panel(
-                        ncols=4, width=270, height=280)
+            child=gridplot(
-        section.markdown("""
+                [bar_totals, bar_performance, bar_missed, bar_novel],
-        #### Sample_id: {}
+                ncols=2, width=350, height=260), title=id_))
-        """.format(id_))
+
-        section.plot(grid)
+    cover_panel = Tabs(tabs=tabs)
+    section.plot(cover_panel)
 
     names = {
         '=': 'ExactMatch:=',
@@ -403,7 +405,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
     # Plot overlaps panel:
     section = report.add_section()
     section.markdown('''
-    ### Query transfrag class assignments
+    ### Query transfrag classes
 
     The classes that are assigned by
     [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml),
@@ -435,15 +437,17 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
             tracking['Overlaps'].values.tolist(),
             tracking['Percent'].values.tolist(), title="{}".format(id_))
 
-        tracking.drop(columns=['sample_id'], inplace=True)
         tracking_dfs.append(tracking)
 
-        tracking['description'] = pd.Series(tracking.Overlaps.apply(
+        tracking['Description'] = pd.Series(tracking.Overlaps.apply(
             lambda x: x.split(':')[0]))
 
-        tracking['code'] = pd.Series(tracking.Overlaps.apply(
+        tracking['Code'] = pd.Series(tracking.Overlaps.apply(
             lambda x: x.split(':')[1]))
 
+        tracking.drop(columns=['sample_id', 'Overlaps'], inplace=True)
+        tracking = tracking[['Code', 'Description', 'Count', 'Percent']]
+
         cols = [TableColumn(field=Ci, title=Ci, width=100)
                 for Ci in tracking.columns]
 

main.nf

@@ -159,7 +159,12 @@ 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]
+        // Convert batch name to stringtie prefix to prevent clashing attribute names
+        // eg "0000000123_cluster..." to 123
+        def prefix =  StringUtils.stripStart(bam.name.split('_')[0],"0")
+        if (!prefix){
+            prefix = "0"
+        }
     """
     stringtie --rf ${G_FLAG} -L -v -A gene_abund.tab -p ${params.threads} ${params.stringtie_opts} -o  ${out_filename} \
         -l $prefix ${bam} 2>/dev/null