This GitHub repository contains a Bash script that contains a series of code blocks that automate key steps in RNA sequencing data analysis, along with additional Python and R scripts for downstream analysis.
-
Sequencing read data in the fastq.gz format
-
Index files for the reference genome of interest, in this case Human Genome hg38
-
Ideally perform your own indexing using HISAT2
-
A .gtf file of annotated features for your indexed genome
-
Splice Site file for your indexed genome to improve alignment accuracy across exon-exon boundaries
-
For the example analysis below the GENCODE Release 49 (GRCh38.p14) of the human genome was used https://www.gencodegenes.org/human/
-
GCA_000001405.15_GRCh38_no_alt_analysis_set obtained from NCBI was used for an index and the corresponding gtf file gencode.v49.annotation.gtf was obtained from GENCODE
-
Splice Site file was generated by using a python script packaged with HISAT2 using the following code:
hisat2_extract_splice_sites.py gencode.v49.annotation.gtf > splice_sites.txtExample repository structure
n.b. some folders are not described as these numbers will be used later for the steps involving deduplication with UMI-tools.
.
├── rnaseq_analysis.sh
├── mergeFeatureCounts.py
├── 00_raw_fastq
├── 01_fastqc_pretrim
├── 02_trimmed_fastq
├── 03_fastqc_posttrim
├── 04
├── 05
├── 06_sam
├── 07_bam
├── 08_bam_sorted
├── 09
├── 10_featurecounts
├── logs
├── multiqc_reports
└── README.md
To begin:
- Add your paired .fastq files to ./00_raw_fastq
- Follow the configuration steps below
- Give rnaseq_analysis.sh execution permissions and run the script from its location in your terminal
This section allows you to control some of the sections of the script including:
- Number of CPU threads used where this can be set by specific packages
- Location and names of the index files
- Directory variables
- Short code block to make the directories if the are not available already
- A control switch to set which step the script starts from
The purpose of this step is to ascertain what level of processing the samples have already undergone e.g. adaptor removal, and whether they need trimming with step 2 below.
if [ "$START_STEP" -le 1 ]; then
echo "=== Step 1: FastQC BEFORE trimming ==="
fastqc -t $THREADS -o $PRETRIM_QC $RAW_DIR/*.fastq.gz
multiqc -o $MULTIQC_DIR $PRETRIM_QC
fiSamples which require trimming are trimmed with fastp. The fastp tool generates its own reports, but FastQC is also run on the post-trimmed samples to directly compare to pre-trimmed results.
if [ "$START_STEP" -le 2 ]; then
echo "=== Step 2: Trimming with fastp ==="
for fq1 in $RAW_DIR/*_R1_001.fastq.gz; do # Change this based on the format of your file
fq2=${fq1/_R1_001.fastq.gz/_R2_001.fastq.gz}
base1=$(basename "$fq1" .fastq.gz)
base2=$(basename "$fq2" .fastq.gz)
sample_prefix=${base1%_R1_001}
fastp -i "$fq1" -I "$fq2" \
-o "$TRIM_DIR/${base1}_val_1.fq.gz" \
-O "$TRIM_DIR/${base2}_val_2.fq.gz" \
-w $THREADS \
--detect_adapter_for_pe \
--length_required 30 \
--html "$LOG_DIR/${sample_prefix}_fastp.html" \
--json "$LOG_DIR/${sample_prefix}_fastp.json"
done
fastqc -t $THREADS -o $POSTTRIM_QC $TRIM_DIR/*.fq.gz
multiqc -o $MULTIQC_DIR $POSTTRIM_QC
fiReads are aligned to the human genome using the index files as described above. This generates SAM files and a report of the read mapping statistics.
if [ "$START_STEP" -le 3 ]; then
echo "=== Step 3: Mapping with HISAT2 ==="
for fq1 in $TRIM_DIR/*_val_1.fq.gz; do
fq2=${fq1/_R1_001_val_1.fq.gz/_R2_001_val_2.fq.gz}
sample=$(basename "$fq1" | sed 's/_R1_001_val_1.fq.gz//')
hisat2 -x $GENOME_INDEX \
--known-splicesite-infile $SPLICE_SITES \
-p $THREADS \
-1 "$fq1" \
-2 "$fq2" \
-S "$SAM_DIR/${sample}.sam" \
2> "$LOG_DIR/${sample}_hisat2.log"
done
fiSAM files are converted to BAM files and indexed (step 5) so these files can be used to feature annotate the reads locations to corresponding genes (step 6).
if [ "$START_STEP" -le 4 ]; then
echo "=== Step 4: SAM to BAM ==="
for sam in $SAM_DIR/*.sam; do
bam=$BAM_DIR/$(basename "$sam" .sam).bam
samtools view -bS "$sam" > "$bam"
done
fiif [ "$START_STEP" -le 5 ]; then
echo "=== Step 5: Sort & Index BAM ==="
for bam in $BAM_DIR/*.bam; do
sorted=$SORTED_BAM_DIR/$(basename "$bam" .bam)_sorted.bam
samtools sort "$bam" -o "$sorted"
samtools index "$sorted"
done
fiif [ "$START_STEP" -le 6 ]; then
echo "=== Step 6: featureCounts ==="
for sorted_bam_file in $SORTED_BAM_DIR/*_sorted.bam; do
sample=$(basename "$sorted_bam_file" _sorted.bam)
featureCounts -T $THREADS -p --countReadPairs \
-t exon -a "$GTF_FILE" \
-o "$COUNTS_DIR/${sample}_featurecounts.txt" \
-g gene_id "$sorted_bam_file"
done
fi- Script:
merge_featureCounts.py - Usage:
python merge_featureCounts.py ./10_featurecounts/*.txt mergedCounts.csv- Script:
DSeq2_analysis.R - Usage:
Execute in an R environment
Publicly available example data were obtained from the European Nucleotide Archive (ENA) under project accession PRJNA1007680.
In this study, the authors used the human prostate cancer cell line PC3 and treated cells with either RSM3 (a peptide inhibitor targeting METTL3) or a control agent. Following RNA isolation, paired-end sequencing libraries were generated and sequenced on an Illumina NovaSeq 6000.
Initial quality control (QC) indicated that the sequencing reads had likely already been preprocessed using fastp, as described in the corresponding GEO accession GSE241315. Therefore, the read preprocessing step (Step 2) was omitted from this pipeline.
A Conda environment YAML file of the required packages Option for single end (SE) reads Option to include de-duplication where umis have been included Further QC steps e.g. BAM file QC with Qualimap