diff --git a/README.md b/README.md
index 7347de7..d8cb0ec 100644
--- a/README.md
+++ b/README.md
@@ -115,13 +115,12 @@ For example:
* wf-isoforms-report.html
- 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: Currently only available for the reference-guilded assembly
-The report includes a searchable filterable table of transcript isoforms
+
+Note transcript isoform table is currently only available for the reference-guilded assembly.
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,11 +132,8 @@ 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
* [nextflow](https://www.nextflow.io/)
diff --git a/bin/merge_gff.py b/bin/merge_gff.py
new file mode 100755
index 0000000..5dd56cd
--- /dev/null
+++ b/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)
diff --git a/bin/report.py b/bin/report.py
index a15431d..c570dee 100755
--- a/bin/report.py
+++ b/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.
- 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)
- if num_bins > 12:
- num_bins = 12
-
- 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 = hist.histogram(
+ [isoforms_per_gene], colors=[Colors.cerulean],
+ title="isoforms per gene")
bar_isos.xaxis.major_label_orientation = math.pi / 2.8
plots.append(bar_isos)
diff --git a/denovo_assembly.nf b/denovo_assembly.nf
index 03abae1..e813e87 100644
--- a/denovo_assembly.nf
+++ b/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
- 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
+ Build coding
*/
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)
"""
@@ -56,64 +53,70 @@ process build_backbones {
UNID=\$(echo ${cluster_fq_bl[1]} | grep -o -E '[0-9]+')
for cluster in $cluster_fq_bl
- do
- clfq=`basename \$cluster`
- cln=\${clfq%.*}
+ do
+ clfq=`basename \$cluster`
+ cln=\${clfq%.*}
- echo Building backbone for cluster: \$cln
- echo "\tSampling input reads for backbone construction."
- sample=\${cln}_sample.fq
- seqkit head --quiet -n 100 \$cluster > \$sample
- seqkit sample --quiet -n 500 -2 -s 100 \$cluster >> \$sample
- echo "\tConstructing spoa consensus."
- spoa_cons=\${cln}_spoa.fa
- spoa -m 5 -n -4 -g -8 -e -6 -q -10 -c -15 -l 1 -r 0 \$sample > \$spoa_cons
+ echo Building backbone for cluster: \$cln
+ echo "\tSampling input reads for backbone construction."
+ sample=\${cln}_sample.fq
+ seqkit head --quiet -n 100 \$cluster > \$sample
+ seqkit sample --quiet -n 500 -2 -s 100 \$cluster >> \$sample
+ echo "\tConstructing spoa consensus."
+ spoa_cons=\${cln}_spoa.fa
+ spoa -m 5 -n -4 -g -8 -e -6 -q -10 -c -15 -l 1 -r 0 \$sample > \$spoa_cons
- echo "\tPolishing the consensus using racon."
+ echo "\tPolishing the consensus using racon."
- # polish 1
- samgz=\${cln}_aln.sam.gz
- racon_cons=\${cln}_racon.fa
- racon_cons1=\${cln}_racon1.fa
+ # polish 1
+ samgz=\${cln}_aln.sam.gz
+ racon_cons=\${cln}_racon.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
- if [ -f \$racon_cons ];
- then
- cat \$racon_cons | seqkit replace -p ".*" -r cluster_\${cln} > \$racon_cons1
- else
- continue
- fi
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$spoa_cons \$sample | gzip - > \$samgz
+ exitcode=0
+ racon -t ${params.threads} --no-trimming -u -w 2000 \$sample \$samgz \$spoa_cons > \$racon_cons || exitcode=1
+ if [[ \$exitcode -eq 0 ]];
+ then
+ # 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
- racon_cons2=\${cln}_racon2.fa
- minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons1 \$sample | gzip - > \$samgz1
- racon -t ${params.threads} -u --no-trimming \$sample \$samgz1 \$racon_cons1 > \$racon_cons2 || true
- if [ -f \$racon_cons2 ];
- then
- echo "success polish 2"
- else
- continue
- fi
+ # polish 2
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons \$sample | gzip - > \$samgz
+ exitcode=0
+ racon -t ${params.threads} -u --no-trimming \$sample \$samgz \$racon_cons > \$tmpcons || exitcode=1
+ if [[ \$exitcode -eq 0 ]];
+ 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_cons2 \$sample | gzip - > \$samgz2
- EXITCODE=0
- racon -t ${params.threads} -u \$sample \$samgz2 \$racon_cons2 >> ${sample_id}_\${UNID}_final_polished_cds.fa || EXITCODE=\$?
- if [ \$EXITCODE -eq 0 ];
- then
- echo "finished"
- fi
+ # polish 3
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons \$sample | gzip - > \$samgz
+ exitcode=0
+ racon -t ${params.threads} -u \$sample \$samgz \$racon_cons > \$tmpcons || exitcode=1
+ if [[ \$exitcode -eq 0 ]];
+ then
+ cat \$tmpcons >> ${sample_id}_\${UNID}_final_polished_cds.fa
+ echo "polishing 3 success"
+ else
+ echo "Polishing failed for \${cln}"
+ fi
- done
+ 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 index $ab;
- ((seqkit bam -s -j ${params.threads} ${ab} 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
+ minimap2 -t ${params.threads} \
+ -ax splice ${params.minimap2_opts} $polished_cds $sorted_reads_dir/sorted_reads.fastq |\
+ samtools view -b - |\
+ samtools sort -o "${sample_id}_reads_aln_sorted.bam";
+ samtools index "${sample_id}_reads_aln_sorted.bam";
+ ((seqkit bam -s -j ${params.threads} "${sample_id}_reads_aln_sorted.bam" 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
"""
}
@@ -168,16 +169,14 @@ 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
nr_bases=\$(seqkit stats -T $fastq|cut -f 5| sed '2q;d')
let batch_size=\$nr_bases/1000/$maxcpus
if [ \$batch_size -lt $minimum_batch_size ];
- then
- batch_size=$minimum_batch_size
+ then
+ batch_size=$minimum_batch_size
fi
else
batch_size=${params.batch_size}
@@ -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} \
- --window-size ${params.window_size} --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} "
+ init_cls_options="--batch-size \$batch_size --kmer-size ${params.kmer_size} --window-size ${params.window_size} \
+ --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} "
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,16 +228,15 @@ 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 |\
- samtools view -q 2 -F 2304 -b - |\
- samtools sort - -o $bam;
+ samtools view -q 2 -F 2304 -b - |\
+ samtools sort - -o $bam;
samtools index $bam;
compute_cluster_quality.py --sizes $final_clusters_dir/clusters_info.tsv \
- --outfile ${qc_dir}/cluster_quality.csv --ont --clusters $final_clusters_dir/clusters.tsv \
- --classes $bam --report ${qc_dir}/cluster_quality.pdf --raw_data_out $qc_dir_raw
+ --outfile ${qc_dir}/cluster_quality.csv --ont --clusters $final_clusters_dir/clusters.tsv \
+ --classes $bam --report ${qc_dir}/cluster_quality.pdf --raw_data_out $qc_dir_raw
"""
}
@@ -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
- .join(make_batches.output.sorted_reads_dir))
+ dump_clusters(
+ clustering.out.root_cluster
+ .join(make_batches.out.sorted_reads_dir))
- build_backbones(dump_clusters.output.final_clusters
- .flatMap(map_sample_ids_cls)
- .groupTuple(size: 10, remainder: true)
- )
+ build_backbones(
+ dump_clusters.out.final_clusters
+ .flatMap(map_sample_ids_cls)
+ .groupTuple(size: 10, remainder: true))
- merge_cds(build_backbones.output.polished_cds
- .flatMap(map_sample_ids_cls)
- .groupTuple()
- )
+ merge_cds(
+ build_backbones.out.polished_cds
+ .flatMap(map_sample_ids_cls)
+ .groupTuple())
- cds_align(merge_cds.out.final_polished_cds.view()
- .join(make_batches.output.sorted_reads_dir))
+ cds_align(
+ 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
- .join(dump_clusters.output.final_clusters_dir), reference)
+ cluster_quality.out.cluster_qc_dir.set { opt_qual_ch }
- cluster_quality.out.cluster_qc_dir
- .set { opt_qual_ch }
-
- cluster_quality.output.cluster_qc_raw
- .set { opt_qual_raw_ch }
+ cluster_quality.out.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
}
diff --git a/environment.yaml b/environment.yaml
index 4921c24..fd72183 100644
--- a/environment.yaml
+++ b/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
diff --git a/main.nf b/main.nf
index e0d2329..a81fcaf 100644
--- a/main.nf
+++ b/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
@@ -105,7 +106,7 @@ process build_minimap_index{
path "genome_index.mmi", emit: index
script:
"""
- minimap2 -t ${params.threads} ${params.minimap_index_opts} -I 1000G -d "genome_index.mmi" ${reference}
+ minimap2 -t ${params.threads} ${params.minimap_index_opts} -I 1000G -d "genome_index.mmi" ${reference}
"""
}
@@ -128,7 +129,7 @@ process split_bam{
if (params["bundle_min_reads"] != false)
"""
seqkit bam -j ${params.threads} -N ${params.bundle_min_reads} ${bam} -o bam_bundles/
- mv bam_bundles/* .
+ mv bam_bundles/* .
for f in *:*; do mv -v "\$f" \$(echo "\$f" | tr ':' '-'); done
"""
else
@@ -145,7 +146,7 @@ process assemble_transcripts{
Take aligned reads in bam format that may be a chunk of a larger alignment file.
Optionally use reference annotation to guide assembly.
- Output gff annotation files in a tuple with `sample_id` for combining into samples late rin the pipeline.
+ Output gff annotation files in a tuple with `sample_id` for combining into samples later in the pipeline.
*/
label 'isoforms'
cpus params.threads
@@ -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'
@@ -215,14 +217,14 @@ process run_gffcompare{
mkdir $out_dir
echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation"
- gffcompare -o ${out_dir}/str_merged -r ${ref_annotation} \
- ${params.gffcompare_opts} ${query_annotation}
+ gffcompare -o ${out_dir}/str_merged -r ${ref_annotation} \
+ ${params.gffcompare_opts} ${query_annotation}
generate_tracking_summary.py --tracking $out_dir/str_merged.tracking \
- --output_dir ${out_dir} --annotation ${ref_annotation}
+ --output_dir ${out_dir} --annotation ${ref_annotation}
- mv *.tmap $out_dir
- mv *.refmap $out_dir
+ mv *.tmap $out_dir
+ mv *.refmap $out_dir
"""
}
}
@@ -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
+ then
+ 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,61 +362,63 @@ 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,
- workflow_params,
- params.denovo,
- summariseConcatReads.out.summary
- .join(m.stats)
- .join(run_gffcompare.out.gffcmp_dir)
- .join(preprocess_reads.out.report)
- .join(merge_gff_bundles.out.gff)
- .toList().transpose().toList()
- )
+ makeReport(
+ software_versions,
+ workflow_params,
+ params.denovo,
+ summariseConcatReads.out.summary
+ .join(m.stats)
+ .join(run_gffcompare.out.gffcmp_dir)
+ .join(preprocess_reads.out.report)
+ .join(merge_gff_bundles.out.gff)
+ .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,
- m.stats,
- get_transcriptome.out.flatMap(map_sample_ids_cls))
- .map {it -> it[1]}
- .concat(makeReport.out.report)
+ .concat(
+ run_gffcompare.output.gffcmp_dir,
+ m.stats,
+ get_transcriptome.out.flatMap(map_sample_ids_cls))
+ .map {it -> it[1]}
+ .concat(makeReport.out.report)
}
if (!use_ref_ann && !params.denovo){
results = preprocess_reads.out.report
- .concat(m.stats,
- get_transcriptome.out.flatMap(map_sample_ids_cls))
- .map {it -> it[1]}
- .concat(makeReport.out.report)
+ .concat(m.stats,
+ get_transcriptome.out.flatMap(map_sample_ids_cls))
+ .map {it -> it[1]}
+ .concat(makeReport.out.report)
}
if (params.denovo){
results = m.cds
.concat(m.stats,
- seq_for_transcriptome_build,
- get_transcriptome.out.flatMap(map_sample_ids_cls),
- merge_gff_bundles.out.gff,
- m.opt_qual_ch.flatMap {it ->
- l = []
- for (x in it[1..-1]){
- l.add(tuple(it[0], x))
- }
- return l
- })
+ seq_for_transcriptome_build,
+ get_transcriptome.out.flatMap(map_sample_ids_cls),
+ merge_gff_bundles.out.gff,
+ m.opt_qual_ch.flatMap {it ->
+ l = []
+ for (x in it[1..-1]){
+ l.add(tuple(it[0], x))
+ }
+ return l
+ })
.map {it -> it[1]}
.concat(makeReport.out.report)
}
diff --git a/reference_assembly.nf b/reference_assembly.nf
index d88767b..ae58d79 100644
--- a/reference_assembly.nf
+++ b/reference_assembly.nf
@@ -16,23 +16,23 @@ 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},
- RegexEnd: "[Aa]{${params.max_poly_run},}",
- Stranded: True,Invert: True, Tsv: "internal_priming_fail.tsv"} """
- """
- seqkit fq2fa ${fastq_reads} -o "reads.fa";
- minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} "reads.fa"\
- | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
- | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\
- | samtools sort -@ ${params.threads} -o "${sample_id}_reads_aln_sorted.bam" - ;
- ((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j ${params.threads} - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
+ def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context},
+ 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";
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} "reads.fa"\
+ | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
+ | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\
+ | samtools sort -@ ${params.threads} -o "${sample_id}_reads_aln_sorted.bam" - ;
+ ((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j ${params.threads} - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
- if [[ -s "internal_priming_fail.tsv" ]];
- then
- tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$4 }}' - > "context_internal_priming_fail_start.fasta"
- tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$6 }}' - > "context_internal_priming_fail_end.fasta"
- fi
- """
+ if [[ -s "internal_priming_fail.tsv" ]];
+ then
+ tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$4 }}' - > "context_internal_priming_fail_start.fasta"
+ tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$6 }}' - > "context_internal_priming_fail_end.fasta"
+ fi
+ """
}
workflow reference_assembly {