diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..978bd3a --- /dev/null +++ b/.gitattributes @@ -0,0 +1,4 @@ +# Keep text byte-for-byte reproducible across Linux and Windows checkouts. +* text=auto eol=lf +*.png binary +*.gz binary diff --git a/.gitignore b/.gitignore index 3b93011..9892427 100644 --- a/.gitignore +++ b/.gitignore @@ -41,3 +41,7 @@ Step_03_UMIAdapterTrimR1_old_v1.py # Backup files *.bak* *.bak_* + +# Local configuration and sample results +EVscope.local.conf +*_EVscope_output/ diff --git a/EVscope.conf b/EVscope.conf index 1543901..98cf459 100644 --- a/EVscope.conf +++ b/EVscope.conf @@ -1,6 +1,6 @@ # ============================================================================== # CONFIGURATION FILE FOR THE EVSCOPE PIPELINE -# Version: 1.0.0 +# Version: 1.1.0 # ============================================================================== # --- A. Core Paths --- @@ -37,7 +37,60 @@ RSEM_BOWTIE2_INDEX="${EVscope_PATH}/references/index/RSEM_bowtie2_index_3659642_ STAR_INDEX="${EVscope_PATH}/references/index/Index_STAR-2.7.10b_sjdbOverhang99_GeneCode45" BWA_INDEX="${EVscope_PATH}/references/index/Index_BWA_V0.7.18/hg38.p14.whole.genome.fa" -# --- F. General Annotation --- +# --- F. STAR Alignment Parameters --- +# Bash arrays are used so every option/value remains a separate command-line token. +# Edit the appropriate profile to change an existing option. Add new, non-duplicate +# options to STAR_EXTRA_ARGS or the kit-specific EXTRA array. +STAR_COMMON_ARGS=( + --readFilesCommand zcat + --outFilterMultimapNmax 100 + --winAnchorMultimapNmax 100 + --outSAMtype BAM SortedByCoordinate +) +STAR_TOTAL_RNA_ARGS=( + --twopassMode Basic + --chimSegmentMin 10 + --chimJunctionOverhangMin 10 + --chimScoreMin 1 + --chimOutType Junctions WithinBAM +) +STAR_SMRNA_SE_ARGS=( + --twopassMode None + --alignEndsType EndToEnd + --alignIntronMax 1 + --outFilterMismatchNmax 1 + --outFilterMultimapScoreRange 0 + --outFilterScoreMinOverLread 0 + --outFilterMatchNminOverLread 0 + --outFilterScoreMin 12 + --outFilterMatchNmin 14 + --chimSegmentMin 10 + --chimJunctionOverhangMin 10 + --chimScoreMin 1 + --chimOutType Junctions WithinBAM +) +STAR_SMRNA_PE_ARGS=( + --twopassMode None + --alignEndsType EndToEnd + --alignIntronMax 1 + --outFilterMismatchNmax 1 + --outFilterMultimapScoreRange 0 + --outFilterScoreMinOverLread 0 + --outFilterMatchNminOverLread 0 + --outFilterScoreMin 24 + --outFilterMatchNmin 28 + --chimSegmentMin 10 + --chimJunctionOverhangMin 10 + --chimScoreMin 1 + --chimOutType Junctions WithinBAM +) +# Extra options applied to every STAR run, or only to the selected kit. +# Example: STAR_SMRNA_EXTRA_ARGS=(--limitBAMsortRAM 30000000000) +STAR_EXTRA_ARGS=() +STAR_TOTAL_RNA_EXTRA_ARGS=() +STAR_SMRNA_EXTRA_ARGS=() + +# --- G. General Annotation --- PIPELINE_BIN_DIR="${EVscope_PATH}/bin" TOTAL_GENE_GTF="${EVscope_PATH}/references/annotations_HG38/HG38_3659642_combined_RNAs_with_gold_standard_piRNAs.gtf" TOTAL_GENE_REFFLAT="${EVscope_PATH}/references/annotations_HG38/HG38_3659642_combined_RNAs_with_gold_standard_piRNAs.refflat" @@ -48,7 +101,7 @@ GENCODE_V45_non_overlapping_exon_BED="${EVscope_PATH}/references/annotations_HG3 HUMAN_RRNA_INTERVAL="${EVscope_PATH}/references/annotations_HG38/GeneCodeV45_combined_Nucl_Mt_rRNA_589.interval_list" ENCODE_BLACKLIST_BED="${EVscope_PATH}/references/annotations_HG38/Encode_hg38-blacklist.v2.bed" -# --- G. Step 18: Merged SAF for meta-gene region counting --- +# --- H. Step 18: Merged SAF for meta-gene region counting --- # Single merged SAF with 7 mutually exclusive regions (801,166 features). # GeneID format: {region}__region_{N} (e.g. exon__region_1, intron__region_42) # Regions: 5UTR, 3UTR, exon, intron, promoter, downstream, intergenic @@ -57,7 +110,7 @@ export STEP18_MERGED_SAF="${EVscope_PATH}/references/annotations_HG38/HG38_all7_ # Region labels for output (order matters for pie chart) export STEP18_REGION_LABELS="5UTR,3UTR,exon,intron,promoter,downstream,intergenic" -# --- H. Step 26a: BED files for RNA Type Density Plot --- +# --- I. Step 26a: BED files for RNA Type Density Plot --- export STEP26_RNATYPE_BEDS=( "${EVscope_PATH}/references/annotations_HG38/HG38_protein_coding.bed" "${EVscope_PATH}/references/annotations_HG38/HG38_tRNAs.bed" @@ -81,7 +134,7 @@ export STEP26_RNATYPE_BEDS=( ) export STEP26_RNATYPE_LABELS="[protein_coding,tRNAs,rRNAs,Y_RNAs,snoRNAs,snRNAs,scaRNAs,vault_RNAs,lncRNAs,TEC_protein_coding,IG_genes,TR_genes,misc-sncRNAs,miRNAs,pseudogenes,ERVs,LINEs,SINEs,piRNAs]" -# --- I. Step 26b: BED files for Meta-Gene Density Plot --- +# --- J. Step 26b: BED files for Meta-Gene Density Plot --- export STEP26_METAGENE_BEDS=( "${EVscope_PATH}/references/annotations_HG38/HG38_promoter_upstream2000bp_noOverlap.bed" "${EVscope_PATH}/references/annotations_HG38/HG38_5UTR_noOverlap.bed" @@ -93,7 +146,7 @@ export STEP26_METAGENE_BEDS=( ) export STEP26_METAGENE_LABELS="[Promoters,5'UTR,Exon,Intron,3'UTR,Downstream,Intergenic]" -# --- J. Deconvolution Refs --- +# --- K. Deconvolution Refs --- GTEX_SMTS_REF="${EVscope_PATH}/references/deconvolution_HG38/human_ref_GTEx_v10_SMTS_AVG_TPM.csv" Monaco2020_ImmuneCell_REF="${EVscope_PATH}/references/deconvolution_HG38/human_ref_Monaco2020_ImmuneCell_AVG_TPM.csv" BRAIN_SC_REF="${EVscope_PATH}/references/deconvolution_HG38/human_ref_brain_single_cell_atlasV1_31superclusters_AVP_CPM.csv" diff --git a/EVscope.sh b/EVscope.sh index 4dbf968..c4fac8b 100755 --- a/EVscope.sh +++ b/EVscope.sh @@ -2,7 +2,7 @@ # SPDX-License-Identifier: MIT # ============================================================================== # EVscope.sh: Modular RNA-seq Analysis Pipeline -# Version: 1.0.0 +# Version: 1.1.0 # Description: # RNA-seq pipeline with comprehensive error handling, strict # variable scoping, robust input validation, and production-grade QC. Features @@ -23,6 +23,7 @@ # Repository: https://github.com/TheDongLab/EVscope # Changelog: # v1.0.0 - Initial release with 27-step modular workflow +# v1.1.0 - Add kit-aware SMARTer smRNA-seq preprocessing and no-UMI mode # ============================================================================== # ------------------------------------------------------------------------------ # SHELL OPTIONS AND SAFETY SETTINGS @@ -42,7 +43,7 @@ if command -v realpath &>/dev/null; then else readonly SCRIPT_DIR="$(cd -P "$(dirname "${BASH_SOURCE[0]}")" && pwd -P)" fi -readonly VERSION="1.0.0" +readonly VERSION="1.1.0" readonly PIPELINE_BASE_DIR="$SCRIPT_DIR" readonly MIN_BASH_VERSION="4.0" readonly MAX_PIPELINE_STEPS=27 @@ -59,6 +60,9 @@ declare -- circ_tool="both" declare -- read_count_mode="uniq" declare -- gDNA_correction="no" declare -- strand="unstrand" +declare -- kit="smarter_total_rna" +declare -- alignment_tag="umi_dedup" +declare -- result_dedup_tag="dedup" declare -- config_file="" declare -- output_dir="" declare -- fastq_read1="" @@ -155,15 +159,16 @@ ${C_BOLD}USAGE:${C_RESET} bash ${SCRIPT_NAME} [OPTIONS] --sample_name --input_fastqs [file2] ${C_BOLD}REQUIRED ARGUMENTS:${C_RESET} --sample_name Unique sample identifier (alphanumeric, underscore, hyphen, dot) - --input_fastqs Input FASTQ file(s): one for SE, two for PE (space- or comma-separated) + --input_fastqs Input FASTQ(s): Total RNA and smRNA accept SE or paired R1/R2 ${C_BOLD}OPTIONAL ARGUMENTS:${C_RESET} --threads CPU threads for parallel operations (default: 1) + --kit Library kit: smarter_total_rna (default) or smarter_smrna --run_steps Steps to execute: "all" or comma-separated list/ranges --skip_steps Steps to skip (same format as --run_steps) --circ_tool circRNA detection tool (default: both) --read_count_mode Read counting strategy (default: uniq) --gDNA_correction Apply genomic DNA correction (default: no) - --strand Library strand orientation (default: unstrand) + --strand Library strand orientation (default: unstrand; smRNA: forward) --config Path to configuration file --resume Resume an existing output directory after metadata checks --force Re-run steps and allow overwriting an existing output directory @@ -198,6 +203,22 @@ sanitize_string() { local input="${1:-}" echo "${input//[^a-zA-Z0-9_.-]/_}" } +normalize_kit() { + local value="${1:-}" + value="${value,,}" + value="${value//-/_}" + case "$value" in + smarter_total_rna|smartseq_total_rna|smartseq_total|total_rna) + printf '%s\n' "smarter_total_rna" + ;; + smarter_smrna|smarter_smrna_seq|smrna|smrna_seq) + printf '%s\n' "smarter_smrna" + ;; + *) + return 1 + ;; + esac +} get_absolute_path() { local path="${1:-}" if [[ -z "$path" ]]; then @@ -334,7 +355,10 @@ check_dependencies() { fi any_step_selected 1 4 && add_unique_item core_dependencies "fastqc" - step_is_selected 3 && { add_unique_item core_dependencies "umi_tools"; add_unique_item core_dependencies "cutadapt"; } + step_is_selected 3 && add_unique_item core_dependencies "cutadapt" + if [[ "$kit" == "smarter_total_rna" ]] && any_step_selected 3 6; then + add_unique_item core_dependencies "umi_tools" + fi step_is_selected 6 && { add_unique_item core_dependencies "STAR"; add_unique_item core_dependencies "samtools"; } any_step_selected 7 8 11 12 13 18 25 && add_unique_item core_dependencies "samtools" step_is_selected 9 && { add_unique_item core_dependencies "bwa"; add_unique_item core_dependencies "perl"; } @@ -464,6 +488,18 @@ validate_config_vars() { exit 1 fi + if step_is_selected 6; then + local declaration + for var in STAR_COMMON_ARGS STAR_TOTAL_RNA_ARGS STAR_SMRNA_SE_ARGS STAR_SMRNA_PE_ARGS \ + STAR_EXTRA_ARGS STAR_TOTAL_RNA_EXTRA_ARGS STAR_SMRNA_EXTRA_ARGS; do + declaration="$(declare -p "$var" 2>/dev/null || true)" + if [[ "$declaration" != "declare -a "* ]]; then + log 4 "ERROR" "STAR configuration variable must be a Bash indexed array: '${var}'" + missing_vars=$((missing_vars + 1)) + fi + done + fi + for var in "${path_vars[@]}"; do path="${!var:-}" if [[ -z "$path" ]]; then @@ -522,6 +558,37 @@ validate_config_vars() { log 2 "INFO" "Selected-step configuration variables validated." return 0 } + +set_star_config_defaults() { + # Backward compatibility for custom configuration files created before v1.1.0. + declare -p STAR_COMMON_ARGS &>/dev/null || STAR_COMMON_ARGS=( + --readFilesCommand zcat --outFilterMultimapNmax 100 + --winAnchorMultimapNmax 100 --outSAMtype BAM SortedByCoordinate + ) + declare -p STAR_TOTAL_RNA_ARGS &>/dev/null || STAR_TOTAL_RNA_ARGS=( + --twopassMode Basic --chimSegmentMin 10 --chimJunctionOverhangMin 10 + --chimScoreMin 1 --chimOutType Junctions WithinBAM + ) + declare -p STAR_SMRNA_SE_ARGS &>/dev/null || STAR_SMRNA_SE_ARGS=( + --twopassMode None --alignEndsType EndToEnd --alignIntronMax 1 + --outFilterMismatchNmax 1 --outFilterMultimapScoreRange 0 + --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 + --outFilterScoreMin 12 --outFilterMatchNmin 14 + --chimSegmentMin 10 --chimJunctionOverhangMin 10 + --chimScoreMin 1 --chimOutType Junctions WithinBAM + ) + declare -p STAR_SMRNA_PE_ARGS &>/dev/null || STAR_SMRNA_PE_ARGS=( + --twopassMode None --alignEndsType EndToEnd --alignIntronMax 1 + --outFilterMismatchNmax 1 --outFilterMultimapScoreRange 0 + --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 + --outFilterScoreMin 24 --outFilterMatchNmin 28 + --chimSegmentMin 10 --chimJunctionOverhangMin 10 + --chimScoreMin 1 --chimOutType Junctions WithinBAM + ) + declare -p STAR_EXTRA_ARGS &>/dev/null || STAR_EXTRA_ARGS=() + declare -p STAR_TOTAL_RNA_EXTRA_ARGS &>/dev/null || STAR_TOTAL_RNA_EXTRA_ARGS=() + declare -p STAR_SMRNA_EXTRA_ARGS &>/dev/null || STAR_SMRNA_EXTRA_ARGS=() +} assert_file_exists() { local filepath="${1:-}" local description="${2:-file}" @@ -569,6 +636,15 @@ assert_valid_bam() { fi return 0 } +link_or_copy_file() { + local source_path="${1:-}" + local destination_path="${2:-}" + assert_file_exists "$source_path" "source file for no-dedup pass-through" + rm -f -- "$destination_path" + if ! ln -- "$source_path" "$destination_path" 2>/dev/null; then + cp -- "$source_path" "$destination_path" + fi +} assert_dir_exists() { local dirpath="${1:-}" local description="${2:-directory}" @@ -635,6 +711,7 @@ current_run_signature() { printf 'version=%s\n' "$VERSION" printf 'step=%s\n' "$step_name" printf 'sample=%s\n' "$sample_name" + printf 'kit=%s\n' "$kit" printf 'is_paired_end=%s\n' "$is_paired_end" printf 'strand=%s\n' "$strand" printf 'gDNA_correction=%s\n' "$gDNA_correction" @@ -669,8 +746,8 @@ write_step_meta() { printf ' "config_file": "%s",\n' "$(json_escape "$config_file")" printf ' "fastq_read1": "%s",\n' "$(json_escape "$fastq_read1")" printf ' "fastq_read2": "%s",\n' "$(json_escape "$fastq_read2")" - printf ' "parameters": {"strand": "%s", "gDNA_correction": "%s", "circ_tool": "%s", "read_count_mode": "%s", "threads": "%s"}\n' \ - "$(json_escape "$strand")" "$(json_escape "$gDNA_correction")" "$(json_escape "$circ_tool")" "$(json_escape "$read_count_mode")" "$(json_escape "$thread_count")" + printf ' "parameters": {"kit": "%s", "strand": "%s", "gDNA_correction": "%s", "circ_tool": "%s", "read_count_mode": "%s", "threads": "%s"}\n' \ + "$(json_escape "$kit")" "$(json_escape "$strand")" "$(json_escape "$gDNA_correction")" "$(json_escape "$circ_tool")" "$(json_escape "$read_count_mode")" "$(json_escape "$thread_count")" printf '}\n' } > "$meta_file" } @@ -802,8 +879,31 @@ declare -A STEP_DESCRIPTIONS=( [24]="MultiQC comprehensive QC summary" [25]="Coverage analysis and BigWig generation (EMapper)" [26]="Coverage density plots (RNA types and meta-gene regions)" - [27]="Final HTML report generation" + [27]="Detailed and concise standalone HTML report generation" ) +apply_kit_profile() { + local strand_was_set="${1:-false}" + case "$kit" in + smarter_total_rna) + alignment_tag="umi_dedup" + result_dedup_tag="dedup" + ;; + smarter_smrna) + alignment_tag="no_dedup" + result_dedup_tag="no_dedup" + if [[ "$strand_was_set" != "true" ]]; then + strand="forward" + fi + STEP_DESCRIPTIONS[2]="SMARTer smRNA library-structure profile (no UMI)" + STEP_DESCRIPTIONS[3]="SMARTer smRNA 5-prime/poly(A) and quality trimming" + STEP_DESCRIPTIONS[6]="STAR alignment without UMI deduplication" + ;; + *) + log 5 "FATAL" "Internal error: unsupported normalized kit '${kit}'" + exit 1 + ;; + esac +} print_pipeline_steps() { log 2 "INFO" "EVscope Pipeline Steps (Version: ${VERSION})" log 2 "INFO" "========================================" @@ -933,10 +1033,10 @@ get_circ_expr_matrix() { circ_matrix="${output_dir}/Step_10_circRNA_Merge/${sample_name}_combined_CIRCexplorer2_CIRI2.tsv" ;; CIRCexplorer2) - circ_matrix="${output_dir}/Step_08_CIRCexplorer2_circRNA/${sample_name}_CIRCexplorer2_dedup_junction_readcounts_CPM.tsv" + circ_matrix="${output_dir}/Step_08_CIRCexplorer2_circRNA/${sample_name}_CIRCexplorer2_${result_dedup_tag}_junction_readcounts_CPM.tsv" ;; CIRI2) - circ_matrix="${output_dir}/Step_09_CIRI2_circRNA/${sample_name}_CIRI2_dedup_junction_readcounts_CPM.tsv" + circ_matrix="${output_dir}/Step_09_CIRI2_circRNA/${sample_name}_CIRI2_${result_dedup_tag}_junction_readcounts_CPM.tsv" ;; esac echo "$circ_matrix" @@ -967,6 +1067,29 @@ run_step_1() { } _step_2_impl() { local step_dir="${output_dir}/Step_02_UMI_Analysis" + if [[ "$kit" == "smarter_smrna" ]]; then + local input_layout="Single-end Read 1" + [[ "$is_paired_end" == "true" ]] && input_layout="Paired-end Read 1 / Read 2" + { + printf 'Parameter\tValue\n' + printf 'Kit\tSMARTer smRNA-Seq Kit for Illumina\n' + printf 'Input layout\t%s\n' "$input_layout" + printf 'Read orientation\tSense\n' + printf 'UMI\tNo\n' + printf 'R1 structure\t3-nt template-switch sequence + RNA insert + poly(A)\n' + printf '5-prime fixed trim\t3 nt\n' + printf '3-prime adapter\tAAAAAAAAAA\n' + if [[ "$is_paired_end" == "true" ]]; then + printf 'R2 structure\t14-nt poly(T) + reverse-complement RNA insert + 3-nt template-switch sequence + Illumina adapter\n' + printf 'R2 5-prime fixed trim\t14 nt\n' + printf 'R2 3-prime adapter\tNNNAGATCGGAAGAGC\n' + fi + printf 'Minimum retained length\t15 nt\n' + printf 'Molecular deduplication\tDisabled\n' + } > "${step_dir}/${sample_name}_library_profile.tsv" + return 0 + fi + local input_fq="${fastq_read2:-${fastq_read1}}" run_python "${EVscope_PATH}/bin/Step_02_plot_fastq2UMI_motif.py" \ -head "${sample_name}" -fq "${input_fq}" -n 14 -r 1000000 -o "$step_dir" @@ -988,6 +1111,44 @@ _step_3_impl() { local r2_umi_trim="${step_dir}/${sample_name}_R2_adapter_UMI_trimmed.fq.gz" local r1_clean="${step_dir}/${sample_name}_R1_clean.fq.gz" local r2_clean="${step_dir}/${sample_name}_R2_clean.fq.gz" + + if [[ "$kit" == "smarter_smrna" ]]; then + # Takara's recommended structure trimming is: -m 15 -u 3 -a AAAAAAAAAA. + # Quality trimming remains a separate pass so the length plot can distinguish + # protocol-derived trimming from the pipeline's general Q20 filtering. + if [[ "$is_paired_end" == "true" ]]; then + cutadapt -u 3 -U 14 \ + -a 'AAAAAAAAAA;min_overlap=3' \ + -A 'NNNAGATCGGAAGAGC;min_overlap=10' \ + --minimum-length 15 -j "$thread_count" \ + --json "${step_dir}/${sample_name}_structure_trimming.cutadapt.json" \ + -o "$r1_trim" -p "$r2_trim" "$fastq_read1" "$fastq_read2" + cutadapt -q 20 -Q 20 --minimum-length 15 -j "$thread_count" \ + --json "${step_dir}/${sample_name}_quality_trimming.cutadapt.json" \ + -o "$r1_clean" -p "$r2_clean" "$r1_trim" "$r2_trim" + run_python "${EVscope_PATH}/bin/Step_03_plot_fastq_read_length_dist.py" \ + --input_fastqs "$r1_trim" "$r1_clean" "$r2_clean" \ + --titles "R1 Structure-Trimmed" "R1 Clean" "R2 Clean" \ + --output_pdf "${step_dir}/${sample_name}_read_length_distribution.pdf" \ + --output_png "${step_dir}/${sample_name}_read_length_distribution.png" + assert_valid_gzip "$r2_clean" "SMARTer smRNA clean R2 FASTQ" + else + cutadapt -u 3 -a AAAAAAAAAA --overlap 3 --minimum-length 15 \ + -j "$thread_count" --json "${step_dir}/${sample_name}_structure_trimming.cutadapt.json" \ + -o "$r1_trim" "$fastq_read1" + cutadapt -q 20 --minimum-length 15 -j "$thread_count" \ + --json "${step_dir}/${sample_name}_quality_trimming.cutadapt.json" \ + -o "$r1_clean" "$r1_trim" + run_python "${EVscope_PATH}/bin/Step_03_plot_fastq_read_length_dist.py" \ + --input_fastqs "$fastq_read1" "$r1_trim" "$r1_clean" \ + --titles "Raw R1" "Structure-Trimmed R1" "Clean R1" \ + --output_pdf "${step_dir}/${sample_name}_read_length_distribution.pdf" \ + --output_png "${step_dir}/${sample_name}_read_length_distribution.png" + fi + assert_valid_gzip "$r1_clean" "SMARTer smRNA clean R1 FASTQ" + return 0 + fi + if [[ "$is_paired_end" == "true" ]]; then umi_tools extract --bc-pattern='NNNNNNNNNNNNNN' \ --stdin="${fastq_read2}" --stdout="$r2_umi" \ @@ -1066,85 +1227,108 @@ run_step_5() { local step_dir="${output_dir}/Step_05_Bacterial_Filter" run_step "$step_dir" "false" _step_5_impl } +configure_star_protocol_args() { + if [[ "$kit" == "smarter_smrna" ]]; then + if [[ "$is_paired_end" == "true" ]]; then + STAR_PROTOCOL_ARGS=("${STAR_COMMON_ARGS[@]}" "${STAR_SMRNA_PE_ARGS[@]}" \ + "${STAR_EXTRA_ARGS[@]}" "${STAR_SMRNA_EXTRA_ARGS[@]}") + else + STAR_PROTOCOL_ARGS=("${STAR_COMMON_ARGS[@]}" "${STAR_SMRNA_SE_ARGS[@]}" \ + "${STAR_EXTRA_ARGS[@]}" "${STAR_SMRNA_EXTRA_ARGS[@]}") + fi + else + STAR_PROTOCOL_ARGS=("${STAR_COMMON_ARGS[@]}" "${STAR_TOTAL_RNA_ARGS[@]}" \ + "${STAR_EXTRA_ARGS[@]}" "${STAR_TOTAL_RNA_EXTRA_ARGS[@]}") + fi +} _step_6_initial_impl() { local step_dir="${output_dir}/Step_06_Alignment_Initial" local r1_clean="${output_dir}/Step_03_UMI_Adaptor_Trim/${sample_name}_R1_clean.fq.gz" local r2_clean="${output_dir}/Step_03_UMI_Adaptor_Trim/${sample_name}_R2_clean.fq.gz" local initial_bam="${step_dir}/${sample_name}_Aligned.sortedByCoord.out.bam" - local dedup_bam="${step_dir}/${sample_name}_Aligned.sortedByCoord_umi_dedup.out.bam" - local r1_dedup="${step_dir}/${sample_name}_R1_umi_dedup.clean.fq.gz" - local r2_dedup="${step_dir}/${sample_name}_R2_umi_dedup.clean.fq.gz" + local processed_bam="${step_dir}/${sample_name}_Aligned.sortedByCoord_${alignment_tag}.out.bam" + local r1_processed="${step_dir}/${sample_name}_R1_${alignment_tag}.clean.fq.gz" + local r2_processed="${step_dir}/${sample_name}_R2_${alignment_tag}.clean.fq.gz" assert_file_exists "$r1_clean" "R1 clean FASTQ" mkdir -p "$step_dir" + configure_star_protocol_args if [[ "$is_paired_end" == "true" ]]; then assert_file_exists "$r2_clean" "R2 clean FASTQ" STAR --genomeDir "$STAR_INDEX" --readFilesIn "$r1_clean" "$r2_clean" \ --outFileNamePrefix "${step_dir}/${sample_name}_" \ - --runThreadN "$thread_count" --twopassMode Basic --runMode alignReads \ - --readFilesCommand zcat --outFilterMultimapNmax 100 --winAnchorMultimapNmax 100 \ - --outSAMtype BAM SortedByCoordinate --chimSegmentMin 10 --chimJunctionOverhangMin 10 \ - --chimScoreMin 1 --chimOutType Junctions WithinBAM --outBAMsortingThreadN "$thread_count" + --runThreadN "$thread_count" --runMode alignReads \ + --outBAMsortingThreadN "$thread_count" \ + "${STAR_PROTOCOL_ARGS[@]}" samtools index -@ "$thread_count" "$initial_bam" assert_valid_bam "$initial_bam" "initial STAR BAM" assert_nonempty_file "${initial_bam}.bai" "initial STAR BAM index" - umi_tools dedup -I "$initial_bam" -S "$dedup_bam" \ - --log="${step_dir}/${sample_name}_umi_dedup.log" --extract-umi-method=read_id --paired - samtools index -@ "$thread_count" "$dedup_bam" - samtools fastq -@ "$thread_count" -1 "$r1_dedup" -2 "$r2_dedup" \ - -0 /dev/null -s /dev/null -n "$dedup_bam" + if [[ "$kit" == "smarter_total_rna" ]]; then + umi_tools dedup -I "$initial_bam" -S "$processed_bam" \ + --log="${step_dir}/${sample_name}_umi_dedup.log" --extract-umi-method=read_id --paired + samtools index -@ "$thread_count" "$processed_bam" + samtools fastq -@ "$thread_count" -1 "$r1_processed" -2 "$r2_processed" \ + -0 /dev/null -s /dev/null -n "$processed_bam" + else + link_or_copy_file "$initial_bam" "$processed_bam" + link_or_copy_file "${initial_bam}.bai" "${processed_bam}.bai" + link_or_copy_file "$r1_clean" "$r1_processed" + link_or_copy_file "$r2_clean" "$r2_processed" + fi else STAR --genomeDir "$STAR_INDEX" --readFilesIn "$r1_clean" \ --outFileNamePrefix "${step_dir}/${sample_name}_" \ - --runThreadN "$thread_count" --twopassMode Basic --runMode alignReads \ - --readFilesCommand zcat --outFilterMultimapNmax 100 --winAnchorMultimapNmax 100 \ - --outSAMtype BAM SortedByCoordinate --chimSegmentMin 10 --chimJunctionOverhangMin 10 \ - --chimScoreMin 1 --chimOutType Junctions WithinBAM --outBAMsortingThreadN "$thread_count" + --runThreadN "$thread_count" --runMode alignReads \ + --outBAMsortingThreadN "$thread_count" \ + "${STAR_PROTOCOL_ARGS[@]}" samtools index -@ "$thread_count" "$initial_bam" assert_valid_bam "$initial_bam" "initial STAR BAM" assert_nonempty_file "${initial_bam}.bai" "initial STAR BAM index" - umi_tools dedup -I "$initial_bam" -S "$dedup_bam" \ - --log="${step_dir}/${sample_name}_umi_dedup.log" --extract-umi-method=read_id - samtools index -@ "$thread_count" "$dedup_bam" - samtools fastq -@ "$thread_count" "$dedup_bam" | gzip > "$r1_dedup" - fi - assert_valid_bam "$dedup_bam" "UMI deduplicated BAM" - assert_nonempty_file "${dedup_bam}.bai" "UMI deduplicated BAM index" - assert_valid_gzip "$r1_dedup" "R1 UMI deduplicated FASTQ" + if [[ "$kit" == "smarter_total_rna" ]]; then + umi_tools dedup -I "$initial_bam" -S "$processed_bam" \ + --log="${step_dir}/${sample_name}_umi_dedup.log" --extract-umi-method=read_id + samtools index -@ "$thread_count" "$processed_bam" + samtools fastq -@ "$thread_count" "$processed_bam" | gzip > "$r1_processed" + else + link_or_copy_file "$initial_bam" "$processed_bam" + link_or_copy_file "${initial_bam}.bai" "${processed_bam}.bai" + link_or_copy_file "$r1_clean" "$r1_processed" + fi + fi + assert_valid_bam "$processed_bam" "kit-processed BAM" + assert_nonempty_file "${processed_bam}.bai" "kit-processed BAM index" + assert_valid_gzip "$r1_processed" "R1 kit-processed FASTQ" if [[ "$is_paired_end" == "true" ]]; then - assert_valid_gzip "$r2_dedup" "R2 UMI deduplicated FASTQ" + assert_valid_gzip "$r2_processed" "R2 kit-processed FASTQ" fi } _step_6_refined_impl() { local step_dir="${output_dir}/Step_06_Alignment_Refined" - local r1_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R1_umi_dedup.clean.fq.gz" - local r2_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R2_umi_dedup.clean.fq.gz" - local final_bam="${step_dir}/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" - assert_file_exists "$r1_dedup" "R1 dedup FASTQ from Step 6 initial" + local r1_processed="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R1_${alignment_tag}.clean.fq.gz" + local r2_processed="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R2_${alignment_tag}.clean.fq.gz" + local final_bam="${step_dir}/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" + assert_file_exists "$r1_processed" "R1 processed FASTQ from Step 6 initial" mkdir -p "$step_dir" + configure_star_protocol_args if [[ "$is_paired_end" == "true" ]]; then - assert_file_exists "$r2_dedup" "R2 dedup FASTQ from Step 6 initial" - STAR --genomeDir "$STAR_INDEX" --readFilesIn "$r1_dedup" "$r2_dedup" \ - --outFileNamePrefix "${step_dir}/${sample_name}_STAR_umi_dedup_" \ - --runThreadN "$thread_count" --twopassMode Basic --runMode alignReads \ - --quantMode GeneCounts --readFilesCommand zcat --outFilterMultimapNmax 100 \ - --winAnchorMultimapNmax 100 --outSAMtype BAM SortedByCoordinate \ - --chimSegmentMin 10 --chimJunctionOverhangMin 10 --chimScoreMin 1 \ - --chimOutType Junctions WithinBAM --outBAMsortingThreadN "$thread_count" + assert_file_exists "$r2_processed" "R2 processed FASTQ from Step 6 initial" + STAR --genomeDir "$STAR_INDEX" --readFilesIn "$r1_processed" "$r2_processed" \ + --outFileNamePrefix "${step_dir}/${sample_name}_STAR_${alignment_tag}_" \ + --runThreadN "$thread_count" --runMode alignReads \ + --quantMode GeneCounts --outBAMsortingThreadN "$thread_count" \ + "${STAR_PROTOCOL_ARGS[@]}" else - STAR --genomeDir "$STAR_INDEX" --readFilesIn "$r1_dedup" \ - --outFileNamePrefix "${step_dir}/${sample_name}_STAR_umi_dedup_" \ - --runThreadN "$thread_count" --twopassMode Basic --runMode alignReads \ - --quantMode GeneCounts --readFilesCommand zcat --outFilterMultimapNmax 100 \ - --winAnchorMultimapNmax 100 --outSAMtype BAM SortedByCoordinate \ - --chimSegmentMin 10 --chimJunctionOverhangMin 10 --chimScoreMin 1 \ - --chimOutType Junctions WithinBAM --outBAMsortingThreadN "$thread_count" + STAR --genomeDir "$STAR_INDEX" --readFilesIn "$r1_processed" \ + --outFileNamePrefix "${step_dir}/${sample_name}_STAR_${alignment_tag}_" \ + --runThreadN "$thread_count" --runMode alignReads \ + --quantMode GeneCounts --outBAMsortingThreadN "$thread_count" \ + "${STAR_PROTOCOL_ARGS[@]}" fi samtools index -@ "$thread_count" "$final_bam" samtools flagstat -@ "$thread_count" "$final_bam" > \ - "${step_dir}/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.flagstat" + "${step_dir}/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.flagstat" assert_valid_bam "$final_bam" "refined STAR BAM" assert_nonempty_file "${final_bam}.bai" "refined STAR BAM index" - assert_nonempty_file "${step_dir}/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.flagstat" "refined STAR flagstat" + assert_nonempty_file "${step_dir}/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.flagstat" "refined STAR flagstat" } _step_6_impl() { mkdir -p "${output_dir}/Step_06_Alignment_Initial" @@ -1158,9 +1342,9 @@ run_step_6() { } _step_7_impl() { local step_dir="${output_dir}/Step_07_Strand_Detection" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" assert_file_exists "$final_bam" "Aligned BAM from Step 6" - local star_log_refined="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Log.final.out" + local star_log_refined="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Log.final.out" local star_log_args=() [[ -f "$star_log_refined" ]] && star_log_args=(--star_log "$star_log_refined") run_python "${EVscope_PATH}/bin/Step_07_bam2strand.py" \ @@ -1173,11 +1357,11 @@ run_step_7() { } _step_8_impl() { local step_dir="${output_dir}/Step_08_CIRCexplorer2_circRNA" - local chimeric="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Chimeric.out.junction" + local chimeric="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Chimeric.out.junction" local bsj_bed="${step_dir}/${sample_name}_back_spliced_junction.bed" local known_circs="${step_dir}/${sample_name}_circularRNA_known.txt" - local final_out="${step_dir}/${sample_name}_CIRCexplorer2_dedup_junction_readcounts_CPM.tsv" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_out="${step_dir}/${sample_name}_CIRCexplorer2_${result_dedup_tag}_junction_readcounts_CPM.tsv" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" assert_file_exists "$chimeric" "Chimeric junctions from Step 6" assert_file_exists "$final_bam" "Aligned BAM from Step 6" CIRCexplorer2 parse -t STAR -b "$bsj_bed" "$chimeric" @@ -1196,15 +1380,15 @@ run_step_8() { } _step_9_impl() { local step_dir="${output_dir}/Step_09_CIRI2_circRNA" - local r1_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R1_umi_dedup.clean.fq.gz" - local r2_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R2_umi_dedup.clean.fq.gz" - local bwa_sam="${step_dir}/${sample_name}_umi_dedup.bwa.sam" + local r1_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R1_${alignment_tag}.clean.fq.gz" + local r2_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R2_${alignment_tag}.clean.fq.gz" + local bwa_sam="${step_dir}/${sample_name}_${alignment_tag}.bwa.sam" local ciri2_out="${step_dir}/${sample_name}_CIRI2_out.tsv" - local final_out="${step_dir}/${sample_name}_CIRI2_dedup_junction_readcounts_CPM.tsv" - assert_valid_gzip "$r1_dedup" "R1 dedup FASTQ from Step 6" + local final_out="${step_dir}/${sample_name}_CIRI2_${result_dedup_tag}_junction_readcounts_CPM.tsv" + assert_valid_gzip "$r1_dedup" "R1 processed FASTQ from Step 6" trap 'rm -f "$bwa_sam"' RETURN if [[ "$is_paired_end" == "true" ]]; then - assert_valid_gzip "$r2_dedup" "R2 dedup FASTQ from Step 6" + assert_valid_gzip "$r2_dedup" "R2 processed FASTQ from Step 6" bwa mem -t "$thread_count" -T 19 "${BWA_INDEX:-}" "$r1_dedup" "$r2_dedup" > "$bwa_sam" else bwa mem -t "$thread_count" -T 19 "${BWA_INDEX:-}" "$r1_dedup" > "$bwa_sam" @@ -1228,8 +1412,8 @@ run_step_9() { } _step_10_impl() { local step_dir="${output_dir}/Step_10_circRNA_Merge" - local cirexp="${output_dir}/Step_08_CIRCexplorer2_circRNA/${sample_name}_CIRCexplorer2_dedup_junction_readcounts_CPM.tsv" - local ciri2="${output_dir}/Step_09_CIRI2_circRNA/${sample_name}_CIRI2_dedup_junction_readcounts_CPM.tsv" + local cirexp="${output_dir}/Step_08_CIRCexplorer2_circRNA/${sample_name}_CIRCexplorer2_${result_dedup_tag}_junction_readcounts_CPM.tsv" + local ciri2="${output_dir}/Step_09_CIRI2_circRNA/${sample_name}_CIRI2_${result_dedup_tag}_junction_readcounts_CPM.tsv" assert_file_exists "$cirexp" "CIRCexplorer2 results from Step 8" assert_file_exists "$ciri2" "CIRI2 results from Step 9" run_python "${EVscope_PATH}/bin/Step_10_circRNA_merge.py" \ @@ -1247,7 +1431,7 @@ run_step_10() { } _step_11_impl() { local step_dir="${output_dir}/Step_11_RNA_Metrics" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" local picard_out="${step_dir}/${sample_name}_picard_metrics.tsv" local insert_out="${step_dir}/${sample_name}_insert_size_metrics.tsv" local insert_pdf="${step_dir}/${sample_name}_insert_size_histogram.pdf" @@ -1281,7 +1465,7 @@ run_step_11() { } _step_12_impl() { local step_dir="${output_dir}/Step_12_featureCounts_Quant" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" assert_file_exists "$final_bam" "Aligned BAM from Step 6" if [[ "$is_paired_end" == "true" ]]; then featureCounts -a "$TOTAL_GENE_GTF" -o "${step_dir}/${sample_name}_featureCounts.tsv" \ @@ -1298,7 +1482,7 @@ run_step_12() { } _step_13_impl() { local step_dir="${output_dir}/Step_13_gDNA_Corrected_Quant" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" local forward_counts="${step_dir}/${sample_name}_featureCounts_s1_forward.tsv" local reverse_counts="${step_dir}/${sample_name}_featureCounts_s2_reverse.tsv" local corrected="${step_dir}/${sample_name}_gDNA_corrected_counts.tsv" @@ -1328,14 +1512,14 @@ run_step_13() { } _step_14_impl() { local step_dir="${output_dir}/Step_14_RSEM_Quant" - local r1_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R1_umi_dedup.clean.fq.gz" - local r2_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R2_umi_dedup.clean.fq.gz" + local r1_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R1_${alignment_tag}.clean.fq.gz" + local r2_dedup="${output_dir}/Step_06_Alignment_Initial/${sample_name}_R2_${alignment_tag}.clean.fq.gz" local rsem_strand="$strand" [[ "$rsem_strand" == "unstrand" ]] && rsem_strand="none" - assert_file_exists "$r1_dedup" "R1 dedup FASTQ from Step 6" + assert_file_exists "$r1_dedup" "R1 processed FASTQ from Step 6" mkdir -p "${step_dir}/tmp" if [[ "$is_paired_end" == "true" ]]; then - assert_file_exists "$r2_dedup" "R2 dedup FASTQ from Step 6" + assert_file_exists "$r2_dedup" "R2 processed FASTQ from Step 6" perl "$RSEM_CALC_EXPR" --paired-end --bowtie2 --strandedness "$rsem_strand" \ --bowtie2-k 2 -p "$thread_count" --no-bam-output --seed 12345 \ "$r1_dedup" "$r2_dedup" "${RSEM_BOWTIE2_INDEX:-}" \ @@ -1368,7 +1552,8 @@ _run_expression_analysis_impl() { --gene_expr "$base_matrix" --circRNA_expr "$circ_matrix" --out_matrix "$combined" else log 3 "WARN" "circRNA matrix not found, using gene expression only" - cp "$base_matrix" "$combined" + run_python "${EVscope_PATH}/bin/Step_15_combine_total_RNA_expr_matrix.py" \ + --gene_expr "$base_matrix" --out_matrix "$combined" fi run_python "${EVscope_PATH}/bin/Step_15_plot_RNA_distribution_1subplot.py" \ --sample_name "${sample_name}" --Expr_matrix "$combined" \ @@ -1423,7 +1608,8 @@ _step_17_impl() { --gene_expr "$base_matrix" --circRNA_expr "$circ_matrix" --out_matrix "$combined" else log 3 "WARN" "circRNA matrix not found, using gene expression only" - cp "$base_matrix" "$combined" + run_python "${EVscope_PATH}/bin/Step_15_combine_total_RNA_expr_matrix.py" \ + --gene_expr "$base_matrix" --out_matrix "$combined" fi run_python "${EVscope_PATH}/bin/Step_15_plot_RNA_distribution_1subplot.py" \ --sample_name "${sample_name}" --Expr_matrix "$combined" \ @@ -1446,7 +1632,7 @@ run_step_17() { # Step 18: Parallel featureCounts for genomic regions - OPTIMIZED _step_18_impl() { local step_dir="${output_dir}/Step_18_Genomic_Regions" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" assert_file_exists "$final_bam" "Aligned BAM from Step 6" # Use pre-built merged SAF (7 regions, 801K features, mutually exclusive) @@ -1615,7 +1801,7 @@ run_step_23() { _step_24_impl() { local step_dir="${output_dir}/Step_24_MultiQC_Summary" local -a multiqc_dirs=() - local -a qc_args=("--output" "${step_dir}/${sample_name}_QC_summary.tsv") + local -a qc_args=("--output" "${step_dir}/${sample_name}_QC_summary.tsv" "--kit" "$kit") local -a potential_dirs=( "${output_dir}/Step_01_Raw_QC" "${output_dir}/Step_04_Trimmed_QC" @@ -1679,16 +1865,16 @@ _step_24_impl() { done < <(find "$ribo_dir" -maxdepth 1 -type f -name "${sample_name}_rRNA_*.fq.gz" -print0 2>/dev/null | sort -z || true) [[ -f "${ribo_dir}/step.done" ]] && qc_args+=("--ribo_fastqs" "${ribo_fqs[@]}") - local star_log_refined="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Log.final.out" + local star_log_refined="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Log.final.out" [[ -f "$star_log_refined" ]] && qc_args+=("--STAR_log" "$star_log_refined") local star_log_initial="${output_dir}/Step_06_Alignment_Initial/${sample_name}_Log.final.out" [[ -f "$star_log_initial" ]] && qc_args+=("--STAR_log_initial" "$star_log_initial") - local strand_file="${output_dir}/Step_07_Strand_Detection/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out_bam2strandness.tsv" + local strand_file="${output_dir}/Step_07_Strand_Detection/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out_bam2strandness.tsv" [[ -f "$strand_file" ]] && qc_args+=("--bam2strand_file" "$strand_file") - local length_strandness_file="${output_dir}/Step_07_Strand_Detection/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out_read_length_stratified_strandness.tsv" + local length_strandness_file="${output_dir}/Step_07_Strand_Detection/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out_read_length_stratified_strandness.tsv" [[ -f "$length_strandness_file" ]] && qc_args+=("--length_strandness_file" "$length_strandness_file") local picard_insert="${output_dir}/Step_11_RNA_Metrics/${sample_name}_insert_size_metrics.tsv" @@ -1750,15 +1936,19 @@ run_step_24() { } _step_25_emapper_impl() { local step_dir="${output_dir}/Step_25_EMapper_BigWig_Quantification/EMapper_output" - local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_umi_dedup_Aligned.sortedByCoord.out.bam" + local final_bam="${output_dir}/Step_06_Alignment_Refined/${sample_name}_STAR_${alignment_tag}_Aligned.sortedByCoord.out.bam" assert_file_exists "$final_bam" "Aligned BAM from Step 6" mkdir -p "$step_dir" + local -a polya_args=() + if [[ "$kit" == "smarter_total_rna" ]]; then + polya_args=(--polyA_tail_detection) + fi run_python "${EVscope_PATH}/bin/Step_25_EMapper.py" \ --num_threads "$thread_count" --input_bam "$final_bam" \ --prefix "${sample_name}" --output_dir "$step_dir" \ --strandness "$strand" --reference_fasta "$HUMAN_GENOME_FASTA" \ --gtf "$TOTAL_GENE_GTF" --positional_assignment EM \ - --gene_disambiguation EM --no-cleanup --polyA_tail_detection + --gene_disambiguation EM --no-cleanup "${polya_args[@]}" } _step_25_expression_impl() { local step_dir="${output_dir}/Step_25_EMapper_BigWig_Quantification/bigwig2expression" @@ -1883,12 +2073,34 @@ _step_27_impl() { local abs_step_dir="${abs_output_dir}/Step_27_HTML_Report" local rmd_input="${abs_bin_dir}/Step_27_html_report.Rmd" local report_file="${sample_name}_final_report.html" + local summary_script="${abs_bin_dir}/Step_27_generate_summary_report.py" + local summary_file="${sample_name}_summary_report.html" local local_rmd="${abs_step_dir}/${sample_name}_report.Rmd" assert_file_exists "$rmd_input" "R Markdown template" + assert_file_exists "$summary_script" "concise summary report generator" mkdir -p "$abs_step_dir" cp "$rmd_input" "$local_rmd" + + local qc_summary="${abs_output_dir}/Step_24_MultiQC_Summary/${sample_name}_QC_summary.tsv" + local expression_matrix="${abs_output_dir}/Step_15_featureCounts_Expression/${sample_name}_combined_expression_matrix_linearRNA_TPM_circRNA_CPM_featureCounts.tsv" + if [[ ! -f "$expression_matrix" ]]; then + expression_matrix="${abs_output_dir}/Step_15_featureCounts_Expression/${sample_name}_Gene_readcounts_normalized_expression_matrix_featureCounts.tsv" + fi + local -a summary_args=( + --sample-name "$sample_name" + --kit "$kit" + --strand "$strand" + --paired-end "$is_paired_end" + --output "${abs_step_dir}/${summary_file}" + ) + [[ -f "$qc_summary" ]] && summary_args+=(--qc-summary "$qc_summary") + [[ -f "$expression_matrix" ]] && summary_args+=(--expression-matrix "$expression_matrix") + run_python "$summary_script" "${summary_args[@]}" + assert_nonempty_file "${abs_step_dir}/${summary_file}" "standalone concise summary report" + export EVscope_OUTPUT_DIR="$abs_output_dir" export EVscope_FIGURES_DIR="$abs_figures_dir" + export EVscope_KIT="$kit" export EVscope_RMD_INPUT="$local_rmd" export EVscope_REPORT_FILE="$report_file" export EVscope_REPORT_DIR="$abs_step_dir" @@ -1916,6 +2128,7 @@ main() { local config_file_input="$config_file" local -a input_fastqs=() local dry_run="false" + local strand_was_set="false" while (( $# > 0 )); do case "$1" in @@ -1927,6 +2140,14 @@ main() { [[ ! "$2" =~ ^[0-9]+$ ]] && { log 5 "FATAL" "--threads must be a positive integer"; exit 1; } (( $2 < 1 )) && { log 5 "FATAL" "--threads must be >= 1"; exit 1; } thread_count="$2"; shift 2 ;; + --kit) + [[ -z "${2:-}" ]] && { log 5 "FATAL" "--kit requires a value"; exit 1; } + local normalized_kit="" + if ! normalized_kit="$(normalize_kit "$2")"; then + log 5 "FATAL" "--kit must be 'smarter_total_rna' or 'smarter_smrna' (received: '$2')" + exit 1 + fi + kit="$normalized_kit"; shift 2 ;; --run_steps) [[ -z "${2:-}" ]] && { log 5 "FATAL" "--run_steps requires a value"; exit 1; } run_steps="$2"; shift 2 ;; @@ -1948,7 +2169,7 @@ main() { --strand) [[ -z "${2:-}" ]] && { log 5 "FATAL" "--strand requires a value"; exit 1; } [[ ! "$2" =~ ^(reverse|forward|unstrand)$ ]] && { log 5 "FATAL" "--strand must be 'reverse', 'forward', or 'unstrand'"; exit 1; } - strand="$2"; shift 2 ;; + strand="$2"; strand_was_set="true"; shift 2 ;; --config) [[ -z "${2:-}" ]] && { log 5 "FATAL" "--config requires a value"; exit 1; } config_file="$2"; config_file_input="$2"; shift 2 ;; @@ -1989,6 +2210,10 @@ main() { esac done + apply_kit_profile "$strand_was_set" + if [[ "$kit" == "smarter_smrna" && "$strand_was_set" == "true" && "$strand" != "forward" ]]; then + log 3 "WARN" "SMARTer smRNA Read 1 is sense; explicit --strand ${strand} will override the recommended forward setting." + fi if [[ "$gDNA_correction" == "yes" && "$strand" == "unstrand" ]]; then log 5 "FATAL" "Logic Error: --gDNA_correction yes requires --strand forward/reverse." exit 1 @@ -2006,7 +2231,6 @@ main() { log 5 "FATAL" "Provided: ${#input_fastqs[@]} files" exit 1 fi - output_dir="${sample_name}_EVscope_output" if [[ -d "$output_dir" ]]; then if [[ -n "$(find "$output_dir" -mindepth 1 -maxdepth 1 -print -quit 2>/dev/null)" && "$resume_mode" != "true" && "$force_mode" != "true" ]]; then @@ -2073,6 +2297,7 @@ main() { # shellcheck source=/dev/null source "$config_file" log 2 "INFO" "Loaded configuration from: ${config_file}" + set_star_config_defaults validate_config_vars check_dependencies check_conda_envs @@ -2101,6 +2326,7 @@ main() { log 2 "INFO" "EVscope Pipeline v${VERSION}" log 2 "INFO" "============================================================" log 2 "INFO" "Sample name: ${sample_name}" + log 2 "INFO" "Library kit: ${kit}" log 2 "INFO" "Output directory: ${output_dir} (${abs_output_dir})" log 2 "INFO" "Paired-end mode: ${is_paired_end}" log 2 "INFO" "Thread count: ${thread_count}" @@ -2144,4 +2370,4 @@ main() { # ============================================================================== # SCRIPT ENTRY POINT # ============================================================================== -main "$@" \ No newline at end of file +main "$@" diff --git a/README.md b/README.md index 29476ed..b09c8de 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ -# EVscope: A Comprehensive Bioinformatics Workflow for Accurate and Robust Analysis of Total RNA Sequencing from Extracellular Vesicles +# EVscope: A Comprehensive Bioinformatics Workflow for Accurate and Robust Analysis of RNA Sequencing from Extracellular Vesicles -**EVscope** is an open-source, modular bioinformatics pipeline designed for the analysis of extracellular vesicle (EV)-enriched total RNA sequencing data. Tailored to EV RNA-seq challenges - low RNA yield, fragmented inserts, diverse RNA biotypes, high multi-mapping, and contamination risk - EVscope processes paired-end or single-end FASTQ files through an end-to-end workflow. It includes quality control, UMI-based deduplication, two-pass STAR alignment, circular RNA detection, expression matrix generation, contamination screening, exploratory source-enrichment analysis, and comprehensive reporting. Optimized for the SMARTer Stranded Total RNA-Seq Kit v3 (Pico Input), EVscope introduces EMapper, an expectation-maximization (EM) module whose main EMapper-specific contribution in EVscope is EM-weighted genome-coordinate BigWig/coverage generation with RNA annotation support; gene-level count concordance with featureCounts/RSEM is used as a sanity check, not as a claim that EMapper is superior to RSEM for conventional gene readcount quantification. +**EVscope** is an open-source, modular bioinformatics pipeline designed for extracellular-vesicle (EV) RNA sequencing data. It provides kit-aware processing for the SMARTer Stranded Total RNA-Seq Kit v3 (Pico Input) and the SMARTer smRNA-Seq Kit for Illumina. Select the protocol with `--kit smarter_total_rna` (the backward-compatible default) or `--kit smarter_smrna`. Tailored to EV RNA-seq challenges - low RNA yield, fragmented inserts, diverse RNA biotypes, high multi-mapping, and contamination risk - EVscope processes supported paired-end or single-end FASTQ inputs through an end-to-end workflow. It includes quality control, kit-appropriate duplicate handling, kit-specific STAR alignment, circular RNA detection, expression matrix generation, contamination screening, exploratory source-enrichment analysis, and comprehensive reporting. EVscope also introduces EMapper, an expectation-maximization (EM) module whose main EMapper-specific contribution is EM-weighted genome-coordinate BigWig/coverage generation with RNA annotation support; gene-level count concordance with featureCounts/RSEM is used as a sanity check, not as a claim that EMapper is superior to RSEM for conventional gene readcount quantification.

EVscope Pipeline Overview @@ -29,22 +29,23 @@ ## Key Features -- **Novel Read-Through Detection**: Trims UMI-derived adapter sequences from Read1 using reverse-complemented Read2 UMIs (`bin/Step_03_UMIAdapterTrimR1.py`). +- **Kit-Aware Preprocessing**: Preserves the existing 14-bp UMI workflow for `smarter_total_rna`, while `smarter_smrna` uses the official SMARTer smRNA 5'/poly(A) trimming profile without UMI extraction or deduplication. +- **Novel Read-Through Detection for Total RNA**: Trims UMI-derived adapter sequences from Read1 using reverse-complemented Read2 UMIs (`bin/Step_03_UMIAdapterTrimR1.py`) when `--kit smarter_total_rna` is selected. - **EM-Weighted BigWig Coverage Profiling**: Uses a genome-wide expectation-maximization algorithm to assign multi-mapped fragments at single-base or binned resolution and generate strand-aware, RPM-normalized BigWig tracks (`bin/Step_25_EMapper.py`). - **Comprehensive RNA Annotation**: Supports 3,659,642 RNAs across 20 biotypes (e.g., protein-coding, lncRNAs, miRNAs, piRNAs, retrotransposons) from GENCODE v45, piRBase v3.0, and RepeatMasker. -- **Dual circRNA Detection**: Integrates CIRCexplorer2 and CIRI2 for robust circular RNA identification, with merged results for enhanced sensitivity (`bin/Step_10_circRNA_merge.py`). +- **Dual circRNA Detection**: Integrates CIRCexplorer2 and CIRI2 for circular RNA identification, with merged results (`bin/Step_10_circRNA_merge.py`); SMARTer smRNA libraries use the no-dedup reads and kit-specific STAR chimeric output. - **Tissue Deconvolution**: Infers EV RNA cellular origins using GTEx v10 and Human Brain Cell Atlas v1.0 references (`bin/Step_22_run_RNA_deconvolution_ARIC.py`). - **Contamination Screening**: Filters bacterial (BBSplit) and microbial (Kraken2) contamination, with optional genomic DNA correction via strand-specific subtraction. -- **Extensive Quality Control**: Validates raw and trimmed FASTQs (FastQC), UMI motifs, and alignment metrics (`bin/Step_24_generate_QC_matrix.py`). +- **Extensive Quality Control**: Validates raw and trimmed FASTQs (FastQC), generates a kit-aware library profile, and summarizes alignment metrics (`bin/Step_24_generate_QC_matrix.py`). - **Expression Quantification**: Produces TPM/CPM matrices using featureCounts and RSEM, with RNA distribution visualizations (`bin/Step_15_plot_RNA_distribution_*.py`). -- **Interactive Reporting**: Generates bigWig tracks, density plots, and a comprehensive HTML report via R Markdown (`bin/Step_27_html_report.Rmd`). +- **Interactive Reporting**: Generates bigWig tracks, density plots, a comprehensive HTML report via R Markdown (`bin/Step_27_html_report.Rmd`), and a separate one-paragraph RNA-seq summary covering read retention, human/E. coli/Mycoplasma mapping, alignment type, genomic regions, strand orientation, and the top 10 genes. - **Reproducibility**: Single-command Bash script with Conda environments, containerization support, and detailed logging. ## Motivation Extracellular vesicles (EVs) are critical mediators of intercellular communication, carrying diverse RNAs that serve as potential biomarkers for diseases like cancer and neurodegeneration. However, EV RNA sequencing faces unique challenges: low RNA abundance, fragmented transcripts, contamination from genomic DNA or bacterial RNA, and the presence of non-polyadenylated RNAs (e.g., miRNAs, lncRNAs). Standard RNA-seq pipelines, designed for cellular RNA, often fail to address these issues, leading to unreliable results due to multi-mapping reads, incomplete RNA annotations, or unfiltered contaminants. -EVscope provides a specialized, end-to-end pipeline optimized for EV-enriched total RNA-seq. Its distinctive contribution is not to replace transcript quantifiers such as RSEM for ordinary gene readcount estimation, but to connect EV-oriented QC, broad RNA annotation, EM-weighted genome-coordinate BigWig tracks, RNA-biotype/meta-gene coverage profiling, and report generation in one reproducible workflow. +EVscope provides a specialized, end-to-end pipeline with profiles for EV-enriched SMARTer Total RNA and SMARTer smRNA sequencing. Its distinctive contribution is not to replace transcript quantifiers such as RSEM for ordinary gene readcount estimation, but to connect kit-aware EV RNA preprocessing, broad RNA annotation, EM-weighted genome-coordinate BigWig tracks, RNA-biotype/meta-gene coverage profiling, and report generation in one reproducible workflow. ## Directory Structure @@ -53,7 +54,7 @@ The EVscope repository is organized as follows: ``` EVscope/ ├── EVscope.conf # Configuration file for tool and reference paths -├── EVscope.sh # Main pipeline script (v1.0.0) +├── EVscope.sh # Main pipeline script (v1.1.0) ├── README.md # This documentation ├── bin/ # Custom scripts for pipeline steps │ ├── Step_02_calculate_ACC_motif_fraction.py # Calculates ACC motif fractions @@ -79,7 +80,8 @@ EVscope/ │ ├── Step_25_EMapper.py # EM-based read coverage estimation │ ├── Step_26_density_plot_over_meta_gene.sh # Density plots for meta-gene regions │ ├── Step_26_density_plot_over_RNA_types.sh # Density plots for RNA types -│ └── Step_27_html_report.Rmd # Generates HTML report +│ ├── Step_27_generate_summary_report.py # Generates standalone concise summary +│ └── Step_27_html_report.Rmd # Generates comprehensive HTML report ├── figures/ # Pipeline visualization │ └── EVscope_pipeline.png # Pipeline overview image ├── references/ # Reference genomes, annotations, and indices @@ -157,7 +159,7 @@ EVscope/ ## Smoke Test -A lightweight repository smoke test is available under `tests/smoke/`. It validates script syntax and selected toy-data transformations without requiring the full human reference bundle. +A lightweight repository smoke test is available under `tests/smoke/`. It validates script syntax, kit-aware CLI behavior, summary-report logic, and selected synthetic-data transformations without requiring the full human reference bundle. The SMARTer smRNA pipeline checks use test doubles for cutadapt, STAR, and samtools; they verify trimming fixtures and command routing, rather than alignment accuracy on biological data. ```bash bash tests/smoke/run_smoke.sh @@ -168,39 +170,49 @@ For a full end-to-end run, download the Zenodo reference bundle and use the SRA ## Usage ### Command Syntax + ```bash -bash EVscope.sh --sample_name --input_fastqs [options] +bash EVscope.sh --sample_name --input_fastqs [--kit ] [options] ``` **Required Arguments**: + - `--sample_name `: Unique sample identifier (used for output files). -- `--input_fastqs `: FASTQ file paths: one file for single-end input or two files for paired-end input. Space-separated (`R1.fq.gz R2.fq.gz`) and comma-separated (`R1.fq.gz,R2.fq.gz`) forms are accepted; R1 must be listed before R2. +- `--input_fastqs `: FASTQ file paths. Both kit profiles accept one single-end file or two paired-end files. Paired-end inputs may be space-separated (`R1.fq.gz R2.fq.gz`) or comma-separated (`R1.fq.gz,R2.fq.gz`), with R1 first. **Optional Arguments**: + | Option | Description | Default | |--------|-------------|---------| | `--threads ` | Number of CPU threads | 1 | +| `--kit ` | Library kit (`smarter_total_rna` or `smarter_smrna`) | `smarter_total_rna` | | `--run_steps ` | Steps to run (e.g., `1,3,5-8`, `all`) | `all` | | `--skip_steps ` | Steps to skip (e.g., `2,4`) | None | | `--circ_tool ` | circRNA detection tool (`CIRCexplorer2`, `CIRI2`, `both`) | `both` | | `--gDNA_correction ` | Apply genomic DNA correction | `no` | -| `--strand ` | Library strandedness (`forward`, `reverse`, `unstrand`) | `unstrand` | +| `--strand ` | Library strandedness (`forward`, `reverse`, `unstrand`); when omitted for `smarter_smrna`, R1 is treated as sense and `forward` is selected | `unstrand` for Total RNA; `forward` for smRNA | | `--config ` | Custom configuration file | `EVscope.conf` | | `--resume` | Resume an existing output directory only when step metadata match current inputs/config/parameters | off | | `--force` | Allow re-running selected steps in an existing output directory | off | | `--dry-run` | Validate inputs/config/dependencies and print execution plan without running steps | off | | `-V, --verbosity ` | Logging level (1=DEBUG, 2=INFO, 3=WARN, 4=ERROR, 5=FATAL) | 2 | +| `-h, --help` | Display help message | - | +| `-v, --version` | Show pipeline version | - | -**Resume behavior**: EVscope refuses to reuse a non-empty output directory unless `--resume` or `--force` is provided. `--resume` only skips steps whose `step.meta.json` matches the current inputs/config/parameters. Legacy output directories created before this metadata file existed may require `--force` or cleanup before rerunning. +**Resume behavior**: EVscope refuses to reuse a non-empty output directory unless `--resume` or `--force` is provided. `--resume` only skips steps whose `step.meta.json` matches the current inputs, configuration, pipeline version, and parameters including `--kit`. EVscope v1.1 adds the kit to the run signature, so output directories created by v1.0 have incompatible signatures and must be regenerated with `--force` (or removed and rerun); they cannot be safely resumed with v1.1. **Report R environment**: Step 27 can render through a configured Conda R environment. Set `REPORT_R_ENV` in `EVscope.conf` to a Conda R environment name (for example, `r_env_4.4.3`) and, if needed, `CONDA_EXE` to the Conda executable. The default public config leaves `REPORT_R_ENV` empty and uses `Rscript` from `PATH`. -| `-h, --help` | Display help message | - | -| `-v, --version` | Show pipeline version | - | -### Example: Full Pipeline +**Local configuration**: Copy `EVscope.conf` to `EVscope.local.conf`, set your installation and reference paths there, and pass `--config EVscope.local.conf`. This local file and sample output directories (`*_EVscope_output/`) are ignored by Git. The public configuration uses placeholder paths; committed FASTQ fixtures contain synthetic reads only. + +**STAR parameter injection**: Step 6 reads STAR options from Bash arrays in `EVscope.conf`. `STAR_COMMON_ARGS` applies to all libraries; `STAR_TOTAL_RNA_ARGS`, `STAR_SMRNA_SE_ARGS`, and `STAR_SMRNA_PE_ARGS` define the protocol profiles. Add new options to `STAR_EXTRA_ARGS`, `STAR_TOTAL_RNA_EXTRA_ARGS`, or `STAR_SMRNA_EXTRA_ARGS`. Keep each option and value as separate array elements, for example `STAR_SMRNA_EXTRA_ARGS=(--limitBAMsortRAM 30000000000)`. To change an option already present, edit it in its profile array instead of adding the same option again. Input/output paths, threads, run mode, sorting threads, and refined-run gene counting remain pipeline-controlled. + +### Example: SMARTer Stranded Total RNA-Seq Kit v3 + ```bash -bash EVscope.sh --sample_name Example_Data \ +bash EVscope.sh --sample_name Total_RNA_Example \ --input_fastqs R1.fq.gz R2.fq.gz \ + --kit smarter_total_rna \ --threads 20 \ --run_steps all \ --gDNA_correction yes \ @@ -208,22 +220,45 @@ bash EVscope.sh --sample_name Example_Data \ --verbosity 2 ``` +Omitting `--kit` gives the same `smarter_total_rna` behavior for backward compatibility. + +### Example: SMARTer smRNA-Seq Kit for Illumina + +```bash +bash EVscope.sh --sample_name smRNA_Example \ + --input_fastqs R1.fq.gz R2.fq.gz \ + --kit smarter_smrna \ + --threads 20 \ + --run_steps all \ + --verbosity 2 +``` + +In this example `--strand` is intentionally omitted, so EVscope uses `forward` because SMARTer smRNA R1 is the sense read. +Steps 8-10 remain available for this kit. EVscope runs circRNA detection on the no-dedup processed reads and retains STAR chimeric-junction output for CIRCexplorer2. Because the library's 3'-polyadenylation workflow is not optimized for intact circRNA capture, interpret circRNA calls cautiously and validate important back-splice junctions independently. + ## Input Data Format -- **FASTQ Files**: Gzipped, paired-end (`R1.fastq.gz`, `R2.fastq.gz`) or single-end. -- **Sequencing Protocol**: Optimized for SMARTer Stranded Total RNA-Seq Kit v3 (Pico Input) with 14-bp UMIs in Read2. +- **FASTQ files**: Gzipped FASTQ input with R1 listed first whenever paired-end input is permitted. - **Quality**: High-quality reads suitable for EV RNA-seq. +| Kit | `--kit` value | Accepted input | Orientation and default strand | Steps 2-3 and duplicate handling | +|-----|---------------|----------------|--------------------------------|--------------------------------| +| SMARTer Stranded Total RNA-Seq Kit v3 (Pico Input) | `smarter_total_rna` | Single-end or paired-end; paired-end requires R1 then R2 | Existing strand behavior; default `unstrand` | Existing UMI/ACC profiling and 14-bp UMI workflow, adapter and Q20 trimming, read-through UMI removal, and UMI-based deduplication | +| SMARTer smRNA-Seq Kit for Illumina | `smarter_smrna` | Single-end R1 or paired-end R1/R2 | R1 is sense; if `--strand` is omitted, EVscope selects `forward` | R1 uses `-u 3` and poly(A) trimming. For PE data, R2 uses a fixed 14-nt 5' poly(T) trim and removes `NNNAGATCGGAAGAGC` from the 3' end. Reads shorter than 15 nt are discarded, Q20 trimming is applied to both mates, no UMI deduplication is performed, and circRNA Steps 8-10 operate on no-dedup reads. | + +The SMARTer smRNA profile follows Takara Bio's official [supplemental product information and data-processing guidance](https://catalog.takara-bio.co.jp/PDFS/SMARTer_smRNA-Seq_Kit_for_Illumina.pdf). Takara's [read-trimming FAQ](https://help.takarabio.com/36168/kb/article/109694/do-i-have-to-perform-any-nucleotide-trimming-of-the-smarter-smrnasequencing-reads-prior-to-mapping) explains the R1 5' three-nucleotide and 3' poly(A) trimming requirements, and its [UMI FAQ](https://help.takarabio.com/36168/kb/article/109681) confirms that SMARTer smRNA libraries do not incorporate UMIs. The PE extension uses the observed complementary R2 structure: 14 technical 5' bases, reverse-complement insert, three template-switch-derived bases, then the Illumina adapter. +Because this protocol introduces an artificial poly(A) tract, EMapper's biological poly(A)-tail detection is disabled under `smarter_smrna`. + ## Pipeline Steps | Step | Description | |------|-------------| | 1 | Raw FASTQ quality control using FastQC | -| 2 | UMI motif analysis and ACC motif fraction calculation | -| 3 | UMI extraction, adapter trimming, and read-through UMI removal | +| 2 | Kit-aware library profiling: UMI motif/ACC analysis for `smarter_total_rna`, or a no-UMI library profile for `smarter_smrna` | +| 3 | Kit-aware preprocessing: existing UMI extraction, adapter trimming, and read-through removal for Total RNA; official 5'/poly(A) trimming followed by Q20 trimming for smRNA | | 4 | Quality control of trimmed FASTQs | | 5 | Bacterial contamination screening (E. coli, Mycoplasma) using BBSplit | -| 6 | Two-pass STAR alignment with UMI deduplication | +| 6 | Kit-aware STAR alignment: Total RNA retains two-pass/splice/chimeric alignment and UMI deduplication; smRNA uses end-to-end, unspliced short-read alignment with absolute match thresholds and no UMI deduplication | | 7 | Library strandedness detection; splice/kb complementary gDNA-contamination QC proxy (primary source: STAR Log.final.out from Step 6) | | 8 | CIRCexplorer2-based circular RNA detection | | 9 | CIRI2-based circular RNA detection using BWA alignments | @@ -242,7 +277,7 @@ bash EVscope.sh --sample_name Example_Data \ | 24 | Comprehensive quality control summary generation | | 25 | Coverage analysis and bigWig generation using EMapper | | 26 | Coverage density plots for RNA types and meta-gene regions | -| 27 | Final interactive HTML report generation | +| 27 | Final interactive HTML report plus a separate concise RNA-seq summary report | **splice/kb interpretation.** EVscope reports splice junction crossings per kilobase of uniquely mapped aligned sequence as a complementary genomic-DNA-contamination QC proxy. For paired-end STAR logs, the denominator uses both mates' uniquely mapped aligned bases per fragment/read-pair, not insert size: `splice/kb = Number of splices: Total x 1000 / (Uniquely mapped reads number x Average mapped length)`. In an in-house EV RNA-seq pilot using EXODUS-M isolation, miRNeasy Advanced extraction and 400 uL plasma input (N=5), the single no-DNase control had splice/kb = 0.13; the TURBO DNase pilot libraries had splice/kb values of 1.67 (0.5 U, 10 min), 0.63 (0.25 U, 10 min), 1.64 (0.5 U, 5 min), and 0.94 (0.25 U, 5 min). Values closer to the no-DNase pilot value may indicate higher residual genomic-DNA-derived contribution, whereas higher values are directionally compatible with increased spliced/transcript-derived signal after DNase treatment. These empirical values are exploratory reference guides, not universal cutoffs and not a direct DNA quantification assay. @@ -253,12 +288,12 @@ Each sample generates an output directory with the following structure: ``` _EVscope_output/ ├── Step_01_Raw_QC/ # FastQC reports for raw reads -├── Step_02_UMI_Analysis/ # UMI motif analysis and ACC fraction -├── Step_03_UMI_Adaptor_Trim/ # UMI extraction, adapter trimming, clean FASTQs +├── Step_02_UMI_Analysis/ # Kit-aware library profile; UMI/ACC only for Total RNA +├── Step_03_UMI_Adaptor_Trim/ # Kit-aware trimming and clean FASTQs ├── Step_04_Trimmed_QC/ # FastQC reports for trimmed reads ├── Step_05_Bacterial_Filter/ # BBSplit bacterial/mycoplasma screening -├── Step_06_Alignment_Initial/ # Initial STAR alignment and UMI-deduplicated FASTQs -├── Step_06_Alignment_Refined/ # Refined STAR alignment and UMI-deduplicated BAM +├── Step_06_Alignment_Initial/ # Initial STAR alignment and kit-appropriate duplicate handling +├── Step_06_Alignment_Refined/ # Refined STAR alignment and final BAM ├── Step_07_Strand_Detection/ # Strandedness and splice/kb metrics ├── Step_08_CIRCexplorer2_circRNA/ # CIRCexplorer2 circRNA results ├── Step_09_CIRI2_circRNA/ # CIRI2 circRNA results @@ -277,10 +312,14 @@ Each sample generates an output directory with the following structure: ├── Step_24_MultiQC_Summary/ # QC summary matrix ├── Step_25_EMapper_BigWig_Quantification/ # EMapper coverage and bigWig ├── Step_26_BigWig_Density_Plot/ # RNA type and meta-gene density plots -├── Step_27_HTML_Report/ # Interactive HTML report +├── Step_27_HTML_Report/ # Detailed and standalone concise HTML reports └── EVscope_pipeline.log # Pipeline execution log ``` +The historical Step 2, Step 3, and Step 6 directory names are retained for output-path compatibility. Under `smarter_smrna`, Step 2 contains the smRNA library profile, Step 3 contains the trimmed/clean R1 and, for PE input, R2, while Step 6 does not perform UMI extraction or UMI-based deduplication. UMI/ACC and read-through-specific files are produced only by `smarter_total_rna`. + +For `smarter_smrna`, STAR is run with end-to-end, unspliced alignment (`--alignEndsType EndToEnd`, `--alignIntronMax 1`) and at most one mismatch. The normalized score/match filters are disabled and replaced with absolute floors derived from the 15-nt post-trimming minimum: 12 alignment-score and 14 matched bases per mate (24/28 for PE). Equally best multi-mapping loci are retained, up to the existing limit of 100, because many short ncRNAs are not uniquely mappable. Two-pass alignment remains disabled, while chimeric-junction output is enabled for the circRNA workflow. + ## Troubleshooting - **Dependency Not Found**: Verify Conda environments with `conda list -n evscope_env`. @@ -290,14 +329,26 @@ Each sample generates an output directory with the following structure: ## FAQ -**Q: Can EVscope process non-SMARTer-seq data?** -A: Yes, modify UMI parameters in `bin/Step_02_*.py` and `bin/Step_03_UMIAdapterTrimR1.py`. +**Q: Which library kits are supported?** +A: Use `--kit smarter_total_rna` for the SMARTer Stranded Total RNA-Seq Kit v3 (the default) or `--kit smarter_smrna` for the SMARTer smRNA-Seq Kit for Illumina. Other library protocols are not represented by these profiles and require separately validated preprocessing. + +**Q: Can `smarter_smrna` process paired-end FASTQs?** +A: Yes. Provide R1 first and R2 second. R1 is trimmed with the official `-u 3`/poly(A) profile. R2 is trimmed by 14 nt at its 5' poly(T) end and by the `NNNAGATCGGAAGAGC` structure at its 3' end. Both mates remain synchronized and are passed to paired-end alignment without UMI deduplication. + +**Q: How does EVscope preprocess SMARTer smRNA reads?** +A: R1 is treated as the sense read. Step 2 records a no-UMI library profile. Step 3 applies `cutadapt -m 15 -u 3 -a AAAAAAAAAA`, followed by the existing Q20 quality trimming. It does not run UMI extraction, UMI read-through removal, or UMI-based deduplication. Steps 8-10 can detect circRNAs from the no-dedup reads, although calls from this 3'-polyadenylation protocol should be interpreted cautiously. + +**Q: Do I need to specify `--strand forward` for SMARTer smRNA?** +A: No. If `--strand` is omitted with `--kit smarter_smrna`, EVscope selects `forward`. An explicitly supplied `--strand` value is honored. + +**Q: Can I resume a v1.0 output directory with EVscope v1.1?** +A: No. v1.1 includes the kit and pipeline version in step signatures. Use `--force` to regenerate the selected steps, or remove the old output directory and start a fresh run. **Q: How do I run specific pipeline steps?** A: Use `--run_steps`, e.g., `--run_steps 1,3,5-8`. **Q: How do I view the final report?** -A: Open `Step_27_HTML_Report/_final_report.html` in a web browser. +A: Open `Step_27_HTML_Report/_final_report.html` for the comprehensive report, or `_summary_report.html` in the same directory for the standalone one-paragraph summary. ## Contributing diff --git a/bin/Step_07_bam2strand.py b/bin/Step_07_bam2strand.py index ae7dc5a..d0f2407 100755 --- a/bin/Step_07_bam2strand.py +++ b/bin/Step_07_bam2strand.py @@ -84,8 +84,10 @@ def clean_sample_name_from_bam_basename(base_name): """Convert EVscope BAM basename to display sample name for plotting.""" suffixes = [ '_STAR_umi_dedup_Aligned.sortedByCoord.out', + '_STAR_no_dedup_Aligned.sortedByCoord.out', '_Aligned.sortedByCoord.out', '_Aligned.sortedByCoord_umi_dedup.out', + '_Aligned.sortedByCoord_no_dedup.out', ] sample_name = base_name for suffix in suffixes: diff --git a/bin/Step_15_combine_total_RNA_expr_matrix.py b/bin/Step_15_combine_total_RNA_expr_matrix.py index 330b799..acd515e 100755 --- a/bin/Step_15_combine_total_RNA_expr_matrix.py +++ b/bin/Step_15_combine_total_RNA_expr_matrix.py @@ -18,7 +18,10 @@ def parse_args() -> argparse.Namespace: """Parse command-line arguments.""" parser = argparse.ArgumentParser(description="Merge gene and circRNA expression tables.") parser.add_argument("--gene_expr", required=True) - parser.add_argument("--circRNA_expr", required=True) + parser.add_argument( + "--circRNA_expr", + help="Optional circRNA expression matrix. Omit it to normalize a gene-only matrix.", + ) parser.add_argument("--out_matrix", required=True) return parser.parse_args() @@ -103,7 +106,7 @@ def write_merged_matrix(rows: Iterable[Dict[str, str]], output_path: str) -> Non def main() -> None: args = parse_args() gene_data = read_gene_data(args.gene_expr) - circ_data = read_circ_data(args.circRNA_expr) + circ_data = read_circ_data(args.circRNA_expr) if args.circRNA_expr else [] write_merged_matrix(gene_data + circ_data, args.out_matrix) diff --git a/bin/Step_24_generate_QC_matrix.py b/bin/Step_24_generate_QC_matrix.py index 6b3bdcb..81d6387 100755 --- a/bin/Step_24_generate_QC_matrix.py +++ b/bin/Step_24_generate_QC_matrix.py @@ -73,11 +73,17 @@ def parse_args(): sys.exit("Error: No input parameters provided.") parser.add_argument("--raw_fastqc_zips", nargs='+', required=False, help="List of raw FastQC zip files from Step 1.") + parser.add_argument( + "--kit", + choices=("smarter_total_rna", "smarter_smrna"), + default="smarter_total_rna", + help="Library kit profile used for kit-aware QC labels (default: smarter_total_rna).", + ) parser.add_argument("--trimmed_fastqs", nargs='+', required=False, help="List of trimmed FASTQ files from Step 3.") parser.add_argument("--ecoli_fastqs", nargs='*', required=False, help="List of E.coli FASTQ files from BBSplit (Step 5). Empty list means Step 5 ran and found 0 reads.") parser.add_argument("--myco_fastqs", nargs='*', required=False, help="List of Mycoplasma FASTQ files from BBSplit (Step 5). Empty list means Step 5 ran and found 0 reads.") parser.add_argument("--ribo_fastqs", nargs='*', required=False, help="List of rRNA FASTQ files from RiboDetector (Step 23). Empty list means Step 23 ran and found 0 reads.") - parser.add_argument("--ACC_motif_fraction", required=False, help="ACC motif fraction TSV file from Step 2.") + parser.add_argument("--ACC_motif_fraction", required=False, help="ACC motif fraction TSV file from Step 2 (SMARTer Total RNA only).") parser.add_argument("--kraken_report", required=False, help="Kraken2 report file from Step 19.") parser.add_argument("--bam2strand_file", required=False, help="bam2strandness output file from Step 7.") parser.add_argument("--length_strandness_file", required=False, help="Read-length-stratified strandness TSV from Step 7.") @@ -358,7 +364,7 @@ def process_featureCounts_file(file_path): def main(): args = parse_args() - ordered_metrics = [] + ordered_metrics = [("Library Kit", args.kit)] total_raw_reads, raw_gcs, raw_lengths = 0, [], [] if args.raw_fastqc_zips: @@ -377,7 +383,9 @@ def main(): avg_raw_gc = mean(raw_gcs) if raw_gcs else 0 avg_raw_read_length = mean(raw_lengths) if raw_lengths else 0 - acc_fraction = safe_call(extract_acc_motif_fraction, args.ACC_motif_fraction, default="NA") + acc_fraction = "NA" + if args.kit == "smarter_total_rna": + acc_fraction = safe_call(extract_acc_motif_fraction, args.ACC_motif_fraction, default="NA") total_trimmed_reads = 0 avg_trimmed_read_length = 0 @@ -395,7 +403,7 @@ def main(): ordered_metrics.append(("Total Raw Read Pairs (Fragments)", int(total_raw_read_pairs) if total_raw_read_pairs > 0 else "NA")) ordered_metrics.append(("%GC of Raw Reads", round(avg_raw_gc, 2) if avg_raw_gc > 0 else "NA")) ordered_metrics.append(("Average Raw Read Length", int(round(avg_raw_read_length)) if avg_raw_read_length > 0 else "NA")) - ordered_metrics.append(("ACC motif fraction from UMI region", acc_fraction)) + ordered_metrics.append(("ACC motif fraction from UMI region (SMARTer Total RNA only)", acc_fraction)) ordered_metrics.append(("Total Trimmed Reads (R1+R2)", int(total_trimmed_reads) if total_trimmed_reads > 0 else "NA")) ordered_metrics.append(("Average Trimmed Read Length", int(round(avg_trimmed_read_length)) if avg_trimmed_read_length > 0 else "NA")) ordered_metrics.append(("Percentage of Reads Remaining after Trimming", percent_reads_after_trimming if total_raw_reads > 0 else "NA")) @@ -426,7 +434,9 @@ def main(): star_metrics = safe_call(parse_star_log, args.STAR_log, default={}) star_input_reads = int(star_metrics.get("Number of input reads (STAR)", 0)) - perc_reads_after_dedup = round((star_input_reads / initial_star_mapped_fragments) * 100, 2) if initial_star_mapped_fragments > 0 else 0 + perc_reads_after_dedup = "NA" + if args.kit == "smarter_total_rna" and initial_star_mapped_fragments > 0: + perc_reads_after_dedup = round((star_input_reads / initial_star_mapped_fragments) * 100, 2) star_unique_reads = int(star_metrics.get("Uniquely mapped reads number (STAR)", 0)) perc_unique_vs_star_input = round((star_unique_reads / star_input_reads) * 100, 2) if star_input_reads > 0 else 0 @@ -447,10 +457,20 @@ def main(): splice_per_kb_source = sanitize_splice_per_kb_source(splice_per_kb_source, args.STAR_log) splice_available = splice_per_kb != "NA" + if args.kit == "smarter_total_rna": + refined_input_metrics = [ + ("Number of Reads after UMI-deduplication (STAR Input)", star_input_reads if star_input_reads > 0 else "NA"), + ("Percentage of UMI-dedup Fragments", perc_reads_after_dedup), + ] + else: + refined_input_metrics = [ + ("Number of Reads in Refined STAR Input (No UMI Deduplication)", star_input_reads if star_input_reads > 0 else "NA"), + ("Percentage of UMI-dedup Fragments (Not Applicable: No UMI Deduplication)", "NA"), + ] + ordered_metrics.extend([ ("Total Fragments Mapped to Human (First STAR)", int(initial_star_mapped_fragments) if initial_star_mapped_fragments > 0 else "NA"), - ("Number of Reads after UMI-deduplication (STAR Input)", star_input_reads if star_input_reads > 0 else "NA"), - ("Percentage of UMI-dedup Fragments", perc_reads_after_dedup if perc_reads_after_dedup > 0 else "NA"), + *refined_input_metrics, ("Average Input Read Length (STAR)", float(star_metrics.get("Average input read length (STAR)", 0))), ("Average Mapped Length per Unique Fragment (STAR)", float(star_metrics.get("Average mapped length (STAR)", 0))), ("Percentage of Mapped Reads on Forward Strand", value_or_na(forward_strand, strand_available)), @@ -550,4 +570,4 @@ def main(): print(f"QC matrix successfully generated: {args.output}") if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/bin/Step_27_generate_summary_report.py b/bin/Step_27_generate_summary_report.py new file mode 100644 index 0000000..8a46b0a --- /dev/null +++ b/bin/Step_27_generate_summary_report.py @@ -0,0 +1,303 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: MIT +"""Generate a concise, standalone RNA-seq summary report for one EVscope sample.""" + +from __future__ import annotations + +import argparse +import csv +import html +import math +from datetime import datetime +from pathlib import Path +from typing import Dict, Iterable, List, Optional, Tuple + + +NA_VALUES = {"", "na", "n/a", "nan", "none", "null"} + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Generate a one-paragraph standalone EVscope RNA-seq summary report." + ) + parser.add_argument("--sample-name", required=True) + parser.add_argument("--qc-summary", help="Step 24 two-column QC summary TSV.") + parser.add_argument("--expression-matrix", help="Gene expression TSV used for the top-10 list.") + parser.add_argument( + "--kit", + choices=("smarter_total_rna", "smarter_smrna"), + default="smarter_total_rna", + ) + parser.add_argument("--strand", choices=("forward", "reverse", "unstrand"), default="unstrand") + parser.add_argument("--paired-end", choices=("true", "false"), required=True) + parser.add_argument("--output", required=True) + return parser.parse_args() + + +def read_qc_metrics(path: Optional[str]) -> Dict[str, str]: + if not path or not Path(path).is_file(): + return {} + metrics: Dict[str, str] = {} + with open(path, newline="", encoding="utf-8-sig") as handle: + reader = csv.DictReader(handle, delimiter="\t") + if not reader.fieldnames or not {"Metric", "Value"}.issubset(reader.fieldnames): + return {} + for row in reader: + metric = (row.get("Metric") or "").strip() + if metric: + metrics[metric] = (row.get("Value") or "").strip() + return metrics + + +def number(value: object) -> Optional[float]: + if value is None: + return None + text = str(value).strip() + if text.lower() in NA_VALUES: + return None + try: + result = float(text.replace(",", "").rstrip("%")) + except ValueError: + return None + return result if math.isfinite(result) else None + + +def metric_number(metrics: Dict[str, str], name: str) -> Optional[float]: + return number(metrics.get(name)) + + +def format_count(value: Optional[float]) -> str: + if value is None: + return "not available" + absolute = abs(value) + for scale, label in ((1_000_000_000, "billion"), (1_000_000, "million")): + if absolute >= scale: + shown = f"{value / scale:.2f}".rstrip("0").rstrip(".") + return f"{shown} {label}" + if float(value).is_integer(): + return f"{int(value):,}" + return f"{value:,.1f}" + + +def format_pct(value: Optional[float]) -> str: + return "not available" if value is None else f"{value:.2f}%" + + +def fragment_count(read_count: Optional[float], paired_end: bool) -> Optional[float]: + if read_count is None: + return None + return read_count / 2.0 if paired_end else read_count + + +def join_phrases(items: List[str]) -> str: + if len(items) < 2: + return "".join(items) + if len(items) == 2: + return " and ".join(items) + return ", ".join(items[:-1]) + ", and " + items[-1] + + +def first_present(row: Dict[str, str], names: Iterable[str]) -> str: + for name in names: + value = (row.get(name) or "").strip() + if value and value.lower() not in NA_VALUES: + return value + return "" + + +def top_genes(path: Optional[str], limit: int = 10) -> List[Tuple[str, float]]: + if not path or not Path(path).is_file(): + return [] + genes: Dict[str, float] = {} + with open(path, newline="", encoding="utf-8-sig") as handle: + reader = csv.DictReader(handle, delimiter="\t") + for row in reader: + gene_type = first_present(row, ("GeneType", "gene_type", "Biotype")).lower() + if gene_type in {"circrnas", "circrna", "artifact"}: + continue + gene = first_present(row, ("GeneSymbol", "GeneName", "gene_name", "GeneID", "gene_id")) + read_count = number(first_present(row, ("ReadCounts", "read_count", "expected_count"))) + if not gene or read_count is None: + continue + genes[gene] = genes.get(gene, 0.0) + read_count + return sorted(genes.items(), key=lambda item: (-item[1], item[0]))[:limit] + + +def expected_orientation(kit: str, strand: str) -> Tuple[str, str]: + if strand in {"forward", "reverse"}: + expected = strand + else: + expected = "reverse" if kit == "smarter_total_rna" else "forward" + if kit == "smarter_total_rna" and expected == "reverse": + label = "expected reverse/dUTP orientation" + elif kit == "smarter_smrna" and expected == "forward": + label = "expected forward/sense orientation" + else: + label = f"configured {expected} orientation" + return expected, label + + +def build_summary( + sample_name: str, + metrics: Dict[str, str], + genes: List[Tuple[str, float]], + kit: str, + strand: str, + paired_end: bool, +) -> str: + unit = "read pairs" if paired_end else "reads" + complete_unit = "complete read pairs" if paired_end else "reads" + clauses: List[str] = [] + + raw = metric_number( + metrics, + "Total Raw Read Pairs (Fragments)" if paired_end else "Total Raw Reads (R1+R2)", + ) + trimmed = fragment_count(metric_number(metrics, "Total Trimmed Reads (R1+R2)"), paired_end) + if raw is not None: + opening = f"Sample {sample_name} initially yielded {format_count(raw)} {unit}" + if trimmed is not None: + opening += f", with {format_count(trimmed)} {unit} remaining after trimming" + clauses.append(opening + ".") + elif trimmed is not None: + clauses.append(f"After trimming, {format_count(trimmed)} {unit} remained for sample {sample_name}.") + + refined_metric = ( + "Number of Reads after UMI-deduplication (STAR Input)" + if kit == "smarter_total_rna" + else "Number of Reads in Refined STAR Input (No UMI Deduplication)" + ) + refined = metric_number(metrics, refined_metric) + unique = metric_number(metrics, "Number of Uniquely Mapped Reads (STAR)") + multi = metric_number(metrics, "Number of Multi-mapped Reads (STAR)") + mapped = None if unique is None and multi is None else (unique or 0.0) + (multi or 0.0) + mapped_pct = mapped * 100.0 / refined if mapped is not None and refined and refined > 0 else None + if refined is not None: + if kit == "smarter_total_rna": + alignment = ( + f"After UMI deduplication and orphan-mate removal, {format_count(refined)} " + f"{complete_unit} remained" + ) + else: + alignment = f"No UMI deduplication was applied; {format_count(refined)} {unit} entered refined alignment" + if mapped is not None: + alignment += f", of which {format_count(mapped)} ({format_pct(mapped_pct)}) mapped to hg38" + clauses.append(alignment + ".") + elif mapped is not None: + clauses.append(f"The refined STAR alignment mapped {format_count(mapped)} {unit} to hg38.") + + if mapped and mapped > 0: + unique_pct = (unique or 0.0) * 100.0 / mapped + multi_pct = (multi or 0.0) * 100.0 / mapped + clauses.append( + f"Among human-mapped {unit}, {format_pct(unique_pct)} were unique and " + f"{format_pct(multi_pct)} were multi-mapped." + ) + + contamination_parts: List[str] = [] + for organism, count_metric, pct_metric in ( + ( + "E. coli", + "Number of Reads Mapped to Escherichia coli (BBSplit)", + "Percentage of Trimmed Reads Mapped to E. coli (BBSplit)", + ), + ( + "Mycoplasma", + "Number of Reads Mapped to Mycoplasma (BBSplit)", + "Percentage of Trimmed Reads Mapped to Mycoplasma (BBSplit)", + ), + ): + count = fragment_count(metric_number(metrics, count_metric), paired_end) + pct = metric_number(metrics, pct_metric) + if count is not None: + pct_text = f" ({format_pct(pct)} of trimmed {unit})" if pct is not None else "" + contamination_parts.append(f"{format_count(count)} {unit}{pct_text} to {organism}") + if contamination_parts: + clauses.append("BBSplit assigned " + join_phrases(contamination_parts) + ".") + + region_parts: List[str] = [] + for adjective, region in (("intronic", "Intronic"), ("exonic", "Exonic"), ("intergenic", "Intergenic")): + count = metric_number(metrics, f"Number of Reads Mapped to {region} Regions") + pct = metric_number(metrics, f"Percentage of Reads Mapped to {region} Regions") + if count is not None: + pct_text = f" ({format_pct(pct)})" if pct is not None else "" + region_parts.append(f"{format_count(count)} {unit}{pct_text} as {adjective}") + if region_parts: + clauses.append(f"Region assignment counted {join_phrases(region_parts)}.") + + forward = metric_number(metrics, "Percentage of Mapped Reads on Forward Strand") + reverse = metric_number(metrics, "Percentage of Mapped Reads on Reverse Strand") + informative = (forward or 0.0) + (reverse or 0.0) + expected, orientation_label = expected_orientation(kit, strand) + expected_value = forward if expected == "forward" else reverse + if expected_value is not None and informative > 0: + expected_pct = expected_value * 100.0 / informative + clauses.append( + f"{format_pct(expected_pct)} of strand-informative alignments matched the {orientation_label}." + ) + + if genes: + gene_text = ", ".join(f"{gene} ({format_count(count)})" for gene, count in genes) + clauses.append(f"The top {len(genes)} genes by read count were {gene_text}.") + else: + clauses.append("A top-gene list was not available because no compatible expression matrix was found.") + + if not clauses: + return f"Summary metrics were not available for sample {sample_name}." + return " ".join(clauses) + + +def render_html(sample_name: str, summary: str) -> str: + safe_sample = html.escape(sample_name) + safe_summary = html.escape(summary) + generated = html.escape(datetime.now().astimezone().strftime("%d %B %Y, %H:%M %Z")) + return f""" + + + + + {safe_sample} — EVscope RNA-seq Summary + + + +

+

EVscope RNA-seq Summary

+
Sample: {safe_sample}
+

{safe_summary}

+
Generated by EVscope on {generated}
+
+ + +""" + + +def main() -> None: + args = parse_args() + metrics = read_qc_metrics(args.qc_summary) + genes = top_genes(args.expression_matrix) + summary = build_summary( + sample_name=args.sample_name, + metrics=metrics, + genes=genes, + kit=args.kit, + strand=args.strand, + paired_end=args.paired_end == "true", + ) + output = Path(args.output) + output.parent.mkdir(parents=True, exist_ok=True) + output.write_text(render_html(args.sample_name, summary), encoding="utf-8") + print(f"Concise summary report generated: {output}") + + +if __name__ == "__main__": + main() diff --git a/bin/Step_27_html_report.Rmd b/bin/Step_27_html_report.Rmd index 71a4a49..dd68a9f 100755 --- a/bin/Step_27_html_report.Rmd +++ b/bin/Step_27_html_report.Rmd @@ -1,5 +1,5 @@ --- -title: "EVscope Comprehensive Analysis Report of Total RNAseqs" +title: "EVscope Comprehensive RNA-seq Analysis Report" license: "MIT" author: - Rambo (Yiyong) Zhao @@ -417,11 +417,12 @@ iframe { /* --- Data Tables Container --- */ .dataTables_wrapper { background: var(--card-bg); + width: 90%; padding: 1.5rem; border-radius: var(--radius-md); box-shadow: var(--shadow-sm); border: 1px solid var(--border-color); - margin: 1rem 0; + margin: 1rem auto; } /* --- DataTables Styling --- */ @@ -753,6 +754,19 @@ if (!dir.exists(absolute_output_dir)) { stop("FATAL: Output directory not found. EVscope_OUTPUT_DIR was not set or is invalid.") } +# --- Library Kit Profile --- +evscope_kit <- Sys.getenv("EVscope_KIT", unset = "smarter_total_rna") +if (!evscope_kit %in% c("smarter_total_rna", "smarter_smrna")) { + warning("Unknown EVscope_KIT value '", evscope_kit, "'; using smarter_total_rna.") + evscope_kit <- "smarter_total_rna" +} +is_smarter_smrna <- identical(evscope_kit, "smarter_smrna") +kit_display_name <- if (is_smarter_smrna) { + "SMARTer smRNA-seq Kit" +} else { + "SMARTer Stranded Total RNA-Seq Kit" +} + # --- Retrieve Figures Directory --- abs_figures_dir <- Sys.getenv("EVscope_FIGURES_DIR") @@ -1037,7 +1051,9 @@ display_all_images <- function(subdir_path, pattern, main_caption) {
-**EVscope** is an open-source bioinformatics pipeline tailored for EV-enriched total RNA-seq datasets. EVscope includes EMapper, a two-round expectation-maximization (EM) module: Round 1 resolves positional ambiguity using alignment scores and local coverage evidence, and Round 2 resolves gene-level ambiguity for overlapping annotations using abundance-based weights with effective-length normalization. EMapper's primary novelty is not a claim of superior conventional gene readcount quantification over tools such as RSEM; instead, it generates EM-weighted genome-coordinate BigWig coverage tracks at single-base or binned resolution, with strand-aware outputs and RNA annotation support for downstream genome-browser inspection, RNA-biotype density plots, and meta-gene coverage profiling. The pipeline systematically integrates 27 major steps, including quality control, library-structure analysis, contamination assessment, read alignment, strandedness detection, UMI-based deduplication, feature/RNA annotation, expression matrix generation, optional genomic DNA (gDNA) sensitivity summaries, exploratory source-enrichment analysis, and visualization in a comprehensive HTML report. EVscope incorporates a comprehensive, updated annotation covering 20 RNA biotypes, including protein-coding genes, lncRNAs, miRNAs, piRNAs, retrotransposons (LINEs, SINEs, ERVs), and additional non-coding RNAs (tRNAs, rRNAs, snoRNAs). It also includes complementary circRNA detection modules. EVscope is implemented as a convenient, single-command Bash pipeline leveraging Conda-managed standard software packages and custom scripts, supporting reproducible deployment and transparent QC reporting. +**EVscope** is an open-source, kit-aware bioinformatics pipeline for EV-enriched RNA sequencing generated with either the SMARTer Stranded Total RNA-Seq Kit or the SMARTer smRNA-seq Kit. The active library profile is selected through EVscope_KIT. Kit-aware preprocessing applies the appropriate library-structure and adapter handling: UMI extraction, UMI read-through trimming, and UMI-based deduplication are used only for the Total RNA profile, while the smRNA profile follows a no-UMI-deduplication workflow. EVscope includes EMapper, a two-round expectation-maximization (EM) module: Round 1 resolves positional ambiguity using alignment scores and local coverage evidence, and Round 2 resolves gene-level ambiguity for overlapping annotations using abundance-based weights with effective-length normalization. EMapper generates EM-weighted genome-coordinate BigWig coverage tracks at single-base or binned resolution, with strand-aware outputs and RNA annotation support for downstream genome-browser inspection, RNA-biotype density plots, and meta-gene coverage profiling. The workflow integrates quality control, library-structure analysis, contamination assessment, read alignment, feature/RNA annotation, expression matrix generation, optional genomic DNA sensitivity summaries, exploratory source-enrichment analysis, and visualization in a comprehensive HTML report. + +**Selected library kit:** `r kit_display_name`
@@ -1069,22 +1085,28 @@ if (check_step_dir("Step_01_Raw_QC")) { ``` -# Step 02: UMI Analysis +# Step 02: Library Structure Analysis ```{r step2, results='asis'} if (check_step_dir("Step_02_UMI_Analysis")) { - display_image("Step_02_UMI_Analysis", "motif_logo\\.png$", "UMI Motif Logo") - display_datatable("Step_02_UMI_Analysis", "_ACC_motif_fraction\\.tsv$", "ACC Motif Fraction") + if (is_smarter_smrna) { + display_datatable("Step_02_UMI_Analysis", "_library_profile\\.tsv$", "SMARTer smRNA-seq Library Profile") + } else { + display_image("Step_02_UMI_Analysis", "motif_logo\\.png$", "UMI Motif Logo") + display_datatable("Step_02_UMI_Analysis", "_ACC_motif_fraction\\.tsv$", "ACC Motif Fraction") + } } ``` -# Step 03: UMI & Adaptor Trimming +# Step 03: Library Preprocessing and Adapter Trimming ```{r step3, results='asis'} if (check_step_dir("Step_03_UMI_Adaptor_Trim")) { display_image("Step_03_UMI_Adaptor_Trim", "_read_length_distribution\\.png$", "Read Length Distribution After Trimming") - display_image("Step_03_UMI_Adaptor_Trim", "_R1_readthrough_UMI_trimming_pie\\.png$", "Read Categories Distribution (No UMI Readthrough, UMI Readthrough, Too Short)") + if (!is_smarter_smrna) { + display_image("Step_03_UMI_Adaptor_Trim", "_R1_readthrough_UMI_trimming_pie\\.png$", "Read Categories Distribution (No UMI Readthrough, UMI Readthrough, Too Short)") + } } ``` @@ -1112,7 +1134,11 @@ if (check_step_dir("Step_06_Alignment_Initial")) { ```{r step6b, results='asis'} if (check_step_dir("Step_06_Alignment_Refined")) { - parse_and_display_star_log("Step_06_Alignment_Refined", "_STAR_umi_dedup_Log\\.final\\.out$", "STAR Refined Alignment Summary (Post-Deduplication)") + if (is_smarter_smrna) { + parse_and_display_star_log("Step_06_Alignment_Refined", "_STAR_no_dedup_Log\\.final\\.out$", "STAR Refined Alignment Summary (No UMI Deduplication)") + } else { + parse_and_display_star_log("Step_06_Alignment_Refined", "_STAR_umi_dedup_Log\\.final\\.out$", "STAR Refined Alignment Summary (Post-UMI-Deduplication)") + } } ``` @@ -1373,4 +1399,4 @@ sessionInfo()
EVscope Comprehensive Analysis Report Generated with R Markdown - Self-contained HTML document -
\ No newline at end of file + diff --git a/repro/README.md b/repro/README.md index 3d7ffd8..ef9724c 100644 --- a/repro/README.md +++ b/repro/README.md @@ -11,6 +11,6 @@ bash tests/smoke/run_smoke.sh sha256sum -c repro/checksums.sha256 ``` -The smoke test covers syntax, CLI help/version, toy FASTQ processing, read-length plotting when matplotlib is available, and a tiny synthetic EMapper/BigWig/bigWig2Expression run when `pysam`, `pyBigWig`, `numba`, and `numpy` are available. +The smoke test covers syntax, kit-aware CLI validation, SMARTer smRNA structure trimming/no-UMI alignment with synthetic reads, STAR argument injection from config, cross-kit resume protection, standalone concise-summary report logic, read-length plotting when matplotlib is available, and a tiny synthetic EMapper/BigWig/bigWig2Expression run when `pysam`, `pyBigWig`, `numba`, and `numpy` are available. -`repro/checksums.sha256` covers tracked release-critical files except itself. Git commit/tag identity remains the authoritative source-tree integrity record. +`repro/checksums.sha256` covers tracked release-critical files except itself. Text files use LF line endings through `.gitattributes` so hashes match across platforms. Git commit/tag identity remains the authoritative source-tree integrity record. diff --git a/repro/checksums.sha256 b/repro/checksums.sha256 index b060cef..20854aa 100644 --- a/repro/checksums.sha256 +++ b/repro/checksums.sha256 @@ -1,19 +1,20 @@ +2934519b0f064594c493c3f47c6e42add62d07126547c265f7eede217dc48c62 .gitattributes 58eb25da24ec299883e2c1b3a60a972476fe3d0c7db4aef00c50966c6dfaa09e .github/workflows/smoke.yml -866eb9228dfc71cfa1cf3aff07ab1beed7ba1ed827f86aa73791bf49eae1da11 .gitignore -bf94ca778d2d1c0295cddcb8a451975d6b4fadbca784cabee1ee0be3b09d956b EVscope.conf -0e3c0bb85dc8303c9cce36bece86ff406c3a61fc11f66e1f218efa921d65b0ed EVscope.sh +0357489076a46a50bc6c32bed3d61c82a2c57205efdf45b57265160eca2adf9f .gitignore +33ba6396cca8ff2301be81b394c64bbf45327b14153c74b6a4bb57787d37b6db EVscope.conf +d990fbfef5f7737e998ebc5f53600e1f3378f2186bda460dee2facb8c14207e6 EVscope.sh 8003003ded85bbc27a43f9719d604a358993818007a6dc7a35eb45e06c749e15 LICENSE -364b98b069a67c435d7ec3d216e48bf57dfb8506f2a747a0dcccfdfb8e7f6188 README.md +1395372d2fef6dc6ae2edc3fd0487924ca36f78ad848dd594dfe76632dc0ac8e README.md a4384863c9f164f3d1a06e782beb766971f4e52da2a9d724a8669a3303bc687f bin/Step_02_calculate_ACC_motif_fraction.py 248f1e3469ff052ab34f191bb2982a2f9b0194fc3cfb06ac458ab931c08a6b7c bin/Step_02_plot_fastq2UMI_motif.py 2c22aad993cce66020d295c35a2a9bba53da89a537a50b575f591e03d56ed3ee bin/Step_03_UMIAdapterTrimR1.py 2a533641c2468d7236712b54dc4dfa57910037ff48ce6eefc66505b9042fb034 bin/Step_03_plot_fastq_read_length_dist.py -53fccf8a0d3d339e3b9fda6959767a1db6ed436de0ed7c3e10baded7442f1cf5 bin/Step_07_bam2strand.py +270775851426458c286b59f31cbbeb74fa8e67633ab442d080e23813e7e57aed bin/Step_07_bam2strand.py b0afe6b993e90c2274ad8987ed1fcad9c3ad1077a13993e13d2d26dddf0c78e1 bin/Step_08_convert_CIRCexplorer2CPM.py 9b62e03d0997f58b5d5e6ea7e0db25f640bdb26089f33045645c7be34f5b6c7e bin/Step_09_convert_CIRI2CPM.py e614cbff8a281e077e75d42ff86422fcdbb9c720275cab1a9c599c3e5075cc91 bin/Step_10_circRNA_merge.py a684163453540a4be692d5c1fde41c9669aa6f22121ffdce50470aa13c24f50d bin/Step_13_gDNA_corrected_featureCounts.py -b397101f601bfdf47aadb097ccf5c181483bfff6357a49b2125408ef5c942d49 bin/Step_15_combine_total_RNA_expr_matrix.py +d8f3cd76775976fc6b4a0b0f7004b9584854abc302bb8a3634ce8b4a8c2d0d80 bin/Step_15_combine_total_RNA_expr_matrix.py d232f038d90e70f94c68e7f53ed9232142d6c9d82ecc4cb20863424717b395c1 bin/Step_15_featureCounts2TPM.py 4ceda668b709f7de6f51b38695130769bf64901deae59e41945b3678c298320f bin/Step_15_plot_RNA_distribution_1subplot.py 4e2487bb3d67e501dd17da90b9598bb4e140d8b9585b7998323129007cb1bdcc bin/Step_15_plot_RNA_distribution_20subplots.py @@ -22,24 +23,35 @@ d232f038d90e70f94c68e7f53ed9232142d6c9d82ecc4cb20863424717b395c1 bin/Step_15_fe a87ad3a3d9f43d70ea8e08a3e68f3c70b009f0c4cf455b3290b6c4b81df8ce42 bin/Step_17_RSEM2expr_matrix.py 227d9fae284ccd69c32d367cfb48b2b05ec46c715644522f3b1113340b851a34 bin/Step_18_plot_reads_mapping_stats.py 2324aabf1128a7bd59254600cab8853a52ef1638acd493bb23dc8ad364a9a6b6 bin/Step_22_run_RNA_deconvolution_ARIC.py -9cd1be5582af3de152c13afcd3fe240cd827afb5f644c8160689b9de3bd2a8f2 bin/Step_24_generate_QC_matrix.py +8168717fc3ddd95b4eac9d566e7ef221ab4ffae682de682e9a0fae570b5a46dd bin/Step_24_generate_QC_matrix.py d13483469a4e47a25104fd7619289ce106d0fcfc25ebb055aed5dbfcd8a52ba2 bin/Step_25_EMapper.py 412b4c1d64f511918031a44896ae7be0f0929ebd25c0261d76e84d540b77ae04 bin/Step_25_bigWig2Expression.py 2b8bb70804ccd527e229c9016b428abf1cca0e5835d4cdd25b0c837352a1ca31 bin/Step_26_density_plot_over_RNA_types.sh 1778752d48d9da3eb8e5f6a99d6c2ecc5b941e686119edebe8835050801c7291 bin/Step_26_density_plot_over_meta_gene.sh -2919740bc39bae68cefa0c1e7079c14407115a1089690e64d77b7fad39aa32eb bin/Step_27_html_report.Rmd -a838bea718359c48a448ef74a874f0d3c5e806e57243f6fa8aae9a6b09e7019b environments/evscope_env.yml +117a2203bc4af999ae793a069275c7c073a9bda049360e891a408f082321d0ea bin/Step_27_generate_summary_report.py +064b4a9a91d2a6c5ba0f45dca2815ff97e24d19ea60fde16f63c54406d1029ec bin/Step_27_html_report.Rmd +9e9d2bee6de40d4d10d73e1f850ad016e35bcf1e1004e675c89879777006883c environments/evscope_env.yml 692b646527aae5dbf832138d4ca4c888a486d30ebc6639147093e8a429ac39f5 environments/kraken2_env.yml 76d6cce96edd915a72d2a41d1016340f9f03173f402edf3772dcf473f09fc078 environments/picard_env.yml 18248c9f29e03c797b339d26280eea968c37f6d764b6a5f8ac9e6a7514fe4a7b figures/EVscope_pipeline.png -4156fbc760db76834fb1ad62d3d9c57c7775f64c6ac2d7efdc9f92dc4006440d repro/README.md +1e65d09e089ff5990606bde7312fe0960e7981aade3cd13446fa11e736017535 repro/README.md 767c5b88b1c8ca7349ba0efb4a71f758ef78bd27d62c4e36cae04b38ae675688 repro/benchmarks/EMapper_MAQC_qPCR/README.md 593f85dc44c8faeb6a5d5c2dcc036bf229e97f52139070259c193a8af753637b repro/figures/README.md -62fa01923420b924e58ac4ac453da12c8b4f03f6cd655f6248d764df026f3d56 repro/manifest.tsv +1733e93f9bcfd90b5922e07627de666ec364cc9cf494e1ccbc72a2fd2351e481 repro/manifest.tsv e37c5a511bfa1d6235abe058cf2941f7ccfdc45fef7d9ff30a56c7387b588434 repro/metadata/software_versions.tsv 61f65dd9598c7421870bf4918332637f83f34abe1f9fc7cc095852ef7a1b44f7 tests/smoke/.gitignore -bf94bd4d4c60a43890acafd0c456ca4d692f009003ba7fee3064887fb39da7ed tests/smoke/README.md +c200023ad0419b2f1f30b63be9a443c46de4a782713f19f34e227f7d4072a68c tests/smoke/README.md +e8dbb54045c744a40b49f1ee069a82703b3fa9acf56f00c19b9c1a809e038412 tests/smoke/data/smarter_smrna_R1.fastq +e835bc8327e9d41dc2fb0a1c741824336644eaaf87186e66f15f7baac748a7bd tests/smoke/data/smarter_smrna_R2.fastq 3bf0f4afb7771bc4916fe8499bb0310beaa9ab377f7d03eb121fc8fd6bad97d4 tests/smoke/data/toy_R2.fastq 4f817aa0b325b273c2e94e57e19aee0e6c0a6e200e7f679687dd0a151c3cc7ce tests/smoke/data/toy_R2.fastq.gz -0a4ad44f7e3244443e7647debc398b462333aca99a0266ccb28147b7b7a49754 tests/smoke/expected/expected_files.txt -2e10439cccafe7348a9c06529f79b33be14839e3c2b3a51774eb666a95d7f681 tests/smoke/run_smoke.sh +ea770e19e9b888df93ec92e7f0945baacc05ead12ab8124b431c5993cd8e04df tests/smoke/expected/expected_files.txt +05dcb4361adcde7000e30fde4de63040222c6348b51072e224e993e010b5533d tests/smoke/expected/smarter_smrna_R2_clean.fastq +7bf9bf3bc51d54589ab699492aff2fcee116662a40840dcc9b1199127261ae30 tests/smoke/expected/smarter_smrna_clean.fastq +9696908a77da99a50f3819a623880b508a45642ad3c3f01a3399d46aad40cb68 tests/smoke/minimal.conf +3f2ac3f7d728632793016d698f75f452053a2e0287b45c26414d0caf65e11f10 tests/smoke/mock_bin/STAR +25d54b2785e68da067602ce73418ab8e423ac0d37a390dbb0d4baf496f9d6e78 tests/smoke/mock_bin/cutadapt +3ef5ebfa90ed8e407eea5b712290b3d073b1bf074a4cc37cbe7885412c0b6d89 tests/smoke/mock_bin/python +ebe67cefd2ff77407eeef562beecf7b7e2886172e83dc5eb0f44b41a2940045a tests/smoke/mock_bin/samtools +01fec2b146a92eb676f72a8d8918bf6db22490e4493025553756f9eac78900bf tests/smoke/run_smoke.sh +df96c28ac7f64a27581718ae72acaaa9495f31b08d9e78034c1a5ce0e926002a tests/test_summary_report.py diff --git a/repro/manifest.tsv b/repro/manifest.tsv index 085c835..e35a511 100644 --- a/repro/manifest.tsv +++ b/repro/manifest.tsv @@ -1,21 +1,34 @@ artifact_id role path_or_accession external_url_or_doi sha256 size_bytes license_or_terms notes -source_tree code . https://github.com/TheDongLab/EVscope NA NA MIT Use exact public Git tag/commit for authoritative source-tree identity; local candidate tag is manuscript-v5-rigorous-local until public release. -README_md documentation README.md NA 364b98b069a67c435d7ec3d216e48bf57dfb8506f2a747a0dcccfdfb8e7f6188 20989 MIT User-facing installation, usage, scope and smoke-test documentation. +source_tree code . https://github.com/TheDongLab/EVscope NA NA MIT Use the exact public Git tag/commit for authoritative source-tree identity. +README_md documentation README.md NA 1395372d2fef6dc6ae2edc3fd0487924ca36f78ad848dd594dfe76632dc0ac8e 29089 MIT User-facing installation, kit selection, usage, scope and smoke-test documentation. LICENSE license LICENSE NA 8003003ded85bbc27a43f9719d604a358993818007a6dc7a35eb45e06c749e15 1082 MIT MIT license for EVscope source code. -EVscope_sh code EVscope.sh NA 0e3c0bb85dc8303c9cce36bece86ff406c3a61fc11f66e1f218efa921d65b0ed 100517 MIT Main Bash entry point. -EVscope_conf configuration EVscope.conf NA bf94ca778d2d1c0295cddcb8a451975d6b4fadbca784cabee1ee0be3b09d956b 6104 MIT Default reference/tool path configuration; users should edit local paths. +EVscope_sh code EVscope.sh NA d990fbfef5f7737e998ebc5f53600e1f3378f2186bda460dee2facb8c14207e6 111426 MIT Main Bash entry point with Total RNA and SMARTer smRNA kit profiles and separate Step 27 reports. +EVscope_conf configuration EVscope.conf NA 33ba6396cca8ff2301be81b394c64bbf45327b14153c74b6a4bb57787d37b6db 7713 MIT Default reference/tool path and injectable STAR-parameter configuration; users should edit local paths. bin_Step_25_EMapper_py code bin/Step_25_EMapper.py NA d13483469a4e47a25104fd7619289ce106d0fcfc25ebb055aed5dbfcd8a52ba2 117075 MIT EMapper EM-weighted BigWig/coverage and gene-summary module. bin_Step_25_bigWig2Expression_py code bin/Step_25_bigWig2Expression.py NA 412b4c1d64f511918031a44896ae7be0f0929ebd25c0261d76e84d540b77ae04 37500 MIT BigWig signal-per-million/TPM-like annotation summary module. bin_Step_26_density_plot_over_RNA_types_sh code bin/Step_26_density_plot_over_RNA_types.sh NA 2b8bb70804ccd527e229c9016b428abf1cca0e5835d4cdd25b0c837352a1ca31 7894 MIT RNA-type BigWig density plotting wrapper. bin_Step_26_density_plot_over_meta_gene_sh code bin/Step_26_density_plot_over_meta_gene.sh NA 1778752d48d9da3eb8e5f6a99d6c2ecc5b941e686119edebe8835050801c7291 8516 MIT Meta-gene BigWig density plotting wrapper. -bin_Step_27_html_report_Rmd code bin/Step_27_html_report.Rmd NA 2919740bc39bae68cefa0c1e7079c14407115a1089690e64d77b7fad39aa32eb 44066 MIT HTML reporting template. +bin_Step_07_bam2strand_py code bin/Step_07_bam2strand.py NA 270775851426458c286b59f31cbbeb74fa8e67633ab442d080e23813e7e57aed 30758 MIT Strandness reporting with kit-aware alignment basename handling. +bin_Step_24_generate_QC_matrix_py code bin/Step_24_generate_QC_matrix.py NA 8168717fc3ddd95b4eac9d566e7ef221ab4ffae682de682e9a0fae570b5a46dd 28878 MIT Kit-aware consolidated QC matrix generation. +bin_Step_27_html_report_Rmd code bin/Step_27_html_report.Rmd NA 064b4a9a91d2a6c5ba0f45dca2815ff97e24d19ea60fde16f63c54406d1029ec 45734 MIT Kit-aware comprehensive HTML reporting template with 90-percent data-table width. +bin_Step_27_generate_summary_report_py code bin/Step_27_generate_summary_report.py NA 117a2203bc4af999ae793a069275c7c073a9bda049360e891a408f082321d0ea 12125 MIT Standalone one-paragraph RNA-seq summary report generator. github_workflows_smoke_yml ci .github/workflows/smoke.yml NA 58eb25da24ec299883e2c1b3a60a972476fe3d0c7db4aef00c50966c6dfaa09e 521 MIT GitHub Actions smoke-test workflow. -tests_smoke_run_smoke_sh test tests/smoke/run_smoke.sh NA 2e10439cccafe7348a9c06529f79b33be14839e3c2b3a51774eb666a95d7f681 3436 MIT Local/CI smoke test covering syntax, CLI, toy FASTQ, tiny EMapper/BigWig when dependencies are available. +tests_smoke_run_smoke_sh test tests/smoke/run_smoke.sh NA 01fec2b146a92eb676f72a8d8918bf6db22490e4493025553756f9eac78900bf 8119 MIT Local/CI smoke test covering kit CLI validation, SMARTer smRNA paired trimming/no-dedup STAR alignment, summary reporting, config injection, resume isolation, and existing lightweight checks. +tests_test_summary_report_py test tests/test_summary_report.py NA df96c28ac7f64a27581718ae72acaaa9495f31b08d9e78034c1a5ce0e926002a 4686 MIT Synthetic unit tests for concise summary metrics, orientation, top-gene parsing, and HTML escaping. tests_smoke_data_toy_R2_fastq_gz test_data tests/smoke/data/toy_R2.fastq.gz NA 4f817aa0b325b273c2e94e57e19aee0e6c0a6e200e7f679687dd0a151c3cc7ce 72 MIT Tiny synthetic FASTQ for smoke testing. -tests_smoke_expected_expected_files_txt test tests/smoke/expected/expected_files.txt NA 0a4ad44f7e3244443e7647debc398b462333aca99a0266ccb28147b7b7a49754 214 MIT Expected smoke-test output paths. -repro_README_md documentation repro/README.md NA 4156fbc760db76834fb1ad62d3d9c57c7775f64c6ac2d7efdc9f92dc4006440d 941 MIT Release-level reproducibility index. -repro_manifest_tsv metadata repro/manifest.tsv NA NA 4904 MIT This manifest; checksum intentionally omitted to avoid self-reference. Validate with git commit/tag and repro/checksums.sha256. -repro_checksums_sha256 checksum repro/checksums.sha256 NA NA 4383 MIT Checksum file; self-checksum intentionally omitted. +tests_smoke_data_smarter_smrna_R1_fastq test_data tests/smoke/data/smarter_smrna_R1.fastq NA e8dbb54045c744a40b49f1ee069a82703b3fa9acf56f00c19b9c1a809e038412 294 MIT Synthetic SMARTer smRNA R1 records covering poly(A), minimum-length filtering and no-poly(A) reads. +tests_smoke_data_smarter_smrna_R2_fastq test_data tests/smoke/data/smarter_smrna_R2.fastq NA e835bc8327e9d41dc2fb0a1c741824336644eaaf87186e66f15f7baac748a7bd 344 MIT Synthetic SMARTer smRNA R2 records covering the 14-nt poly(T), reverse-complement insert and 3' adapter structure. +tests_smoke_expected_smarter_smrna_clean_fastq test tests/smoke/expected/smarter_smrna_clean.fastq NA 7bf9bf3bc51d54589ab699492aff2fcee116662a40840dcc9b1199127261ae30 122 MIT Expected clean FASTQ for the SMARTer smRNA fixture. +tests_smoke_expected_smarter_smrna_R2_clean_fastq test tests/smoke/expected/smarter_smrna_R2_clean.fastq NA 05dcb4361adcde7000e30fde4de63040222c6348b51072e224e993e010b5533d 122 MIT Expected clean R2 FASTQ for the paired SMARTer smRNA fixture. +tests_smoke_minimal_conf test tests/smoke/minimal.conf NA 9696908a77da99a50f3819a623880b508a45642ad3c3f01a3399d46aad40cb68 1261 MIT Minimal configuration with STAR profile and injection arrays for kit-aware pipeline smoke tests. +tests_smoke_mock_cutadapt test tests/smoke/mock_bin/cutadapt NA 25d54b2785e68da067602ce73418ab8e423ac0d37a390dbb0d4baf496f9d6e78 5933 MIT Deterministic paired-aware cutadapt test double for synthetic protocol trimming. +tests_smoke_mock_python test tests/smoke/mock_bin/python NA 3ef5ebfa90ed8e407eea5b712290b3d073b1bf074a4cc37cbe7885412c0b6d89 747 MIT Plotting test double used by the kit branch smoke test. +tests_smoke_mock_STAR test tests/smoke/mock_bin/STAR NA 3f2ac3f7d728632793016d698f75f452053a2e0287b45c26414d0caf65e11f10 808 MIT STAR test double used to verify no-dedup alignment routing. +tests_smoke_mock_samtools test tests/smoke/mock_bin/samtools NA ebe67cefd2ff77407eeef562beecf7b7e2886172e83dc5eb0f44b41a2940045a 640 MIT samtools test double used to verify no-dedup alignment routing. +tests_smoke_expected_expected_files_txt test tests/smoke/expected/expected_files.txt NA ea770e19e9b888df93ec92e7f0945baacc05ead12ab8124b431c5993cd8e04df 621 MIT Expected smoke-test output paths. +repro_README_md documentation repro/README.md NA 1e65d09e089ff5990606bde7312fe0960e7981aade3cd13446fa11e736017535 1196 MIT Release-level reproducibility index. +repro_manifest_tsv metadata repro/manifest.tsv NA NA 8068 MIT This manifest; checksum intentionally omitted to avoid self-reference. Validate with git commit/tag and repro/checksums.sha256. +repro_checksums_sha256 checksum repro/checksums.sha256 NA NA 5616 MIT Checksum file; self-checksum intentionally omitted. repro_metadata_software_versions_tsv metadata repro/metadata/software_versions.tsv NA e37c5a511bfa1d6235abe058cf2941f7ccfdc45fef7d9ff30a56c7387b588434 588 MIT Representative software/version metadata. repro_benchmarks_EMapper_MAQC_qPCR_README_md documentation repro/benchmarks/EMapper_MAQC_qPCR/README.md NA 767c5b88b1c8ca7349ba0efb4a71f758ef78bd27d62c4e36cae04b38ae675688 1064 MIT Benchmark/qPCR documentation and final-release artifact requirements. repro_figures_README_md documentation repro/figures/README.md NA 593f85dc44c8faeb6a5d5c2dcc036bf229e97f52139070259c193a8af753637b 1023 MIT Per-figure reproducibility organization notes. diff --git a/tests/smoke/README.md b/tests/smoke/README.md index abf9bad..6fe8b3d 100644 --- a/tests/smoke/README.md +++ b/tests/smoke/README.md @@ -1,9 +1,34 @@ # EVscope smoke tests -These tests validate repository installability and selected script-level behavior without requiring the full human reference bundle. They are intended for CI and release checks. Full end-to-end examples require the Zenodo reference bundle documented in the main README. +These tests validate repository installability, kit-aware command-line behavior, and selected script-level transformations without requiring the full human reference bundle. They are intended for CI and release checks. + +The v1.1 smoke coverage includes: + +- Bash and Python syntax plus CLI help/version checks. +- The canonical `--kit` values `smarter_total_rna` and `smarter_smrna`, including rejection of unsupported values. +- CLI help for the backward-compatible `smarter_total_rna` default. +- SMARTer smRNA paired-end validation, with R1 treated as sense and an implicit `forward` strand when `--strand` is not supplied. +- Kit-aware Steps 2, 3, and 6 behavior: the Total RNA UMI/ACC workflow remains available, while the smRNA profile applies the R1 `-u 3`/poly(A) trim, the R2 14-nt poly(T)/3' adapter trim, Q20 filtering of both mates, and paired no-UMI/no-deduplication alignment. +- STAR common/profile parameters and global/kit-specific extra parameters are injected from the configuration arrays without unsafe word splitting. +- SMARTer smRNA STAR profiles retain chimeric-junction output so circRNA Steps 8-10 can run without UMI deduplication. +- Step metadata records the normalized kit and effective strand, and cross-kit `--resume` is rejected. +- Selected toy FASTQ, read-length plotting, and EMapper/BigWig transformations when their optional dependencies are installed. +- Standalone Step 27 concise-summary generation, including human/contaminant mapping, genomic-region, strand-orientation, and top-gene reporting logic. + +The SMARTer smRNA pipeline checks use deterministic test doubles for cutadapt, STAR, and samtools, plus placeholder plots inside that mocked branch. They verify the synthetic trimming fixtures, command arguments, output routing, and metadata; they do not validate production alignment or circRNA accuracy. The separate optional plotting and EMapper checks run the real Python modules. + +The smoke suite does not run the complete 27-step workflow or render the comprehensive R Markdown report. Full end-to-end tests require the Zenodo reference bundle and R environment documented in the main README. All committed sequencing fixtures are synthetic; generated logs, reports, metadata, and temporary test environments stay under the ignored `work/` directory. Run: ```bash bash tests/smoke/run_smoke.sh ``` + +Temporary outputs are written under the ignored `tests/smoke/work/` directory. Optional checks report that they were skipped when their lightweight dependencies are unavailable. + +For a release candidate, also validate the tracked release files from the repository root: + +```bash +sha256sum -c repro/checksums.sha256 +``` diff --git a/tests/smoke/data/smarter_smrna_R1.fastq b/tests/smoke/data/smarter_smrna_R1.fastq new file mode 100644 index 0000000..3ab8ebb --- /dev/null +++ b/tests/smoke/data/smarter_smrna_R1.fastq @@ -0,0 +1,12 @@ +@smrna_22 +GGCACGTCCGTAGTCGATCGTACGCAAAAAAAAAAAAAAAAGATCGGAAGAGC ++ +IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII +@smrna_14 +TTATGCCTAGTCGATCGAAAAAAAAAAAAAAAAGATCGGAAGAGC ++ +IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII +@smrna_no_polya +CCCGGCTACGTCGATCGTACGTCAG ++ +IIIIIIIIIIIIIIIIIIIIIIIII diff --git a/tests/smoke/data/smarter_smrna_R2.fastq b/tests/smoke/data/smarter_smrna_R2.fastq new file mode 100644 index 0000000..99b6790 --- /dev/null +++ b/tests/smoke/data/smarter_smrna_R2.fastq @@ -0,0 +1,12 @@ +@smrna_22 +TTTTTTTTTTTTTTGCGTACGATCGACTACGGACGTGCCAGATCGGAAGAGC ++ +IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII +@smrna_14 +TTTTTTTTTTTTTTCGATCGACTAGGCATAAAGATCGGAAGAGC ++ +IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII +@smrna_no_polya +TTTTTTTTTTTTTTCTGACGTACGATCGACGTAGCCGGGAGATCGGAAGAGC ++ +IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII diff --git a/tests/smoke/expected/expected_files.txt b/tests/smoke/expected/expected_files.txt index 329e39c..ec6ebb6 100644 --- a/tests/smoke/expected/expected_files.txt +++ b/tests/smoke/expected/expected_files.txt @@ -1,4 +1,8 @@ work/acc.tsv +work/smrna/smarter_smoke_EVscope_output/Step_02_UMI_Analysis/smarter_smoke_library_profile.tsv +work/smrna/smarter_smoke_EVscope_output/Step_03_UMI_Adaptor_Trim/smarter_smoke_R1_clean.fq.gz +work/smrna/smarter_smoke_EVscope_output/Step_03_UMI_Adaptor_Trim/smarter_smoke_R2_clean.fq.gz +work/smrna/smarter_smoke_EVscope_output/Step_06_Alignment_Refined/smarter_smoke_STAR_no_dedup_Aligned.sortedByCoord.out.bam work/read_length.pdf work/read_length.png work/emapper/out/tiny_readcounts.txt diff --git a/tests/smoke/expected/smarter_smrna_R2_clean.fastq b/tests/smoke/expected/smarter_smrna_R2_clean.fastq new file mode 100644 index 0000000..b0fa5b8 --- /dev/null +++ b/tests/smoke/expected/smarter_smrna_R2_clean.fastq @@ -0,0 +1,8 @@ +@smrna_22 +GCGTACGATCGACTACGGACGT ++ +IIIIIIIIIIIIIIIIIIIIII +@smrna_no_polya +CTGACGTACGATCGACGTAGCC ++ +IIIIIIIIIIIIIIIIIIIIII diff --git a/tests/smoke/expected/smarter_smrna_clean.fastq b/tests/smoke/expected/smarter_smrna_clean.fastq new file mode 100644 index 0000000..684ef47 --- /dev/null +++ b/tests/smoke/expected/smarter_smrna_clean.fastq @@ -0,0 +1,8 @@ +@smrna_22 +ACGTCCGTAGTCGATCGTACGC ++ +IIIIIIIIIIIIIIIIIIIIII +@smrna_no_polya +GGCTACGTCGATCGTACGTCAG ++ +IIIIIIIIIIIIIIIIIIIIII diff --git a/tests/smoke/minimal.conf b/tests/smoke/minimal.conf new file mode 100644 index 0000000..35c384d --- /dev/null +++ b/tests/smoke/minimal.conf @@ -0,0 +1,30 @@ +# Minimal configuration used only by the kit-aware smoke tests. +EVscope_PATH="$(cd -P "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd -P)" +CORE_PYTHON_ENV="" +STAR_INDEX="${EVscope_PATH}/tests/smoke/data" +STAR_COMMON_ARGS=( + --readFilesCommand zcat + --outFilterMultimapNmax 100 + --winAnchorMultimapNmax 100 + --outSAMtype BAM SortedByCoordinate +) +STAR_TOTAL_RNA_ARGS=(--twopassMode Basic --chimSegmentMin 10) +STAR_SMRNA_SE_ARGS=( + --twopassMode None --alignEndsType EndToEnd --alignIntronMax 1 + --outFilterMismatchNmax 1 --outFilterMultimapScoreRange 0 + --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 + --outFilterScoreMin 12 --outFilterMatchNmin 14 + --chimSegmentMin 10 --chimJunctionOverhangMin 10 + --chimScoreMin 1 --chimOutType Junctions WithinBAM +) +STAR_SMRNA_PE_ARGS=( + --twopassMode None --alignEndsType EndToEnd --alignIntronMax 1 + --outFilterMismatchNmax 1 --outFilterMultimapScoreRange 0 + --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 + --outFilterScoreMin 24 --outFilterMatchNmin 28 + --chimSegmentMin 10 --chimJunctionOverhangMin 10 + --chimScoreMin 1 --chimOutType Junctions WithinBAM +) +STAR_EXTRA_ARGS=(--mockGlobalInjected yes) +STAR_TOTAL_RNA_EXTRA_ARGS=() +STAR_SMRNA_EXTRA_ARGS=(--mockSmrnaInjected yes) diff --git a/tests/smoke/mock_bin/STAR b/tests/smoke/mock_bin/STAR new file mode 100644 index 0000000..8790a82 --- /dev/null +++ b/tests/smoke/mock_bin/STAR @@ -0,0 +1,27 @@ +#!/usr/bin/env bash +# SPDX-License-Identifier: MIT +set -euo pipefail + +if [[ "${1:-}" == "--version" ]]; then + echo "STAR_2.7.11b (EVscope smoke mock)" + exit 0 +fi + +prefix="" +args=("$@") +for ((i = 0; i < ${#args[@]}; i++)); do + if [[ "${args[$i]}" == "--outFileNamePrefix" ]]; then + prefix="${args[$((i + 1))]}" + break + fi +done +[[ -n "$prefix" ]] +mkdir -p "$(dirname "$prefix")" +printf 'mock BAM\n' > "${prefix}Aligned.sortedByCoord.out.bam" +printf 'Number of input reads | 2\nUniquely mapped reads number | 2\nAverage input read length | 22\nAverage mapped length | 22\n' > "${prefix}Log.final.out" +printf 'mock junction\n' > "${prefix}Chimeric.out.junction" +printf 'mock splice junction\n' > "${prefix}SJ.out.tab" + +if [[ -n "${MOCK_STAR_LOG:-}" ]]; then + printf '%s\n' "$*" >> "$MOCK_STAR_LOG" +fi diff --git a/tests/smoke/mock_bin/cutadapt b/tests/smoke/mock_bin/cutadapt new file mode 100644 index 0000000..342647a --- /dev/null +++ b/tests/smoke/mock_bin/cutadapt @@ -0,0 +1,170 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: MIT +"""Small cutadapt test double for EVscope's SMARTer smRNA smoke test. + +It implements only the options used by the two Step 3 cutadapt calls. The +production pipeline continues to use the real cutadapt executable. +""" + +import gzip +import json +import os +import sys +from pathlib import Path + + +def open_text(path, mode): + if str(path).endswith(".gz"): + return gzip.open(path, mode, encoding="utf-8", newline="") + return open(path, mode, encoding="utf-8", newline="") + + +def parse_args(argv): + options = { + "cut": 0, + "cut2": 0, + "adapter": None, + "adapter2": None, + "minimum_length": 0, + "quality": None, + "quality2": None, + } + output = None + output2 = None + json_path = None + positional = [] + value_options = { + "-u": "cut", + "-U": "cut2", + "-a": "adapter", + "-A": "adapter2", + "-m": "minimum_length", + "--minimum-length": "minimum_length", + "-q": "quality", + "-Q": "quality2", + } + ignored_value_options = {"--overlap", "-j"} + i = 0 + while i < len(argv): + arg = argv[i] + if arg in value_options: + i += 1 + options[value_options[arg]] = argv[i] + elif arg in ignored_value_options: + i += 1 + elif arg == "-o": + i += 1 + output = argv[i] + elif arg == "-p": + i += 1 + output2 = argv[i] + elif arg == "--json": + i += 1 + json_path = argv[i] + elif arg.startswith("-"): + raise SystemExit(f"mock cutadapt: unsupported option {arg}") + else: + positional.append(arg) + i += 1 + + expected_inputs = 2 if output2 is not None else 1 + if output is None or len(positional) != expected_inputs: + raise SystemExit(f"mock cutadapt: expected {expected_inputs} input(s) and matching output(s)") + options["cut"] = int(options["cut"]) + options["cut2"] = int(options["cut2"]) + options["minimum_length"] = int(options["minimum_length"]) + options["quality"] = int(options["quality"]) if options["quality"] is not None else None + options["quality2"] = int(options["quality2"]) if options["quality2"] is not None else None + return positional, (output, output2), json_path, options + + +def read_fastq(path): + with open_text(path, "rt") as handle: + while True: + header = handle.readline() + if not header: + return + sequence = handle.readline().rstrip("\r\n") + plus = handle.readline() + quality = handle.readline().rstrip("\r\n") + if not plus or len(sequence) != len(quality): + raise SystemExit("mock cutadapt: invalid FASTQ input") + yield header.rstrip("\r\n"), sequence, plus.rstrip("\r\n"), quality + + +def trim_record(record, cut, adapter, quality_cutoff, minimum_length): + header, sequence, plus, quality = record + if cut > 0: + sequence = sequence[cut:] + quality = quality[cut:] + + if adapter: + adapter = adapter.split(";", 1)[0] + wildcard_prefix = len(adapter) - len(adapter.lstrip("N")) + search_adapter = adapter[wildcard_prefix:] + adapter_start = sequence.find(search_adapter) + if adapter_start >= wildcard_prefix: + adapter_start -= wildcard_prefix + if adapter_start >= 0: + sequence = sequence[:adapter_start] + quality = quality[:adapter_start] + + if quality_cutoff is not None: + while quality and ord(quality[-1]) - 33 < quality_cutoff: + sequence = sequence[:-1] + quality = quality[:-1] + + if len(sequence) < minimum_length: + return None + return header, sequence, plus, quality + + +def main(): + if "--version" in sys.argv[1:]: + print("cutadapt 4.9 (EVscope smoke mock)") + return + + input_paths, output_paths, json_path, options = parse_args(sys.argv[1:]) + invocation_log = os.environ.get("MOCK_CUTADAPT_LOG") + if invocation_log: + with open(invocation_log, "a", encoding="utf-8") as handle: + handle.write(" ".join(sys.argv[1:]) + "\n") + + written = 0 + output_path, output_path2 = output_paths + Path(output_path).parent.mkdir(parents=True, exist_ok=True) + if output_path2 is None: + with open_text(output_path, "wt") as output_handle: + for record in read_fastq(input_paths[0]): + trimmed = trim_record( + record, options["cut"], options["adapter"], + options["quality"], options["minimum_length"] + ) + if trimmed is None: + continue + output_handle.write("{}\n{}\n{}\n{}\n".format(*trimmed)) + written += 1 + else: + Path(output_path2).parent.mkdir(parents=True, exist_ok=True) + with open_text(output_path, "wt") as output_handle1, open_text(output_path2, "wt") as output_handle2: + for record1, record2 in zip(read_fastq(input_paths[0]), read_fastq(input_paths[1])): + trimmed1 = trim_record( + record1, options["cut"], options["adapter"], + options["quality"], options["minimum_length"] + ) + trimmed2 = trim_record( + record2, options["cut2"], options["adapter2"], + options["quality2"], options["minimum_length"] + ) + if trimmed1 is None or trimmed2 is None: + continue + output_handle1.write("{}\n{}\n{}\n{}\n".format(*trimmed1)) + output_handle2.write("{}\n{}\n{}\n{}\n".format(*trimmed2)) + written += 1 + + if json_path: + Path(json_path).write_text(json.dumps({"mock": True, "written": written}) + "\n", encoding="utf-8") + + +if __name__ == "__main__": + main() diff --git a/tests/smoke/mock_bin/python b/tests/smoke/mock_bin/python new file mode 100644 index 0000000..dafaf06 --- /dev/null +++ b/tests/smoke/mock_bin/python @@ -0,0 +1,32 @@ +#!/usr/bin/env bash +# SPDX-License-Identifier: MIT +set -euo pipefail + +# Step 3 plotting is not what this branch test exercises. CI also runs the real +# plotting smoke test below when matplotlib is installed. +if [[ "${1:-}" == *"Step_03_plot_fastq_read_length_dist.py" ]]; then + shift + output_pdf="" + output_png="" + while (( $# > 0 )); do + case "$1" in + --output_pdf) + output_pdf="$2" + shift 2 + ;; + --output_png) + output_png="$2" + shift 2 + ;; + *) + shift + ;; + esac + done + [[ -n "$output_pdf" && -n "$output_png" ]] + printf 'EVscope smoke plot placeholder\n' > "$output_pdf" + printf 'EVscope smoke plot placeholder\n' > "$output_png" + exit 0 +fi + +exec python3 "$@" diff --git a/tests/smoke/mock_bin/samtools b/tests/smoke/mock_bin/samtools new file mode 100644 index 0000000..4a418b5 --- /dev/null +++ b/tests/smoke/mock_bin/samtools @@ -0,0 +1,32 @@ +#!/usr/bin/env bash +# SPDX-License-Identifier: MIT +set -euo pipefail + +if [[ "${1:-}" == "--version" ]]; then + echo "samtools 1.21 (EVscope smoke mock)" + exit 0 +fi + +command_name="${1:-}" +shift || true +if [[ -n "${MOCK_SAMTOOLS_LOG:-}" ]]; then + printf '%s %s\n' "$command_name" "$*" >> "$MOCK_SAMTOOLS_LOG" +fi + +case "$command_name" in + index) + bam="${!#}" + printf 'mock BAI\n' > "${bam}.bai" + ;; + quickcheck) + bam="${!#}" + [[ -s "$bam" ]] + ;; + flagstat) + echo '2 + 0 in total (QC-passed reads + QC-failed reads)' + ;; + *) + echo "mock samtools: unsupported command ${command_name}" >&2 + exit 1 + ;; +esac diff --git a/tests/smoke/run_smoke.sh b/tests/smoke/run_smoke.sh index 799fa98..e0f752a 100755 --- a/tests/smoke/run_smoke.sh +++ b/tests/smoke/run_smoke.sh @@ -5,14 +5,25 @@ repo_dir=$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/../.." && pwd) cd "$repo_dir" mkdir -p tests/smoke/work -rm -rf tests/smoke/work/emapper +rm -rf tests/smoke/work/emapper tests/smoke/work/smrna +chmod +x tests/smoke/mock_bin/* -bash -n EVscope.sh -for sh in bin/*.sh; do bash -n "$sh"; done +tr -d '\r' < EVscope.sh | bash -n +for sh in bin/*.sh; do tr -d '\r' < "$sh" | bash -n; done python3 -m py_compile bin/*.py -bash EVscope.sh --help >/dev/null +python3 tests/test_summary_report.py +help_output=$(bash EVscope.sh --help) +grep -q -- '--kit ' <<< "$help_output" +grep -q 'smarter_total_rna' <<< "$help_output" +grep -q 'smarter_smrna' <<< "$help_output" bash EVscope.sh --version >/dev/null +if bash EVscope.sh --kit unsupported_kit >/dev/null 2>tests/smoke/work/invalid_kit.stderr; then + echo "unsupported --kit value was unexpectedly accepted" >&2 + exit 1 +fi +grep -q -- "--kit must be 'smarter_total_rna' or 'smarter_smrna'" tests/smoke/work/invalid_kit.stderr + python3 bin/Step_02_calculate_ACC_motif_fraction.py \ --input_fastq tests/smoke/data/toy_R2.fastq.gz \ --positions 1-3 \ @@ -20,6 +31,78 @@ python3 bin/Step_02_calculate_ACC_motif_fraction.py \ test -s tests/smoke/work/acc.tsv +smrna_work="$repo_dir/tests/smoke/work/smrna" +mkdir -p "$smrna_work" +( + cd "$smrna_work" + PATH="$repo_dir/tests/smoke/mock_bin:$PATH" \ + MOCK_CUTADAPT_LOG="$smrna_work/cutadapt.calls" \ + MOCK_STAR_LOG="$smrna_work/star.calls" \ + MOCK_SAMTOOLS_LOG="$smrna_work/samtools.calls" \ + bash "$repo_dir/EVscope.sh" \ + --sample_name smarter_smoke \ + --kit smarter_smrna \ + --input_fastqs "$repo_dir/tests/smoke/data/smarter_smrna_R1.fastq" "$repo_dir/tests/smoke/data/smarter_smrna_R2.fastq" \ + --run_steps 2-3,6 \ + --config "$repo_dir/tests/smoke/minimal.conf" +) + +smrna_output="$smrna_work/smarter_smoke_EVscope_output" +clean_fastq="$smrna_output/Step_03_UMI_Adaptor_Trim/smarter_smoke_R1_clean.fq.gz" +clean_fastq2="$smrna_output/Step_03_UMI_Adaptor_Trim/smarter_smoke_R2_clean.fq.gz" +profile_tsv="$smrna_output/Step_02_UMI_Analysis/smarter_smoke_library_profile.tsv" +step_meta="$smrna_output/Step_03_UMI_Adaptor_Trim/step.meta.json" +alignment_meta="$smrna_output/Step_06_Alignment_Refined/step.meta.json" +analysis_bam="$smrna_output/Step_06_Alignment_Refined/smarter_smoke_STAR_no_dedup_Aligned.sortedByCoord.out.bam" +gzip -dc "$clean_fastq" | diff -u tests/smoke/expected/smarter_smrna_clean.fastq - +gzip -dc "$clean_fastq2" | diff -u tests/smoke/expected/smarter_smrna_R2_clean.fastq - +awk -F '\t' '$1 == "UMI" && $2 == "No" {found=1} END {exit !found}' "$profile_tsv" +awk -F '\t' '$1 == "Input layout" && $2 == "Paired-end Read 1 / Read 2" {found=1} END {exit !found}' "$profile_tsv" +grep -q '"kit": "smarter_smrna"' "$step_meta" +grep -q '"strand": "forward"' "$step_meta" +grep -Fq -- '-u 3 -U 14 -a AAAAAAAAAA;min_overlap=3' "$smrna_work/cutadapt.calls" +grep -Fq -- '-A NNNAGATCGGAAGAGC;min_overlap=10' "$smrna_work/cutadapt.calls" +grep -Fq -- '-q 20 -Q 20 --minimum-length 15' "$smrna_work/cutadapt.calls" +test ! -e "$smrna_output/Step_03_UMI_Adaptor_Trim/smarter_smoke_R1_umi_tools.fq.gz" +test -s "$analysis_bam" +grep -q '"kit": "smarter_smrna"' "$alignment_meta" +grep -q '"strand": "forward"' "$alignment_meta" +test -s "$smrna_output/Step_06_Alignment_Initial/smarter_smoke_R1_no_dedup.clean.fq.gz" +test -s "$smrna_output/Step_06_Alignment_Initial/smarter_smoke_R2_no_dedup.clean.fq.gz" +test ! -e "$smrna_output/Step_06_Alignment_Initial/smarter_smoke_R1_umi_dedup.clean.fq.gz" +test "$(wc -l < "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--twopassMode None' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--alignEndsType EndToEnd' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--alignIntronMax 1' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--outFilterMismatchNmax 1' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--outFilterScoreMinOverLread 0' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--outFilterMatchNminOverLread 0' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--outFilterScoreMin 24' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--outFilterMatchNmin 28' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--mockGlobalInjected yes' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--mockSmrnaInjected yes' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--chimSegmentMin 10' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--chimJunctionOverhangMin 10' "$smrna_work/star.calls")" -eq 2 +test "$(grep -Fc -- '--chimOutType Junctions WithinBAM' "$smrna_work/star.calls")" -eq 2 +test -s "$smrna_output/Step_06_Alignment_Refined/smarter_smoke_STAR_no_dedup_Chimeric.out.junction" + +if ( + cd "$smrna_work" + PATH="$repo_dir/tests/smoke/mock_bin:$PATH" \ + bash "$repo_dir/EVscope.sh" \ + --sample_name smarter_smoke \ + --kit smarter_total_rna \ + --input_fastqs "$repo_dir/tests/smoke/data/smarter_smrna_R1.fastq" "$repo_dir/tests/smoke/data/smarter_smrna_R2.fastq" \ + --run_steps 2 \ + --resume \ + --config "$repo_dir/tests/smoke/minimal.conf" +) >/dev/null 2>"$smrna_work/cross_kit_resume.stderr" +then + echo "cross-kit --resume unexpectedly reused SMARTer smRNA results" >&2 + exit 1 +fi +grep -q 'metadata does not match current inputs/config/parameters' "$smrna_work/cross_kit_resume.stderr" + if python3 - <<'PY' import importlib.util, sys sys.exit(0 if importlib.util.find_spec('matplotlib') else 1) diff --git a/tests/test_summary_report.py b/tests/test_summary_report.py new file mode 100644 index 0000000..55ad6ec --- /dev/null +++ b/tests/test_summary_report.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +"""Tests for the standalone Step 27 concise summary report.""" + +import importlib.util +import tempfile +import unittest +from pathlib import Path + + +REPO_DIR = Path(__file__).resolve().parents[1] +SCRIPT_PATH = REPO_DIR / "bin" / "Step_27_generate_summary_report.py" +SPEC = importlib.util.spec_from_file_location("evscope_summary_report", SCRIPT_PATH) +SUMMARY = importlib.util.module_from_spec(SPEC) +assert SPEC and SPEC.loader +SPEC.loader.exec_module(SUMMARY) + + +class SummaryReportTests(unittest.TestCase): + def test_paired_total_rna_summary_contains_requested_metrics(self): + metrics = { + "Total Raw Read Pairs (Fragments)": "42650000", + "Total Trimmed Reads (R1+R2)": "78361898", + "Number of Reads after UMI-deduplication (STAR Input)": "1030000", + "Number of Uniquely Mapped Reads (STAR)": "440000", + "Number of Multi-mapped Reads (STAR)": "485000", + "Number of Reads Mapped to Escherichia coli (BBSplit)": "9081000", + "Percentage of Trimmed Reads Mapped to E. coli (BBSplit)": "11.59", + "Number of Reads Mapped to Mycoplasma (BBSplit)": "7766000", + "Percentage of Trimmed Reads Mapped to Mycoplasma (BBSplit)": "9.91", + "Number of Reads Mapped to Intronic Regions": "510000", + "Percentage of Reads Mapped to Intronic Regions": "55.22", + "Number of Reads Mapped to Exonic Regions": "58000", + "Percentage of Reads Mapped to Exonic Regions": "6.30", + "Number of Reads Mapped to Intergenic Regions": "278000", + "Percentage of Reads Mapped to Intergenic Regions": "30.05", + "Percentage of Mapped Reads on Forward Strand": "49.90", + "Percentage of Mapped Reads on Reverse Strand": "50.10", + } + genes = [(f"GENE{i}", float(1100 - i)) for i in range(1, 11)] + result = SUMMARY.build_summary( + "sample_A", metrics, genes, "smarter_total_rna", "unstrand", True + ) + + self.assertIn("42.65 million read pairs", result) + self.assertIn("39.18 million read pairs remaining after trimming", result) + self.assertIn("1.03 million complete read pairs", result) + self.assertIn("925,000 (89.81%) mapped to hg38", result) + self.assertIn("47.57% were unique and 52.43% were multi-mapped", result) + self.assertIn("11.59% of trimmed read pairs", result) + self.assertIn("9.91% of trimmed read pairs", result) + self.assertIn("55.22%", result) + self.assertIn("6.30%", result) + self.assertIn("30.05%", result) + self.assertIn("50.10% of strand-informative alignments", result) + self.assertIn("top 10 genes by read count", result) + self.assertIn("GENE10", result) + + def test_expression_parser_excludes_circrna_and_aggregates_symbols(self): + with tempfile.TemporaryDirectory() as tmpdir: + matrix = Path(tmpdir) / "expression.tsv" + matrix.write_text( + "GeneID\tGeneSymbol\tGeneType\tReadCounts\tNorm_Expr\n" + "ENSG1\tTTN\tprotein_coding\t100\t1\n" + "ENSG2\tTTN\tprotein_coding\t25\t1\n" + "circ1\tcirc-X\tcircRNAs\t1000\t1\n" + "ENSG3\tRERE\tprotein_coding\t80\t1\n", + encoding="utf-8", + ) + genes = SUMMARY.top_genes(str(matrix)) + + self.assertEqual(genes, [("TTN", 125.0), ("RERE", 80.0)]) + + def test_smrna_summary_reports_no_dedup_and_forward_orientation(self): + metrics = { + "Total Raw Reads (R1+R2)": "1000", + "Total Trimmed Reads (R1+R2)": "900", + "Number of Reads in Refined STAR Input (No UMI Deduplication)": "900", + "Number of Uniquely Mapped Reads (STAR)": "700", + "Number of Multi-mapped Reads (STAR)": "100", + "Percentage of Mapped Reads on Forward Strand": "72", + "Percentage of Mapped Reads on Reverse Strand": "18", + } + result = SUMMARY.build_summary( + "smrna_sample", metrics, [], "smarter_smrna", "forward", False + ) + + self.assertIn("No UMI deduplication was applied", result) + self.assertIn("80.00% of strand-informative alignments", result) + self.assertIn("expected forward/sense orientation", result) + + def test_html_contains_one_summary_paragraph(self): + report = SUMMARY.render_html("sample<&>", "A short & safe summary.") + self.assertEqual(report.count('

'), 1) + self.assertIn("sample<&>", report) + self.assertIn("A short & safe summary.", report) + + +if __name__ == "__main__": + unittest.main()