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..4b356a4 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** 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 @@ -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. + +**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..b88bde6 100755 --- a/bin/Step_27_html_report.Rmd +++ b/bin/Step_27_html_report.Rmd @@ -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..e2379c2 100644 --- a/repro/README.md +++ b/repro/README.md @@ -13,4 +13,4 @@ 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. -`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..599fd15 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 +760041c46f7042dc1eb0cc0c8e22a20ef3ac314ddfbf9fc85bee87f8c736553a 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,20 +23,21 @@ 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 +772cf2bc3ab3b4b00b87452cf9061d819613e94ed410a2308255b93c1fffdc02 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 +c9fe5113e00d09485eca2b3ecc9d9f021a7ca927e92e6d4d492c70d4b047bbef repro/README.md 767c5b88b1c8ca7349ba0efb4a71f758ef78bd27d62c4e36cae04b38ae675688 repro/benchmarks/EMapper_MAQC_qPCR/README.md 593f85dc44c8faeb6a5d5c2dcc036bf229e97f52139070259c193a8af753637b repro/figures/README.md -62fa01923420b924e58ac4ac453da12c8b4f03f6cd655f6248d764df026f3d56 repro/manifest.tsv +34e1d1cd6150053c85afcc33b847dfe8e67e1da6b1f0292566230b1e6c3d242d repro/manifest.tsv e37c5a511bfa1d6235abe058cf2941f7ccfdc45fef7d9ff30a56c7387b588434 repro/metadata/software_versions.tsv 61f65dd9598c7421870bf4918332637f83f34abe1f9fc7cc095852ef7a1b44f7 tests/smoke/.gitignore bf94bd4d4c60a43890acafd0c456ca4d692f009003ba7fee3064887fb39da7ed tests/smoke/README.md diff --git a/repro/manifest.tsv b/repro/manifest.tsv index 085c835..a01b1ce 100644 --- a/repro/manifest.tsv +++ b/repro/manifest.tsv @@ -1,21 +1,24 @@ 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 760041c46f7042dc1eb0cc0c8e22a20ef3ac314ddfbf9fc85bee87f8c736553a 28800 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 772cf2bc3ab3b4b00b87452cf9061d819613e94ed410a2308255b93c1fffdc02 45743 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_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. +repro_README_md documentation repro/README.md NA c9fe5113e00d09485eca2b3ecc9d9f021a7ca927e92e6d4d492c70d4b047bbef 1031 MIT Release-level reproducibility index. +repro_manifest_tsv metadata repro/manifest.tsv NA NA 5628 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 4569 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.