From 2c4b4fe79cf78a3e3ed96b6fa3d8f84509079df0 Mon Sep 17 00:00:00 2001 From: Kivanc Corut Date: Sun, 16 Apr 2023 17:48:32 -0400 Subject: [PATCH] feat: add rules for generating igv reports (igv-report) * add a rule for generating IGV reports from source read mapping results * add a rule for generating IGV reports from contig mapping results --- .../report/generate_igv_report_contigs.rst | 1 + workflow/report/generate_igv_report_reads.rst | 4 + workflow/rules/igv_report.smk | 95 +++++++++++++++++++ 3 files changed, 100 insertions(+) create mode 100644 workflow/report/generate_igv_report_contigs.rst create mode 100644 workflow/report/generate_igv_report_reads.rst create mode 100644 workflow/rules/igv_report.smk diff --git a/workflow/report/generate_igv_report_contigs.rst b/workflow/report/generate_igv_report_contigs.rst new file mode 100644 index 0000000..5b9aa45 --- /dev/null +++ b/workflow/report/generate_igv_report_contigs.rst @@ -0,0 +1 @@ +IGV report of the mapping results of assembled contigs. Source reads of trait-associated k-mers were assembled using `Spades `_. Resulted contigs were mapped to the reference genome using `minimap2 `_. Mapping results were filtered based minimum mapping quality score of {{ snakemake.config["params"]["filter_alignment"]["min_map_score"] }}. Alignment BAM files were converted to BED format using `bedtools `_. IGV HTML report was generated using `igv-reports `_. IGV report consists of alignment BAM and BED files, the reference genome FASTA and the reference gene annotations in GTF/GFF3. \ No newline at end of file diff --git a/workflow/report/generate_igv_report_reads.rst b/workflow/report/generate_igv_report_reads.rst new file mode 100644 index 0000000..4f773d1 --- /dev/null +++ b/workflow/report/generate_igv_report_reads.rst @@ -0,0 +1,4 @@ +IGV report of the mapping results of source reads of associated k-mers. Source reads of trait-associated k-mers were mapped to the reference genome using bowtie2_. Mapping results were filtered based minimum mapping quality score of {{ snakemake.config["params"]["filter_alignment"]["min_map_score"] }}. Alignment BAM files were converted to BED format using `bedtools `_. BED files then merged using `bedtools merge `_ with the parameters `-s -c 6 -o distinct`. IGV HTML report was generated using `igv-reports `_. IGV report consists of alignment BAM and BED files, the reference genome FASTA and the reference gene annotations in GTF/GFF3. + +.. _bowtie: https://bowtie-bio.sourceforge.net/index.shtml +.. _bowtie2: https://bowtie-bio.sourceforge.net/bowtie2/index.shtml \ No newline at end of file diff --git a/workflow/rules/igv_report.smk b/workflow/rules/igv_report.smk new file mode 100644 index 0000000..e111a87 --- /dev/null +++ b/workflow/rules/igv_report.smk @@ -0,0 +1,95 @@ +# ================================================================================================= +# Generate IGV HTML report for contigs alignment +# ================================================================================================= + +rule align_contigs_igv_report: + input: + fasta="resources/ref/genome/genome.fasta", + fai="resources/ref/genome/genome.fasta.fai", + vcf="results/align_contigs/{phenos_filt}/alignment/{phenos_filt}_contigs_aligned.filter.sorted.bed", + tracks= [ + get_genome_annotation(), + "results/align_contigs/{phenos_filt}/alignment/{phenos_filt}_contigs_aligned.filter.sorted.bam", + ], + output: + igv_report= report( + "results/igv_reports/align_contigs/{phenos_filt}/{phenos_filt}_contigs_aligned.igv-report.html", + caption= "../report/generate_igv_report_contigs.rst", + category= "IGV Reports - Contig Mapping", + subcategory="{phenos_filt}", + ) + params: + extra= config["params"]["igv_report"]["extra"] # optional params + log: + "logs/igv_report/align_contigs/{phenos_filt}.igv-report.log" + wrapper: + "v1.25.0/bio/igv-reports" + +# ========================================================================================================= +# Aggregate igv_report outputs +# ========================================================================================================= + +rule aggregate_align_contigs_igv_reports: + input: + aggregate_input_align_contigs_igv_report + output: + "results/igv_reports/align_contigs.igv_report.done" + log: + "logs/igv_report/aggregate_align_contigs_igv_reports.log" + message: + "Checking if generating igv_report is done..." + shell: + """ + touch {output} 2> {log} + """ + +# ================================================================================================= +# Generate IGV HTML report for reads alignments +# ================================================================================================= + +rule align_reads_igv_report: + input: + fasta="resources/ref/genome/genome.fasta", + fai="resources/ref/genome/genome.fasta.fai", + vcf="results/align_reads_with_kmers/{phenos_filt}/{phenos_filt}.align_reads_with_kmers.filter.sorted.merged.bed", + tbi = "results/align_reads_with_kmers/{phenos_filt}/{phenos_filt}.align_reads_with_kmers.filter.sorted.merged.bed.gz.tbi", + tracks= [ + get_genome_annotation(), + "results/align_reads_with_kmers/{phenos_filt}/{phenos_filt}.align_reads_with_kmers.filter.sorted.bam", + ], + done = "results/align_reads_with_kmers/{phenos_filt}/bam_to_bed.done" + output: + igv_report= report( + "results/igv_reports/align_reads/{phenos_filt}/{phenos_filt}_reads_aligned.igv-report.html", + caption= "../report/generate_igv_report_reads.rst", + category= "IGV Reports - Read Mapping", + subcategory="{phenos_filt}", + ) + params: + extra= config["params"]["igv_report"]["extra"] # optional params + threads: + config["params"]["igv_report"]["threads"] + log: + "logs/igv_report/align_reads/{phenos_filt}.igv-report.log" + wrapper: + "v1.25.0/bio/igv-reports" + +# ========================================================================================================= +# Aggregate igv_report outputs +# ========================================================================================================= + +rule aggregate_align_reads_igv_reports: + input: + aggregate_input_align_reads_igv_report + output: + "results/igv_reports/align_reads.igv_report.done" + log: + "logs/igv_report/aggregate_align_reads_igv_reports.log" + message: + "Checking if generating igv_report is done..." + shell: + """ + touch {output} 2> {log} + """ + +# ========================================================================================================= \ No newline at end of file