eval "$(/opt/anaconda/anaconda3/bin/conda shell.bash hook)"
conda activate
# Set number of threads
threads=42
# Set paths
ref_genome="../data/merged_out.fasta.TBtools.fa"
input_dir="../output/02-RNAseq"
output_dir="../output/04-align-v1"
# Loop over each FASTQ file
for fq in ${input_dir}/*.fastq.gz; do
# Extract sample name (e.g., Pum_barcode02 from Pum_barcode02_combined.fastq.gz)
sample=$(basename "$fq" .fastq.gz)
echo "Processing $sample..."
# Run minimap2 and samtools sort
minimap2 -ax splice -k14 --MD -t $threads "$ref_genome" "$fq" | \
samtools sort -@ $threads -o "${output_dir}/${sample}.sorted.bam"
# Index the sorted BAM
samtools index "${output_dir}/${sample}.sorted.bam"
done## Processing Pum_barcode02_combined...
## [M::mm_idx_gen::47.587*1.77] collected minimizers
## [M::mm_idx_gen::52.725*2.66] sorted minimizers
## [M::main::52.725*2.66] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::53.392*2.64] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::53.798*2.63] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::154.327*23.10] mapped 558459 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Pum_barcode02_combined.fastq.gz
## [M::main] Real time: 154.512 sec; CPU: 3565.014 sec; Peak RSS: 22.469 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Pum_barcode12_1_combined...
## [M::mm_idx_gen::34.380*1.74] collected minimizers
## [M::mm_idx_gen::36.267*3.08] sorted minimizers
## [M::main::36.267*3.08] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::36.902*3.05] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::37.300*3.02] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::153.440*29.59] mapped 914269 sequences
## [M::worker_pipeline::154.685*29.36] mapped 111278 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Pum_barcode12_1_combined.fastq.gz
## [M::main] Real time: 154.953 sec; CPU: 4542.344 sec; Peak RSS: 22.453 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Pum_barcode13_1_combined...
## [M::mm_idx_gen::24.726*1.76] collected minimizers
## [M::mm_idx_gen::43.691*2.49] sorted minimizers
## [M::main::43.691*2.49] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::44.867*2.45] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::45.616*2.42] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::171.556*20.40] mapped 928290 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Pum_barcode13_1_combined.fastq.gz
## [M::main] Real time: 171.744 sec; CPU: 3499.192 sec; Peak RSS: 20.439 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Pum_barcode13_combined...
## [M::mm_idx_gen::41.288*1.80] collected minimizers
## [M::mm_idx_gen::61.005*2.34] sorted minimizers
## [M::main::61.005*2.34] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::62.247*2.32] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::63.018*2.30] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::203.703*11.53] mapped 572128 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Pum_barcode13_combined.fastq.gz
## [M::main] Real time: 204.021 sec; CPU: 2348.543 sec; Peak RSS: 19.505 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Qui_barcode01_1_combined...
## [M::mm_idx_gen::46.678*1.76] collected minimizers
## [M::mm_idx_gen::56.262*2.55] sorted minimizers
## [M::main::56.262*2.55] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::56.933*2.53] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::57.345*2.52] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::141.200*22.06] mapped 933072 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Qui_barcode01_1_combined.fastq.gz
## [M::main] Real time: 141.340 sec; CPU: 3115.589 sec; Peak RSS: 19.964 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Qui_barcode01_combined...
## [M::mm_idx_gen::21.172*1.77] collected minimizers
## [M::mm_idx_gen::23.078*3.87] sorted minimizers
## [M::main::23.078*3.87] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::23.717*3.79] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::24.113*3.75] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::128.579*21.39] mapped 573113 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Qui_barcode01_combined.fastq.gz
## [M::main] Real time: 128.712 sec; CPU: 2750.190 sec; Peak RSS: 22.330 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Qui_barcode02_1_combined...
## [M::mm_idx_gen::21.349*1.76] collected minimizers
## [M::mm_idx_gen::23.273*3.86] sorted minimizers
## [M::main::23.273*3.86] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::23.916*3.78] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::24.315*3.74] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::642.756*8.82] mapped 837192 sequences
## [M::worker_pipeline::653.263*8.86] mapped 219281 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Qui_barcode02_1_combined.fastq.gz
## [M::main] Real time: 653.993 sec; CPU: 5787.312 sec; Peak RSS: 27.366 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Qui_barcode12_combined...
## [M::mm_idx_gen::21.996*1.78] collected minimizers
## [M::mm_idx_gen::23.914*3.82] sorted minimizers
## [M::main::23.914*3.82] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::24.602*3.74] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::25.005*3.70] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::116.317*28.71] mapped 794489 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Qui_barcode12_combined.fastq.gz
## [M::main] Real time: 116.626 sec; CPU: 3340.009 sec; Peak RSS: 22.345 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Rio_barcode03_1_combined...
## [M::mm_idx_gen::46.681*1.77] collected minimizers
## [M::mm_idx_gen::65.912*2.29] sorted minimizers
## [M::main::65.912*2.29] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::66.965*2.27] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::67.614*2.26] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::202.209*23.45] mapped 917008 sequences
## [M::worker_pipeline::203.632*23.30] mapped 135071 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Rio_barcode03_1_combined.fastq.gz
## [M::main] Real time: 203.940 sec; CPU: 4744.169 sec; Peak RSS: 21.564 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Rio_barcode03_combined...
## [M::mm_idx_gen::46.728*1.76] collected minimizers
## [M::mm_idx_gen::55.564*2.54] sorted minimizers
## [M::main::55.564*2.54] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::56.259*2.53] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::56.665*2.51] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::173.727*19.02] mapped 631809 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Rio_barcode03_combined.fastq.gz
## [M::main] Real time: 173.913 sec; CPU: 3304.477 sec; Peak RSS: 22.401 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Rio_barcode14_1_combined...
## [M::mm_idx_gen::23.844*1.78] collected minimizers
## [M::mm_idx_gen::25.946*3.79] sorted minimizers
## [M::main::25.946*3.79] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::26.602*3.72] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::27.009*3.68] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::142.849*31.10] mapped 825271 sequences
## [M::worker_pipeline::143.475*30.97] mapped 47391 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Rio_barcode14_1_combined.fastq.gz
## [M::main] Real time: 143.722 sec; CPU: 4443.001 sec; Peak RSS: 21.693 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...
## Processing Rio_barcode14_combined...
## [M::mm_idx_gen::20.981*1.78] collected minimizers
## [M::mm_idx_gen::22.861*3.88] sorted minimizers
## [M::main::22.861*3.88] loaded/built the index for 10 target sequence(s)
## [M::mm_mapopt_update::23.489*3.80] mid_occ = 1506
## [M::mm_idx_stat] kmer size: 14; skip: 5; is_hpc: 0; #seq: 10
## [M::mm_idx_stat::23.878*3.76] distinct minimizers: 43379551 (28.91% are singletons); average occurrences: 9.236; average spacing: 2.947
## [M::worker_pipeline::106.025*28.91] mapped 603887 sequences
## [M::main] Version: 2.17-r941
## [M::main] CMD: minimap2 -ax splice -k14 --MD -t 42 ../data/merged_out.fasta.TBtools.fa ../output/02-RNAseq/Rio_barcode14_combined.fastq.gz
## [M::main] Real time: 106.169 sec; CPU: 3065.235 sec; Peak RSS: 21.946 GB
## [bam_sort_core] merging from 0 files and 42 in-memory blocks...find . -name "*.bam" | parallel 'bedtools genomecov -ibam {} -bg > {.}.bedGraph'samtools faidx Apulchra-genome.fa cut -f1,2 Apulchra-genome.fa.fai > genome.chrom.sizes
for bed in *.bedGraph; do base=\((basename "\)bed” .bedGraph) bedGraphToBigWig “\(bed" genome.chrom.sizes "\){base}.bw” done