#!/bin/bash # RNA-Seq # 2022-02 # ============================================================================================================================================================= # == 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 24 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=2 #SBATCH --mem=8GB # Choose between long or short queue: A short job will start quickly after submitting, but will be terminated after 30 minutes # Normally short is long enough #SBATCH --time=24:00:00 ##SBATCH -p short # Personal stuff; Change this to your personal settings #SBATCH --job-name=featureCounts ##SBATCH --mail-type=END ##SBATCH --mail-user=anne.de.jong@rug.nl #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 #export PERL5LIB=$SOFTWAREDIR/molgentools_mirror/lib # ============================================================================================================================================================= # == Functions == # ============================================================================================================================================================= function featureCount { 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 $resultsdir/$name.sorted.bam featureCounts -s 2 -M -t $gff_type -g $gff_id -a $GFF3 -o $resultsdir/$name.2.counts $resultsdir/$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 $resultsdir/$name.sorted.bam -g $genome.chromsizes -strand + > $resultsdir/$name.plus.coverage #genomeCoverageBed -d -ibam $resultsdir/$name.sorted.bam -g $genome.chromsizes -strand - > $resultsdir/$name.minus.coverage genomeCoverageBed -bg -ibam $resultsdir/$name.sorted.bam -g $genome.chromsizes -strand + > $resultsdir/$name.plus.bedgraph genomeCoverageBed -bg -ibam $resultsdir/$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 $resultsdir/$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 $resultsdir/$name.sorted.bam -o $resultsdir/$name.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 genome=/scratch/p127804/PROJECTS/Douwe/genomes/Co-culture_NRBB57_CNCM_I-4319/co-culture GFF3=$genome.gff3 resultsdir=/scratch/p127804/PROJECTS/Douwe/Zhe/results_Zhe #featureCount ZH018_01 #featureCount ZH018_02 #featureCount ZH018_03 #featureCount ZH018_04 #featureCount ZH018_05 #featureCount ZH018_06 #featureCount ZH018_07 #featureCount ZH018_08 # #MergeTables my_bamCoverage ZH018_01 my_bamCoverage ZH018_02 my_bamCoverage ZH018_03 my_bamCoverage ZH018_04 my_bamCoverage ZH018_05 my_bamCoverage ZH018_06 my_bamCoverage ZH018_07 my_bamCoverage ZH018_08