#!/bin/bash # RNA-Seq # 2025-04 # ============================================================================================================================================================= # == SLURM parameters == # ============================================================================================================================================================= # Decide how much computer power you need # The programs used in this script can only be executed on 1 NODE # You can use a max of 128 task per node. But the lower the number of tasks requested the shorter the queue when submitting a job #SBATCH --nodes=1 #SBATCH --ntasks-per-node=32 #SBATCH --mem=8GB #SBATCH --time=24:00:00 #SBATCH --job-name=featureCounts #SBATCH --output=slurm-featureCounts-%j.log # ============================================================================================================================================================= # == Load modules == # ============================================================================================================================================================= module load SAMtools/1.16.1-GCC-11.3.0 module load BEDTools/2.30.0-GCC-11.3.0 module load Biopython/1.81-foss-2022b module load Subread/2.0.4-GCC-11.3.0 module load deepTools/3.5.1-foss-2021b # Global parameters needed for this script # setting global parameters cpu=$SLURM_JOB_CPUS_PER_NODE SOFTWAREDIR=/home1/p127804/software AnneTools=/userapps/hb-molgen/AnneTools # ============================================================================================================================================================= # == Functions == # ============================================================================================================================================================= function featureCount { GFF3=$genome.gff name=$1 # Once the BAM file is ready, this can be combined with genome annotation to calculate how many reads are mapped per feature # Features are given in the .GFF or .GFF3 file format # Note that sequence the RNA library prep can be done stranded or unstranded. Be sure to select the proper option # -s = 0 (unstranded), 1 (stranded) and 2 (reversely stranded) printf "\n========== featureCount: "$name" ======================================\n" # 1. Add number of reads per feature featureCounts -s 1 -M -t $gff_type -g $gff_id -a $GFF3 -o $resultsdir/$name.1.counts $bamdir/$name.sorted.bam featureCounts -s 2 -M -t $gff_type -g $gff_id -a $GFF3 -o $resultsdir/$name.2.counts $bamdir/$name.sorted.bam # 2. Remove the comment created by featureCounts sed -i '1d' $resultsdir/$name.counts # 3. Add RPK and RPKM values python3 $AnneTools/rnaseq/Counts_Convert.py -sessiondir $resultsdir -query $name.1.counts -out $name.1.rpkm python3 $AnneTools/rnaseq/Counts_Convert.py -sessiondir $resultsdir -query $name.2.counts -out $name.2.rpkm } function CoverageBed { name=$1 # This is an optional function to make BigWig and BED files that can used in various other programs # Depends in chromsizes file which can be generated using /data/molgentools/rnaseq/make_chromsizes.pl printf "\n========== CoverageBed: "$name" ======================================\n" stranded=true if [ ! -e $genome.chromsizes ]; then echo "ERROR, $genome.chromsizes Not found!" elif [ $stranded ]; then printf "genomeCoverageBed for $name\n" #genomeCoverageBed -d -ibam $bamdir/$name.sorted.bam -g $genome.chromsizes -strand + > $resultsdir/$name.plus.coverage #genomeCoverageBed -d -ibam $bamdir/$name.sorted.bam -g $genome.chromsizes -strand - > $resultsdir/$name.minus.coverage genomeCoverageBed -bg -ibam $bamdir/$name.sorted.bam -g $genome.chromsizes -strand + > $resultsdir/$name.plus.bedgraph genomeCoverageBed -bg -ibam $bamdir/$name.sorted.bam -g $genome.chromsizes -strand - > $resultsdir/$name.minus.bedgraph sort -k1,1 -k2,2n $resultsdir/$name.plus.bedgraph > $resultsdir/$name.plus.sort.bedgraph sort -k1,1 -k2,2n $resultsdir/$name.minus.bedgraph > $resultsdir/$name.min.sort.bedgraph else genomeCoverageBed -bg -ibam $bamdir/$name.sorted.bam -g $genome.chromsizes > $resultsdir/$name.bedgraph sort -k1,1 -k2,2n $resultsdir/$name.bedgraph > $resultsdir/$name.sort.bedgraph fi } function my_bamCoverage { name=$1 printf "\n========== DeepTools: "$name" ======================================\n" #bamCoverage -b $bamdir/$name.sorted.bam -o $resultsdir/$name.bw --binSize 10 --normalizeUsing RPKM --smoothLength 30 bamCoverage -b $bamdir/$name.sorted.bam --filterRNAstrand=forward -o $resultsdir/$name.forward.bw --binSize 10 --normalizeUsing RPKM --smoothLength 30 bamCoverage -b $bamdir/$name.sorted.bam --filterRNAstrand=reverse -o $resultsdir/$name.reverse.bw --binSize 10 --normalizeUsing RPKM --smoothLength 30 } function MergeTables { # 1. Merge tables: Combine counts or rpkm values from each file to one table python3 $AnneTools/Merge_Tables.py -sessiondir $resultsdir -regex '.*1.counts$' -col 6 -o 00.merged.1.counts.txt python3 $AnneTools/Merge_Tables.py -sessiondir $resultsdir -regex '.*2.counts$' -col 6 -o 00.merged.2.counts.txt python3 $AnneTools/Merge_Tables.py -sessiondir $resultsdir -regex '.*1.rpkm$' -col 8 -o 00.merged.1.rpkm.txt python3 $AnneTools/Merge_Tables.py -sessiondir $resultsdir -regex '.*2.rpkm$' -col 8 -o 00.merged.2.rpkm.txt # 2. Add Column with Coefficient of Variance (CV) per feature python3 $AnneTools/Add_CV.py -sessiondir $resultsdir -query 00.merged.1.counts.txt -out 00.merged.1.counts.CV.txt python3 $AnneTools/Add_CV.py -sessiondir $resultsdir -query 00.merged.2.counts.txt -out 00.merged.2.counts.CV.txt python3 $AnneTools/Add_CV.py -sessiondir $resultsdir -query 00.merged.1.rpkm.txt -out 00.merged.1.rpkm.CV.txt python3 $AnneTools/Add_CV.py -sessiondir $resultsdir -query 00.merged.2.rpkm.txt -out 00.merged.2.rpkm.CV.txt # 3. Clean the table from low values, if needed #$SOFTWAREDIR/molgentools_mirror/tables/remove_sum_zero.pl 00.merged.counts.txt 00.merged_clean_sum100.table } function my_mkdir { if [ ! -d $1 ]; then mkdir $1 ; fi } # ============================================================================================================================================================= # == Main script == # ============================================================================================================================================================= # Most logical struture is to make 3 folders in the SESSIONDIR # /fastq ==> Copy here all fastq reads # /genomes ==> All genome info # /scripts ==> Location of this script # a /results folder will be created by this script # 1.# Study the GFF file and decide what you want to use for feature counts # gff_type is the 3th column of the GFF file. Optionally, change this for your needs # gff_id is one of the fields from the 9th column of the GFF file. Optionally, change this for your needs gff_type='gene' gff_id='locus_tag' # 2. featureCount bamdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/results resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/ mkdir -p $resultsdir genome=/scratch/p127804/PROJECTS/Douwe/genomes/Bifidobacterium_breve_NRBB57 resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_breve_NRBB57 mkdir -p $resultsdir featureCount DS027_29 featureCount DS027_30 MergeTables my_bamCoverage DS027_29 my_bamCoverage DS027_30 genome=/scratch/p127804/PROJECTS/Douwe/genomes/Bifidobacterium_breve_NRBB01 resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_breve_NRBB01 mkdir -p $resultsdir featureCount DS027_31 featureCount DS027_32 MergeTables my_bamCoverage DS027_31 my_bamCoverage DS027_32 genome=/scratch/p127804/PROJECTS/Douwe/genomes/Bifidobacterium_bifidum_CNCM_I_4319 resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_bifidum_CNCM_I_4319 mkdir -p $resultsdir featureCount DS027_33 featureCount DS027_34 MergeTables my_bamCoverage DS027_33 my_bamCoverage DS027_34 genome=/scratch/p127804/PROJECTS/Douwe/genomes/Bifidobacterium_longum_subsp_infantis_Bll190 resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_longum_subsp_infantis_Bll190 mkdir -p $resultsdir featureCount DS027_35 featureCount DS027_36 featureCount DS027_37 featureCount DS027_38 MergeTables my_bamCoverage DS027_35 my_bamCoverage DS027_36 my_bamCoverage DS027_37 my_bamCoverage DS027_38 genome=/scratch/p127804/PROJECTS/Douwe/genomes/2025_Bifido_mix/Bifidobacterium_mixA resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_mixA mkdir -p $resultsdir featureCount DS027_39 featureCount DS027_40 MergeTables my_bamCoverage DS027_39 my_bamCoverage DS027_40 genome=/scratch/p127804/PROJECTS/Douwe/genomes/2025_Bifido_mix/Bifidobacterium_mixB resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_mixB mkdir -p $resultsdir featureCount DS027_41 featureCount DS027_42 featureCount DS027_43 featureCount DS027_44 MergeTables my_bamCoverage DS027_41 my_bamCoverage DS027_42 my_bamCoverage DS027_43 my_bamCoverage DS027_44 genome=/scratch/p127804/PROJECTS/Douwe/genomes/2025_Bifido_mix/Bifidobacterium_mixC resultsdir=/scratch/p127804/PROJECTS/Douwe/Douwe_202512_028/Bifidobacterium_mixC mkdir -p $resultsdir featureCount DS027_25 featureCount DS027_26 featureCount DS027_27 featureCount DS027_28 featureCount DS027_45 featureCount DS027_46 featureCount DS027_47 featureCount DS027_48 MergeTables my_bamCoverage DS027_25 my_bamCoverage DS027_26 my_bamCoverage DS027_27 my_bamCoverage DS027_28 my_bamCoverage DS027_45 my_bamCoverage DS027_46 my_bamCoverage DS027_47 my_bamCoverage DS027_48