RNASeqpipeline_MULTIPLE.sh

u#!/bin/bash

SECONDS=0

# change working directory
cd /home/aysegul/RNAseqAnalysis/RNASeq_pipeline/


# STEP 1: Run fastqc
#fastqc RNAseqAnalysis/fastq_files/*.fq.gz -o RNAseqAnalysis/fastq_files

# run trimmomatic to trim reads with poor quality

for f in ‘ls -1 *_1.fq.gz | sed 's/\_1.fq.gz//'’; do echo java -jar ~/apps/Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 $f\_1.fq.gz $f\_2.fq.gz $f\_1_paired.fq.gz $f\_1_unpaired.fq.gz $f\_2_paired.fq.gz $f\_2_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 >> cmd_file; done
#echo "Trimmomatic finished running!"



# STEP 2: Run HISAT2
# mkdir HISAT2
# get the genome indices
# wget https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz


# run alignment
#for f in ‘ls -1 *_1_paired*.fq.gz | sed 's/\_1_paired.fq.gz//'’; do echo hisat2 -q --rna-strandness FR -x RNAseqAnalysis/HISAT2/grch38/genome -1 ${f}_1_paired.fq.gz -2  ${f}_2_paired.fq.gz | samtools sort -o ${f}.bam
>> cmd_file; done
#echo "HISAT2 finished running!"



# STEP 3: Run featureCounts - Quantification
# get gtf
# wget http://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.gtf.gz
#mkdir bams file

#featureCounts -p -s 2 -T 8 -a Homo_sapiens.GRCh38.110.gtf -o count.txt bams/*.bam
#count.out

duration=$SECONDS
echo "$(($duration / 60)) minutes and $(($duration % 60)) seconds elapsed."