#!/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.1.counts sed -i '1d' $resultsdir/$name.2.counts # 3. Add RPK and RPKM values python3 $AnneTools/Counts_Convert.py -sessiondir $resultsdir -query $name.1.counts -out $name.1.rpkm python3 $AnneTools/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 # ----------------- Bifidobacterium_adolescentis_MM0366 -------------------- workdir=/scratch/p127804/PROJECTS/Douwe/Ortensia_2026_031 bamdir=$workdir/results_MB0123 resultsdir=$workdir/featureCounts_MB0123 mkdir -p $resultsdir genome=/scratch/p127804/PROJECTS/Douwe/genomes/Bifidobacterium_adolescentis_MB0123 # featureCount OG031_01 # featureCount OG031_02 # featureCount OG031_03 # featureCount OG031_04 # featureCount OG031_05 # featureCount OG031_06 # featureCount OG031_07 # featureCount OG031_08 # featureCount OG031_09 # featureCount OG031_10 # featureCount OG031_11 # featureCount OG031_12 # featureCount OG031_13 # featureCount OG031_14 # featureCount OG031_15 # featureCount OG031_16 # MergeTables # # my_bamCoverage OG031_01 # my_bamCoverage OG031_02 # my_bamCoverage OG031_03 # my_bamCoverage OG031_04 # my_bamCoverage OG031_05 # my_bamCoverage OG031_06 # my_bamCoverage OG031_07 # my_bamCoverage OG031_08 # my_bamCoverage OG031_09 # my_bamCoverage OG031_10 # my_bamCoverage OG031_11 # my_bamCoverage OG031_12 # my_bamCoverage OG031_13 # my_bamCoverage OG031_14 # my_bamCoverage OG031_15 # my_bamCoverage OG031_16 # ----------------- Bifidobacterium_adolescentis_MM0366 -------------------- bamdir=$workdir/results_MM0366 resultsdir=$workdir/featureCounts_MM0366 mkdir -p $resultsdir bamdir=/scratch/p127804/PROJECTS/Douwe/Ortensia_2026_031/results_MM0366 genome=/scratch/p127804/PROJECTS/Douwe/genomes/Bifidobacterium_adolescentis_MM0366 featureCount OG031_17 featureCount OG031_18 featureCount OG031_19 featureCount OG031_20 featureCount OG031_21 featureCount OG031_22 featureCount OG031_23 featureCount OG031_24 MergeTables my_bamCoverage OG031_17 my_bamCoverage OG031_18 my_bamCoverage OG031_19 my_bamCoverage OG031_20 my_bamCoverage OG031_21 my_bamCoverage OG031_22 my_bamCoverage OG031_23 my_bamCoverage OG031_24