wf-transcriptomes-v202/subworkflows/reference_assembly.nf

61 lines
2.4 KiB
Plaintext

process map_reads{
/*
Map reads to reference using minimap2.
Filter reads by mapping quality.
Filter internally-primed reads.
*/
label "isoforms"
cpus params.threads
memory "31 GB"
publishDir path: "${params.out_dir}/${publish_prefix_bams}", mode: 'copy', pattern: "${sample_id}_reads_aln_sorted.bam*", overwrite: true
input:
tuple val(sample_id), path (fastq_reads), path(index), path(reference)
val publish_prefix_bams
output:
tuple val(sample_id),
path("${sample_id}_reads_aln_sorted.bam"),
path("${sample_id}_reads_aln_sorted.bam.bai"),
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"} """
def mm2_threads = Math.max(task.cpus - 3, 1)
"""
minimap2 -t ${mm2_threads} -ax splice ${params.minimap2_opts} ${index} ${fastq_reads}\
| samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
| seqkit bam -j 1 -x -T '${ContextFilter}' -\
| samtools sort --write-index -@ 1 -o "${sample_id}_reads_aln_sorted.bam##idx##${sample_id}_reads_aln_sorted.bam.bai" - ;
((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j 1 - 2>&1) | tee "${sample_id}_read_aln_stats.tsv" ) || true
# Add sample id header and column; remove last column (File)
cat "${sample_id}_read_aln_stats.tsv" \
| sed "s/^/${sample_id} /" \
| sed "1 s/^${sample_id}/sample_id/" \
| awk 'NF{NF-=1};1' \
> tmp
mv tmp "${sample_id}_read_aln_stats.tsv"
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 {
take:
index
reference
fastq_reads
publish_prefix_bams
main:
map_reads(fastq_reads.combine(index).combine(reference), publish_prefix_bams)
emit:
bam = map_reads.out.bam
stats = map_reads.out.stats
}