#!/bin/bash # RNA-Seq Single End RNA-Seq # 2020-05-13 #SBATCH --nodes=1 #SBATCH --ntasks-per-node=16 #SBATCH --mem=64GB #SBATCH --time=24:00:00 ##SBATCH -p short #SBATCH --job-name=RNA_Seq_SE # parameters cpu=16 SCRATCHDIR=/scratch/p276244 SOFTWAREDIR=/data/pg-molgen/software # ============================================================================================================================================================= # == Load modules == # ============================================================================================================================================================= module load SAMtools/1.9-GCC-6.4.0-2.28 module load BEDTools/2.25.0-foss-2016a module load FastQC/0.11.7-Java-1.8.0_162 module load Trimmomatic/0.36-Java-1.8.0_92 export PERL5LIB=$SOFTWAREDIR/molgentools_mirror/lib module load BioPerl/1.6.924-intel-2016a-Perl-5.20.3 module load Bowtie2/2.3.4.3-foss-2018a # option for the IonProton bowtie_options=' -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 ' # ============================================================================================================================================================= # == Functions == # ============================================================================================================================================================= function trimmomatic { printf "\n========== Trimmomatic: "$name" ======================================\n" java -jar $EBROOTTRIMMOMATIC/trimmomatic-0.36.jar SE -phred33 -threads $cpu $datadir/$forward $datadir/$forward_trimmed ILLUMINACLIP:$SOFTWAREDIR/trimmomatic/TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 } function my_fastqc { # The forward reads are used for quality control printf "\n========== FastQC: "$name" ======================================\n" my_mkdir $resultdir my_mkdir $resultdir/fastqc my_mkdir $resultdir/fastqc/$name.trimmed fastqc -t $cpu -o $resultdir/fastqc/$name.trimmed -f fastq $datadir/$forward_trimmed } function my_kraken { # kraken is used to identify the taxonomic position of the sample printf "\n========== Kraken: "$name" ======================================\n" krakendir=$resultdir/kraken my_mkdir $krakendir $SOFTWAREDIR/kraken/kraken2/kraken2 --threads $cpu --unclassified -out $krakendir/$name.unclassified.fasta --classified-out $krakendir/$name.classified.fasta $resultdir/$name.400bp.fasta --report $krakendir/$name.report > $krakendir/$name.kraken.table } function my_diamond { # Direct mapping of fastq reads to the protein database printf "\n========== DIAMOND: "$name" ======================================\n" # Example: https://metagenomics-workshop.readthedocs.io/en/latest/annotation/taxonomic_annotation.html #DIAMOND_DB=/data/pg-molgen/databases/UniRef90/UniRef90.dmnd DIAMOND_DB=/data/pg-molgen/databases/RefSeq/RefSeq.dmnd diamondresultdir=$resultdir/DIAMOND_RNA my_mkdir $diamondresultdir # Merge forward and reverse reads cat $datadir/$forward_trimmed $datadir/$reverse_trimmed > $datadir/$name.merged.fastq.gz $SOFTWAREDIR/diamond/diamond blastx --threads $cpu --tmpdir $SCRATCHDIR --query $datadir/$name.merged.fastq.gz --db $DIAMOND_DB --daa $diamondresultdir/$name.diamond.daa --un $diamondresultdir/$name.unaligned.fastq rm $datadir/$name.merged.fastq.gz # The results are stored in a binary format so to see it in plain text you need to convert with diamond view: $SOFTWAREDIR/diamond/diamond view -a $diamondresultdir/$name.diamond.daa > $diamondresultdir/$name.diamond.tab # use /data/pg-molgen/databases/id.idmapping.dat to map functions to UniProt } function my_metaspades { printf "\n========== metaSpades: "$name" ======================================\n" # --meta this flag is required for metagenomic sample data # --rna this flag is required for RNA-Seq data spades.py --rna -t $cpu -k 77,99 --pe1-1 $forward_trimmed --pe1-2 $reverse_trimmed -o SPAdes.$name.result $SOFTWAREDIR/molgentools_mirror/tools/rename_contig_NODES.pl -minlen 400 -mincov 1 -i ./SPAdes.$name.result/contigs.fasta -prefix $name.400bp } function my_bowtie2 { # map the RNA-Seq reads on all the MOCAT2 assembly of the DNA-Seq data # Gene calling was done using MetaGeneMark within MOCAT2 printf "\n\n========= Bowtie2: "$name" =================================\n" > $resultdir/02.$name.bowtie2.log (bowtie2 $bowtie_options --local -U $datadir/$forward_trimmed -p $cpu -S $resultdir/$name.sam -x $genome) 2>>$resultdir/02.$name.bowtie2.log # bowtie2 $bowtie_options --local -U $datadir/$forward_trimmed -p $cpu -S $resultdir/$name.sam -x $genome } function process_bam { printf "\n========== DIAMOND: "$name" ======================================\n" samtools view -S -b -o $resultdir/$name.bam $resultdir/$name.sam # == new in one go: samtools view -bS $resultdir/$name.sam | samtools sort - $resultdir/$name.bam rm $resultdir/$name.sam # save space by removing the SAM file samtools sort $resultdir/$name.bam > $resultdir/$name.sorted.bam samtools index $resultdir/$name.sorted.bam samtools depth $resultdir/$name.sorted.bam > $resultdir/$name.depth } function CoverageBed { # Depends in chromsizes file which can be generated using /data/molgentools/rnaseq/make_chromsizes.pl printf "\n========== DIAMOND: "$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 $resultdir/$name.sorted.bam -g $genome.chromsizes -strand + > $resultdir/$name.plus.coverage genomeCoverageBed -d -ibam $resultdir/$name.sorted.bam -g $genome.chromsizes -strand - > $resultdir/$name.minus.coverage genomeCoverageBed -bg -ibam $resultdir/$name.sorted.bam -g $genome.chromsizes -strand + > $resultdir/$name.plus.bedgraph genomeCoverageBed -bg -ibam $resultdir/$name.sorted.bam -g $genome.chromsizes -strand - > $resultdir/$name.minus.bedgraph # test sort -k1,1 -k2,2n $resultdir/$name.plus.bedgraph > $resultdir/$name.plus.sort.bedgraph sort -k1,1 -k2,2n $resultdir/$name.minus.bedgraph > $resultdir/$name.min.sort.bedgraph ## ===================> TODO add_minus_2_bedgraph("$resultdir/$name.min.sort.bedgraph"); $SOFTWAREDIR/kentUtils/kentUtils/bin/linux.x86_64/bedGraphToBigWig $resultdir/$name.plus.sort.bedgraph $genome.chromsizes $resultdir/$name.plus.bw $SOFTWAREDIR/kentUtils/kentUtils/bin/linux.x86_64/bedGraphToBigWig $resultdir/$name.min.sort.bedgraph $genome.chromsizes $resultdir/$name.min.bw else genomeCoverageBed -bg -ibam $resultdir/$name.sorted.bam -g $genome.chromsizes > $resultdir/$name.bedgraph sort -k1,1 -k2,2n $resultdir/$name.bedgraph > $resultdir/$name.sort.bedgraph $SOFTWAREDIR/kentUtils/kentUtils/bin/linux.x86_64/bedGraphToBigWig $resultdir/$name.sort.bedgraph $genome.chromsizes $resultdir/$name.bw fi } function featureCount { # for expression calling of genes or RNAs 2 programs can be used. Currently 'featureCounts' is the prefered methof printf "\n========== featureCount: "$name" ======================================\n" # OLD: $SOFTWAREDIR/molgentools_mirror/rnaseq/parse_coverage_v4.pl -depth $resultdir/$name -locus refseq -gff $genome.gff -gfftype \'$gff_type\' -o $resultdir/$name.AM.rpkm # -s = 0 (unstranded), 1 (stranded) and 2 (reversely stranded $SOFTWAREDIR/subread/bin/featureCounts -T $cpu -s 1 -M -t 'CDS' -g 'Name' -a $genome.gff -o $resultdir/$name.name.counts $resultdir/$name.sorted.bam # $SOFTWAREDIR/subread/bin/featureCounts -T $cpu -s 1 -M -t 'gene' -g 'locus_tag' -a $genome.gff -o $resultdir/$name.locus_tag.counts $resultdir/$name.sorted.bam # $SOFTWAREDIR/subread/bin/featureCounts -T $cpu -s 1 -M -t 'gene' -g 'Name' -a $genome.gff -o $resultdir/$name.name.counts $resultdir/$name.sorted.bam # $SOFTWAREDIR/subread/bin/featureCounts -T $cpu -s 1 -M -t 'tRNA' -g 'product' -a $genome.gff -o $resultdir/$name.tRNA.counts $resultdir/$name.sorted.bam # $SOFTWAREDIR/subread/bin/featureCounts -T $cpu -s 1 -M -t 'rRNA' -g 'product' -a $genome.gff -o $resultdir/$name.rRNA.counts $resultdir/$name.sorted.bam } function my_bowtie2_build { # build Bowties index files module load Bowtie2/2.3.4.3-foss-2018a bowtie2-build /data/p127804/PROJECTS/JanMaarten/genomes/D3.fna /data/p127804/PROJECTS/JanMaarten/genomes/D3 bowtie2-build /data/p127804/PROJECTS/JanMaarten/genomes/D17.fna /data/p127804/PROJECTS/JanMaarten/genomes/D17 } function my_chromsizes { # Chromsizes export PERL5LIB=$SOFTWAREDIR/molgentools_mirror/lib module load BioPerl/1.6.924-intel-2016a-Perl-5.20.3 $SOFTWAREDIR/molgentools_mirror/rnaseq/make_chromsizes.pl -i /data/pg-molgen/barbara/RNA-Seq/genomes/MG1363_24CH.fna -o /data/pg-molgen/barbara/RNA-Seq/genomes/MG1363_24CH.chromsizes } function MergeTables { # -s output folder [default=current folder] # -tables folder with tables to merge # -regex regular expression the select filenames # -columns comma separed list if columns to parse to the output table e.g. 0,3,7 $SOFTWAREDIR/molgentools_mirror/tables/merge_tables.pl -s $resultdir -tables $resultdir -regex 'locus_tag.counts$' -columns '0,6' -o 00.merged.counts.locus_tag $SOFTWAREDIR/molgentools_mirror/tables/merge_tables.pl -s $resultdir -tables $resultdir -regex 'ID.counts$' -columns '0,6' -o 00.merged.counts.ID $SOFTWAREDIR/molgentools_mirror/tables/merge_tables.pl -s $resultdir -tables $resultdir -regex 'name.counts$' -columns '0,6' -o 00.merged.counts.Name $SOFTWAREDIR/molgentools_mirror/tables/merge_tables.pl -s $resultdir -tables $resultdir -regex 'tRNA.counts$' -columns '0,6' -o 00.merged.counts.tRNA $SOFTWAREDIR/molgentools_mirror/tables/merge_tables.pl -s $resultdir -tables $resultdir -regex 'rRNA.counts$' -columns '0,6' -o 00.merged.counts.rRNA # clean the table from low values #$SOFTWAREDIR/molgentools_mirror/tables/remove_sum_zero.pl 00.merged.table 00.merged_clean_sum100.table } function my_mkdir { if [ ! -d $1 ]; then mkdir $1 ; fi } # ============================================================================================================================================================= # == Main script == # ============================================================================================================================================================= # 1. before starting build the index files of the genome(s) # a. Upload the .fna file # b. Upload the .gff file # c. Check if the header(s) of the .fna file is exactly matching the first column of the .gff file # d. Build the bowtie index ( run only once ) # my_bowtie2_build # 2. folders datadir=/data/p127804/PROJECTS/JanMaarten/primbio_202005 resultdir=/data/p127804/PROJECTS/JanMaarten/results_202005_repeatset3 genomedir=/data/p127804/PROJECTS/JanMaarten/genomes # 3. Submit all samples to the main function # the main function needs 3 parameters; Sample_name Foward_Reads Genome_base_Name function main { # this is the main function for processing each sample name=$1 forward=$2 genome=$3 forward_trimmed=$name.forward.fq # trimmomatic # my_fastqc # my_bowtie2 # process_bam # #CoverageBed featureCount } Bs168=/data/p127804/PROJECTS/JanMaarten/genomes/Bacillus_subtilis_subsp_subtilis_str_168 D3=/data/p127804/PROJECTS/JanMaarten/genomes/D3 D17=/data/p127804/PROJECTS/JanMaarten/genomes/D17 main D3A R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4239_D3A_IonXpressRNA_001.fastq $D3 main D3B R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4240_D3B_IonXpressRNA_002.fastq $D3 main D3C R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4241_D3C_IonXpressRNA_003.fastq $D3 main D17A R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4242_D17A_IonXpressRNA_004.fastq $D17 main D17B R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4243_D17B_IonXpressRNA_005.fastq $D17 main D17C R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4244_D17C_IonXpressRNA_006.fastq $D17 main D3_CIPA R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4245_D3_CIP_A_IonXpressRNA_007.fastq $D3 main D3_CIPB R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4246_D3_CIP_B_IonXpressRNA_008.fastq $D3 main D3_CIPC R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4247_D3_CIP_C_IonXpressRNA_009.fastq $D3 main D17_CIPA R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4248_D17_CIP_A_IonXpressRNA_010.fastq $D17 main D17_CIPB R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4249_D17_CIP_B_IonXpressRNA_011.fastq $D17 main D17_CIPC R_2020_04_30_14_36_59_user_GSS5PR-0198-76-Dijl_12samples_4250_D17_CIP_C_IonXpressRNA_012.fastq $D17 # Bowtie2 report # cat $resultdir/02.*.bowtie2.log > $resultdir/02.bowtie2_all.log # $SOFTWAREDIR/molgentools_mirror/rnaseq/RNAseq_bowtie_report_peregrine.pl -i $resultdir/02.bowtie2_all.log > $resultdir/02.bowtie2_all.log.table ## Manual option for Bowtie2 report ## SOFTWAREDIR=/data/pg-molgen/software ## cat 02.*.bowtie2.log > 02.bowtie2_all.log ## $SOFTWAREDIR/molgentools_mirror/rnaseq/RNAseq_bowtie_report_peregrine.pl -i 02.bowtie2_all.log > 02.bowtie2_all.log.table # 4. Merge Count Tables MergeTables