#!/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=99.slurm-featureCounts-%j.log # ============================================================================================================================================================= # == Load modules == # ============================================================================================================================================================= module load SAMtools/1.10-GCC-9.3.0 module load BEDTools/2.29.2-GCC-9.3.0 module load Biopython/1.78-foss-2020a-Python-3.8.2 module load Subread/2.0.2-GCC-9.3.0 # Global parameters needed for this script SOFTWAREDIR=/data/pg-molgen/software 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.bam featureCounts -s 2 -M -t $gff_type -g $gff_id -a $GFF3 -o $resultsdir/$name.2.counts $resultsdir/$name.bam # 2. Remove the comment created by featureCounts sed -i '1d' $resultsdir/$name.counts # 3. Add RPK and RPKM values python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/Counts_Convert.py -sessiondir $resultsdir -query $name.1.counts -out $name.1.rpkm python3 $SOFTWAREDIR/molgentools_mirror/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 "$genome.chromsizes Not found, creating one:" $SOFTWAREDIR/molgentools_mirror/rnaseq/make_chromsizes.pl -i $genome.fna -o $genome.chromsizes #$SOFTWAREDIR/molgentools_mirror/rnaseq/make_chromsizes.pl -i /data/pg-molgen/rnaseq/chenxi/MG1363.fna -o /data/pg-molgen/rnaseq/chenxi/MG1363.chromsizes 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 /data/pg-molgen/software/kentUtils/kentUtils/bin/linux.x86_64/bedGraphToBigWig $resultsdir/$name.plus.sort.bedgraph $genome.chromsizes $resultsdir/$name.plus.bw /data/pg-molgen/software/kentUtils/kentUtils/bin/linux.x86_64/bedGraphToBigWig $resultsdir/$name.min.sort.bedgraph $genome.chromsizes $resultsdir/$name.min.bw 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 /data/pg-molgen/software/kentUtils/kentUtils/bin/linux.x86_64/bedGraphToBigWig $resultsdir/$name.sort.bedgraph $genome.chromsizes $resultsdir/$name.bw fi } function my_chromsizes { # A simple routine to create $genome.chromsizes which is needed for in the coverage export PERL5LIB=$SOFTWAREDIR/molgentools_mirror/lib module load BioPerl/1.7.7-GCCcore-9.3.0 $SOFTWAREDIR/molgentools_mirror/rnaseq/make_chromsizes.pl -i $genome.fna -o $genome.chromsizes } function MergeTables { # 1. Merge tables: Combine counts or rpkm values from each file to one table python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/Merge_Tables.py -sessiondir $resultsdir -regex '.*1.counts$' -col 6 -o 00.merged.1.counts.txt python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/Merge_Tables.py -sessiondir $resultsdir -regex '.*2.counts$' -col 6 -o 00.merged.2.counts.txt python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/Merge_Tables.py -sessiondir $resultsdir -regex '.*1.rpkm$' -col 8 -o 00.merged.1.rpkm.txt python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/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 $SOFTWAREDIR/molgentools_mirror/rnaseq/Add_CV.py -sessiondir $resultsdir -query 00.merged.1.counts.txt -out 00.merged.1.counts.CV.txt python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/Add_CV.py -sessiondir $resultsdir -query 00.merged.2.counts.txt -out 00.merged.2.counts.CV.txt python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/Add_CV.py -sessiondir $resultsdir -query 00.merged.1.rpkm.txt -out 00.merged.1.rpkm.CV.txt python3 $SOFTWAREDIR/molgentools_mirror/rnaseq/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. before starting build the index files of the genome(s) # 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. Folder names; resultsdir=/data/p127804/PROJECTS/wilfried/20230212_008/results # 3. Feature counts # genome=/data/p127804/PROJECTS/wilfried/genomes/Bamyloliquefasciens_B1895_genome-plus-pBamB1895_startRapat75 # GFF3=$genome.gff # featureCount WM008_01 # featureCount WM008_02 # CoverageBed WM008_01 # CoverageBed WM008_02 genome=/data/p127804/PROJECTS/wilfried/genomes/Bpumilus_NRS576_genome-plus-p576_startRapAt75 GFF3=$genome.gff featureCount WM008_03 featureCount WM008_04 CoverageBed WM008_03 CoverageBed WM008_04 genome=/data/p127804/PROJECTS/wilfried/genomes/Bthuringiensis_kurstaki_HD73_chromosome-plus-pAW63_startRapAt75 GFF3=$genome.gff featureCount WM008_05 featureCount WM008_06 CoverageBed WM008_05 CoverageBed WM008_06 genome=/data/p127804/PROJECTS/wilfried/genomes/Bsp_NMCC4_genome-plus-pBspNMCC4_startRapAt75 GFF3=$genome.gff featureCount WM008_07 featureCount WM008_08 CoverageBed WM008_07 CoverageBed WM008_08 genome=/data/p127804/PROJECTS/wilfried/genomes/Blicheniformis_B4092_genomic-plus-pBli4092_reconstructed_startRapAt75 GFF3=$genome.gff featureCount WM008_09 featureCount WM008_10 CoverageBed WM008_09 CoverageBed WM008_10 genome=/data/p127804/PROJECTS/wilfried/genomes/Batrophaeus_NRS213_genomic-plus-pBatNRS213_startRapAt75 GFF3=$genome.gff featureCount WM008_11 featureCount WM008_12 CoverageBed WM008_11 CoverageBed WM008_12 genome=/data/pg-molgen/rnaseq/genomes/Bacillus_subtilis_subsp_subtilis_str_168_ASM904v1_pLS20 GFF3=$genome.gff featureCount WM008_13 featureCount WM008_14 CoverageBed WM008_13 CoverageBed WM008_14 genome=/data/pg-molgen/rnaseq/genomes/Bacillus_subtilis_subsp_subtilis_str_168_ASM904v1_pLS20 GFF3=$genome.gff featureCount WM008_15 featureCount WM008_16 CoverageBed WM008_15 CoverageBed WM008_16