RiboParser 0.1.14
A pipeline for ribosome profiling data analysis
Requires Python: >=3.12

<!--
 * @Author: 'rensc' 'rensc0718@163.com'
 * @Date: 2024-10-15 11:44:58
 * @LastEditors: 'rensc' 'rensc0718@163.com'
 * @LastEditTime: 2025-02-18 03:50
 * @FilePath: \RiboParser\README.md 
 * 
-->


# RiboParser 

```
Ren, S., Li, Y. & Zhou, Z. 
RiboParser/RiboShiny: An integrated platform for comprehensive analysis and visualization of ribo-seq data. 
Journal of Genetics and Genomics (2025) 
doi:10.1016/j.jgg.2025.04.010.
```

To streamline understanding and application, we will analyze publicly accessible project data, breaking down each analytical step to illustrate the complete workflow.

This process encompasses both general analysis steps and specialized analysis and visualization techniques, facilitated by `RiboParser` and `RiboShiny`.

The specific steps involved are:

1. Software installation
2. Reference file creation
3. Raw data download
4. Raw data cleaning
5. Data alignment
6. Sequencing quality analysis
7. Gene-level analysis
8. Codon-level analysis

The results of this data analysis can be further analyzed and visualized in `RiboShiny`.



## 1. Software configuration

First, you need to install the software of `miniconda` or `micromamba`.

### 1.1 create the environment with `conda`

```bash
conda create -n ribo
conda activate ribo
```

### 1.2 Install `riboparser` and the entire environment directly via `miniconda` or `micromamba`.
```bash
# conda
conda install riboparser -c rensc

# mamba
micromamba install riboparser -c rensc

```


### 1.3 Install `riboparser` and software dependencies step by step.

Install the software dependencies.

```bash
conda install bowtie samtools cutadapt star bedtools subread rsem gffread sra-tools \
 ucsc-genepredtogtf ucsc-gtftogenepred ucsc-gff3togenepred ucsc-bedgraphtobigwig ucsc-bedsort \
 -c bioconda

conda install pigz -c conda-forge
```


`pip` install `RiboParser`

```bash
pip install riboparser
```

Alternatively, we can download the version from GitHub, re-setup and then install it.

```bash
cd RiboParser
# install the dependency
pip install build

# build the riboparser package
python -m build

# install the riboparser
pip install .

```

### 1.4 run the test
Test software for dependency, installation, and operation issues.

```bash
rpf_Check -h
rpf_CST -h
```

## 2. Prepare reference files

### 2.1 An example of the complete project directory is shown below

The complete data analysis includes reference preparation, rawdata, RNA-seq data analysis, and Ribo-seq data analysis.

```bash
$ cd && cd ./sce/
$ tree -d

.
├── 1.reference
│   ├── genome
│   ├── mrna
│   ├── ncrna
│   ├── norm
│   ├── rrna
│   ├── rsem-index
│   ├── star-index
│   └── trna
├── 2.rawdata
│   ├── ribo-seq
│   └── rna-seq
├── 3.rna-seq
│   ├── 1.cleandata
│   ├── 2.bowtie
│   ├── 3.star
│   ├── 4.quantification
│   └── 5.riboparser
│       ├── 01.qc
│       ├── 02.digestion
│       ├── 03.offset
│       ├── 04.density
│       ├── 05.merge
│       ├── 06.periodicity
│       ├── 07.metaplot
│       ├── 08.coverage
│       ├── 09.correlation
│       ├── 10.shuffle
│       └── 11.retrieve
└── 4.ribo-seq
    ├── 1.cleandata
    ├── 2.bowtie
    ├── 3.star
    ├── 4.quantification
    └── 5.riboparser
        ├── 01.qc
        ├── 02.digestion
        ├── 03.offset
        ├── 04.density
        ├── 05.merge
        ├── 06.periodicity
        ├── 07.metaplot
        ├── 08.coverage
        ├── 09.correlation
        ├── 10.quantification
        ├── 11.pausing_score
        ├── 12.codon_occupancy
        ├── 13.codon_decoding_time
        ├── 14.codon_selection_time
        ├── 15.coefficient_of_variation
        ├── 16.meta_codon
        ├── 17.shuffle
        ├── 18.retrieve
        └── 19.frame_shift
```

### 2.2 Prepare the reference genome index
#### 2.2.1 Create the directory

Create folders to hold different types of reference sequence files.

```bash
$ mkdir -p ./sce/1.reference/
$ cd ./sce/1.reference/
$ mkdir cdna genome gtf mrna ncrna rrna trna norm rsem-index
```

#### 2.2.2 Download reference files from NCBI

Use the most common data analysis file format, the genome sequence in fasta format, and the reference file in GTF or GFF3 format.

```bash
# genome sequence
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz

# GTF or GFF3
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.gtf.gz
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.gff.gz

# cDNA sequence
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_rna.fna.gz

# feature table
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_feature_table.txt.gz

# decompression
$ gunzip *.gz

$ gffread -g GCF_000146045.2_R64_genomic.fna GCF_000146045.2_R64_genomic.gff -F -w cdna.fa
```

#### 2.2.3 Create the `genome` index using `bowtie`

```bash
$ bowtie-build ../GCF_000146045.2_R64_genomic.fna ./genome/genome --threads 12 &>> ./genome/genome_build.log
```

#### 2.2.4 Create an `mRNA` index using `bowtie`

A custom script is used here to extract the corresponding sequence information 
from the `fasta` file based on the sequence name.

```bash
$ retrieve_seq -h

usage: retrieve_seq [-h] [-v] -i INPUT -n NAME [-u UNMAPPED] -o OUTPUT

This script is used to retrieve the fasta sequence by name.

options:
  -h, --help     show this help message and exit
  -v, --version  show programs version number and exit
  -i INPUT       input the fasta file
  -n NAME        gene ids in txt format
  -u UNMAPPED    output the unmapped gene ids
  -o OUTPUT      prefix of output file name (default results_peaks.txt)
```


```bash
# filter the mrna sequence
$ grep -i 'gbkey=mRNA' ./cdna.fa | cut -d ' ' -f 1 | cut -c 2- > ./mrna/mrna.ids
$ retrieve_seq -i ./cdna.fa -n ./mrna/mrna.ids -o ./mrna/mrna.fa &>> ./mrna/mrna_build.log
# build the mrna index
$ bowtie-build ./mrna/mrna.fa ./mrna/mrna --threads 12 &>> ./mrna/mrna_build.log
```

#### 2.2.5 Create an `rRNA` index using `bowtie`
```bash
# filter the rrna sequence
$ grep -i 'gbkey=rRNA' ./cdna.fa | cut -d ' ' -f 1 | cut -c 2- > ./rrna/rrna.ids
$ retrieve_seq -i ./cdna.fa -n ./rrna/rrna.ids -o ./rrna/rrna.fa &>> ./rrna/rrna_build.log
# build the rrna index
$ bowtie-build ./rrna/rrna.fa ./rrna/rrna --threads 12 &>> ./rrna/rrna_build.log
```

#### 2.2.6 Create an `tRNA` index using `bowtie`
```bash
# filter the trna sequence
$ grep -i 'gbkey=tRNA' ./cdna.fa | cut -d ' ' -f 1 | cut -c 2- > ./trna/trna.ids
$ retrieve_seq -i ./cdna.fa -n ./trna/trna.ids -o ./trna/trna.fa &>> ./trna/trna_build.log
# build the trna index
$ bowtie-build ./trna/trna.fa ./trna/trna --threads 12 &>> ./trna/trna_build.log
```


#### 2.2.7 Create an `ncRNA` index using `bowtie`
```bash
# filter the ncrna sequence
$ grep -iE 'gbkey=ncRNA|gbkey=lnc_RNA|gbkey=miRNA|gbkey=snoRNA|gbkey=snRNA|gbkey=misc_RNA' ./cdna.fa | cut -d ' ' -f 1 | cut -c 2- > ./ncrna/ncrna.ids
$ retrieve_seq -i ./cdna.fa -n ./ncrna/ncrna.ids -o ./ncrna/ncrna.fa &>> ./ncrna/ncrna_build.log
# build the ncrna index
$ bowtie-build ./ncrna/ncrna.fa ./ncrna/ncrna --threads 12 &>> ./ncrna/ncrna_build.log
```

#### 2.2.8 Standardized `gtf` or `gff3` files

- Explanation of `rpf_Reference` 

```bash
$ rpf_Reference -h

usage: rpf_Reference [-h] -g GENOME -t GTF -o OUTPUT [-u UTR] [-c] [-l] [-w]

This script is used to build the references for the RiboParser.

options:
  -h, --help  show this help message and exit
  -u UTR      add the pseudo UTR to the leaderless transcripts (default: 0 nt).
  -c          only retain the protein coding transcripts (default: False).
  -l          only retain the longest protein coding transcripts, it is recommended to select 
              the longest transcript for subsequent analysis. (default: False).
  -w          output whole message (default: False).

Required arguments:
  -g GENOME   the input file name of genome sequence
  -t GTF      the input file name of gtf file
  -o OUTPUT   the prefix of output file. (prefix + _norm.gtf)
```

- create the references from GTF and Genome fasta

```bash
$ rpf_Reference \
 -g ../GCF_000146045.2_R64_genomic.fna \
 -t ../GCF_000146045.2_R64_genomic.gff \
 -u 30 -o ./norm/gene &>> ./norm/norm_build.log
```

#### 2.2.9 Create a `genome` index using `STAR`

```bash
$ STAR \
 --genomeSAindexNbases 11 \
 --runThreadN 12 \
 --runMode genomeGenerate \
 --genomeDir ./star-index \
 --genomeFastaFiles GCF_000146045.2_R64_genomic.fna \
 --sjdbGTFfile ./norm/gene.norm.gtf
```

#### 2.2.10 Create a `transcriptome` index using `rsem`

```bash
$ rsem-prepare-reference \
 -p 12 \
 --gtf ../norm/gene.norm.gtf ../GCF_000146045.2_R64_genomic.fna ./rsem-index/rsem
```


## 3. Data preprocessing and alignment

In order to introduce the analysis process and usage of `RiboParser`, RNA-seq and Ribo-seq data from dataset `GSE67387` are used as examples here.

```shell
# dataset
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67387

# reference
Nedialkova DD, Leidel SA. Optimization of Codon Translation Rates via tRNA Modifications Maintains Proteome Integrity. Cell 2015 Jun 18;161(7):1606-18. 
PMID: 26052047
```


### 3.1 Fundamental analysis of RNA-seq data in `GSE67387` dataset

#### 3.1.1 Download RNA-seq raw data

Use `prefetch` in `sra-tools` to download the raw SRA-format data and extract it into `fastq` format files.

```bash
$ mkdir -p ./sce/2.rawdata/rna-seq/
$ cd ./sce/2.rawdata/rna-seq/

#################################################
# download rna-seq
$ prefetch -o SRR1944925.sra SRR1944925
$ prefetch -o SRR1944926.sra SRR1944926
$ prefetch -o SRR1944927.sra SRR1944927
$ prefetch -o SRR1944928.sra SRR1944928
$ prefetch -o SRR1944929.sra SRR1944929
$ prefetch -o SRR1944930.sra SRR1944930
$ prefetch -o SRR1944931.sra SRR1944931
$ prefetch -o SRR1944932.sra SRR1944932
$ prefetch -o SRR1944933.sra SRR1944933
$ prefetch -o SRR1944934.sra SRR1944934
$ prefetch -o SRR1944935.sra SRR1944935

# decompression
for sra in *.sra
do

fastq-dump $sra
pigz *fastq

done
```

#### 3.1.2 RNA-seq data cleaning

Because the data from the gse project is cleaned, it does not include adapter and index sequences. So the following is just to show the general steps, do not need to run.

1. RNA-seq data cleaning

```bash
$ cd
$ mkdir -p ./sce/3.rna-seq/1.cleandata/
$ cd ./sce/3.rna-seq/1.cleandata/

#################################################
# run the cutadapt
for fq in ../../2.rawdata/rna-seq/*fastq.gz
do
cutadapt --match-read-wildcards \
 -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGC \
 -m 25 -O 6 -j 12 \
 -o `\basename $fq fastq.gz`clean.fastq.gz $fq &>> $fq".log"

done
```

#### 3.1.3 Align clean data to different types of reference files

To assess library quality and eliminate the impact of reads originating from different non-coding RNAs (ncRNAs) on subsequent analysis, we employed `bowtie` to classify the reads form sequencing data.

Under normal circumstances, especially for RNA-seq libraries constructed using the oligo(dT) method, most reads originate from mRNA. Therefore, for RNA-seq analysis, this step is generally not necessary. It is suitable for use in libraries constructed by the rRNA-depletion method.

1. Align the RNA-seq data to references

```bash
$ cd
$ mkdir -p ./sce/3.rna-seq/2.bowtie/
$ cd ./sce/3.rna-seq/2.bowtie/

#################################################
# set database
rrna='../../sce/1.reference/rrna/rrna'
trna='../../sce/1.reference/trna/trna'
ncrna='../../sce/1.reference/ncrna/ncrna'
mrna='../../sce/1.reference/mrna/mrna'
chrom='../../sce/1.reference/genome/genome'

threads=12
mismatch=1

# alignment reads to reference
for fq in ../1.cleandata/*fastq.gz
do
fqname=`\basename $fq .fastq.gz`

## rrna
bowtie -p $threads -v $mismatch --un="$fqname".norrna.fq --al="$fqname".rrna.fq \
 -x $rrna $fq -S "$fqname".rrna.sam 2>> "$fqname".log

## trna
bowtie -p $threads -v $mismatch --un="$fqname".notrna.fq --al="$fqname".trna.fq \
 -x $trna "$fqname".norrna.fq -S "$fqname".trna.sam 2>> "$fqname".log

## ncrna
bowtie -p $threads -v $mismatch --un="$fqname".noncrna.fq --al="$fqname".ncrna.fq \
 -x $ncrna "$fqname".notrna.fq -S "$fqname".ncrna.sam 2>> "$fqname".log

## mrna
bowtie -p $threads -v $mismatch --un="$fqname".nomrna.fq --al="$fqname".mrna.fq \
 -x $mrna "$fqname".noncrna.fq -S "$fqname".mrna.sam 2>> "$fqname".log

## genome
bowtie -p $threads -v $mismatch --un="$fqname".nogenome.fq --al="$fqname".genome.fq \
 -x $chrom "$fqname".nomrna.fq -S "$fqname".genome.sam 2>> "$fqname".log

## compress fastq
pigz *fq

## compress sam
for sam in *.sam
do

samtools view -h -F 4 $sam | samtools sort -@ $threads -o `\basename $sam sam`bam
rm $sam

done

done
```


2. Statistical alignment results for all references.

- Explanation of `merge_bwt_log`

```bash
$ merge_bwt_log -h

Step1: Checking the input Arguments.
usage: merge_bwt_log [-h] -l LIST [LIST ...] -o OUTPUT [-n NAME]

This script is used to statstic the mapped reads from log files

options:
  -h, --help            show this help message and exit
  -n NAME, --name NAME  set the name of each database (default: rRNA,tRNA,ncRNA,mRNA,Genome).

Required arguments:
  -l LIST [LIST ...], --list LIST [LIST ...]
                        List for bowtie mapping log files (e.g., '*log').
  -o OUTPUT             prefix of output file name.
```

- Statistical alignment results

```bash
#################################################
# merge all log files
merge_bwt_log -n rRNA,tRNA,ncRNA,mRNA,Genome -l *log -o RNA_seq &>> merge_bowtie.log
```

#### 3.1.4 Aligning mRNA reads using `STAR`

Following the removal of ncRNA reads, the remaining clean reads were realigned to the yeast genome using `STAR`.

1. Aligning mRNA reads (RNA-seq) using `STAR`
```bash
$ cd
$ mkdir -p ./sce/3.rna-seq/3.star/
$ cd ./sce/3.rna-seq/3.star/

#################################################
# set the option and database
genome='../../1.reference/star-index/'
threads=12

#################################################
# map the all rna-seq reads to genome and transcriptome region
for fastq in ../2.bowtie/*.noncrna.fq.gz
do

## get file name
output=$(basename $fastq .noncrna.fq.gz)

#################################################
## run the alignment
STAR --runThreadN $threads \
 --readFilesCommand zcat \
 --genomeDir $genome \
 --readFilesIn $fastq \
 --outFileNamePrefix $output \
 --outSAMtype BAM Unsorted \
 --outFilterType BySJout \
 --quantMode TranscriptomeSAM GeneCounts \
 --outReadsUnmapped Fastx \
 --outSAMattributes All \
 --alignEndsType Local \
 --outFilterMultimapNmax 3 \
 --outFilterMismatchNmax 1 \
 --alignIntronMax 10000 \
 --outFilterMatchNmin 20
# --outWigType wiggle --outWigNorm RPM

pigz *mate1

#################################################
## sort the bam file
samtools sort -@ $threads $output"Aligned.out.bam" -o $output"Aligned.sortedByCoord.out.bam"
samtools index -@ $threads $output"Aligned.sortedByCoord.out.bam"
rm $output"Aligned.out.bam"

done
```

#### 3.1.5 Estimating gene expression levels with either `RSEM` or `featureCounts`

Both `RSEM` and `featureCounts` can be employed to quantify gene expression levels. For the purpose of this analysis, we will utilize `RSEM` as a representative tool.

1. Estimating transcript abundance using RNA-seq data

```bash
$ cd
$ mkdir -p ./sce/3.rna-seq/4.quantification/
$ cd ./sce/3.rna-seq/4.quantification/

#################################################
# quantify the gene expression
for bam in ../3.star/*Aligned.toTranscriptome.out.bam
do
rsem-calculate-expression -p 10 \
 --no-bam-output --alignments \
 -q $bam ../../1.reference/rsem-index/rsem `\basename $bam Aligned.toTranscriptome.out.bam`
# rsem-calculate-expression -p 10 \
# --paired-end --no-bam-output --alignments \
# -q $bam ../../1.reference/rsem-index/rsem `\basename $bam Aligned.toTranscriptome.out.bam`

done
```

2. Integrating RNA-seq quantification values for all samples

- Explanation of `merge_rsem`

```bash
$ merge_rsem -h

usage: merge_rsem [-h] -l LIST [LIST ...] -o OUTPUT [-c {expected_count,TPM,FPKM}]

This script is used to merge specified columns from result files

options:
  -h, --help            show this help message and exit
  -c {expected_count,TPM,FPKM}, --column {expected_count,TPM,FPKM}
                        Column name to merge (e.g., 'expected_count').

Required arguments:
  -l LIST [LIST ...], --list LIST [LIST ...]
                        List for result files (e.g., '*results').
  -o OUTPUT             output file name.
```

- Integrating RNA-seq quantification values

```bash
#################################################
# merge the gene expression
merge_rsem -c expected_count -l *.genes.results -o gene.expected_count.txt &>> merge_rsem.log
merge_rsem -c TPM -l *.genes.results -o gene.TPM.txt &>> merge_rsem.log
merge_rsem -c FPKM -l *.genes.results -o gene.FPKM.txt &>> merge_rsem.log

#################################################
# merge the isoforms expression
merge_rsem -c expected_count -l *.isoforms.results -o isoforms.expected_count.txt &>> merge_rsem.log
merge_rsem -c TPM -l *.isoforms.results -o isoforms.TPM.txt &>> merge_rsem.log
merge_rsem -c FPKM -l *.isoforms.results -o isoforms.FPKM.txt &>> merge_rsem.log
```


### 3.2 Fundamental analysis of Ribo-seq data in `GSE67387` dataset

#### 3.2.1 Download Ribo-seq raw data

Use `prefetch` in `sra-tools` to download the raw SRA-format data and extract it into `fastq` format files.

```bash
$ cd
$ mkdir -p ./sce/2.rawdata/ribo-seq/
$ cd ./sce/2.rawdata/ribo-seq/

#################################################
# download ribo-seq
prefetch -o SRR1944912.sra SRR1944912
prefetch -o SRR1944913.sra SRR1944913
prefetch -o SRR1944914.sra SRR1944914
prefetch -o SRR1944915.sra SRR1944915
prefetch -o SRR1944916.sra SRR1944916
prefetch -o SRR1944917.sra SRR1944917
prefetch -o SRR1944918.sra SRR1944918
prefetch -o SRR1944919.sra SRR1944919
prefetch -o SRR1944920.sra SRR1944920
prefetch -o SRR1944921.sra SRR1944921
prefetch -o SRR1944922.sra SRR1944922
prefetch -o SRR1944923.sra SRR1944923

# decompression
for sra in *.sra
do

fastq-dump $sra
pigz *fastq

done
```

#### 3.2.2 Ribo-seq data cleaning

```bash
$ cd
$ mkdir -p ./sce/4.ribo-seq/1.cleandata/
$ cd ./sce/4.ribo-seq/1.cleandata/

#################################################
# run the cutadapt
for fq in ../../2.rawdata/ribo-seq/*fastq.gz
do
cutadapt --match-read-wildcards \
 -a AAAAAAAA \
 -m 25 -O 6 -j 10 \
 -o `\basename $fq fastq.gz`clean.fastq.gz $fq &>> $fq".log"
done
```


#### 3.2.3 Align clean data to different types of reference files

To assess library quality and eliminate the impact of reads originating from different non-coding RNAs (ncRNAs) on subsequent analysis, we employed `bowtie` to classify the reads form sequencing data.

1. Aligning Ribo-seq data
```bash
$ cd
$ mkdir -p ./sce/4.ribo-seq/2.bowtie/
$ cd ./sce/4.ribo-seq/2.bowtie/

#################################################
# set database
rrna='../../sce/1.reference/rrna/rrna'
trna='../../sce/1.reference/trna/trna'
ncrna='../../sce/1.reference/ncrna/ncrna'
mrna='../../sce/1.reference/mrna/mrna'
chrom='../../sce/1.reference/genome/genome'

threads=12
mismatch=1

# alignment reads to reference
for fq in ../1.cleandata/*fastq.gz
do
fqname=`\basename $fq .fastq.gz`

## rrna
bowtie -p $threads -v $mismatch --un="$fqname".norrna.fq --al="$fqname".rrna.fq \
 -x $rrna $fq -S "$fqname".rrna.sam 2>> "$fqname".log

## trna
bowtie -p $threads -v $mismatch --un="$fqname".notrna.fq --al="$fqname".trna.fq \
 -x $trna "$fqname".norrna.fq -S "$fqname".trna.sam 2>> "$fqname".log

## ncrna
bowtie -p $threads -v $mismatch --un="$fqname".noncrna.fq --al="$fqname".ncrna.fq \
 -x $ncrna "$fqname".notrna.fq -S "$fqname".ncrna.sam 2>> "$fqname".log

## mrna
bowtie -p $threads -v $mismatch --un="$fqname".nomrna.fq --al="$fqname".mrna.fq \
 -x $mrna "$fqname".noncrna.fq -S "$fqname".mrna.sam 2>> "$fqname".log

## genome
bowtie -p $threads -v $mismatch --un="$fqname".nogenome.fq --al="$fqname".genome.fq \
 -x $chrom "$fqname".nomrna.fq -S "$fqname".genome.sam 2>> "$fqname".log

## compress fastq
pigz *fq

## compress sam
for sam in *.sam
do

samtools view -h -F 4 $sam | samtools sort -@ $threads -o `\basename $sam sam`bam
rm $sam

done

done
```

2. Statistical alignment results for all databases.
```bash
#################################################
# merge all log files
merge_bwt_log -n rRNA,tRNA,ncRNA,mRNA,Genome -l *log -o sce &>> merge_bowtie.log
```


#### 3.2.4 Aligning mRNA reads (Ribo-seq) using `STAR`

Following the removal of ncRNA reads, the remaining clean reads were realigned to the yeast genome using `STAR`.

```bash
$ cd
$ mkdir -p ./sce/4.ribo-seq/3.star/
$ cd ./sce/4.ribo-seq/3.star/

#################################################
# set the option and database
genome='../../1.reference/star-index/'
threads=12

#################################################
# map the all rna-seq reads to genome and transcriptome region
for fastq in ../2.bowtie/*.noncrna.fq.gz
do

## get file name
output=$(basename $fastq .noncrna.fq.gz)

#################################################
## run the alignment
STAR --runThreadN $threads \
 --readFilesCommand zcat \
 --genomeDir $genome \
 --readFilesIn $fastq \
 --outFileNamePrefix $output \
 --outSAMtype BAM Unsorted \
 --outFilterType BySJout \
 --quantMode TranscriptomeSAM GeneCounts \
 --outReadsUnmapped Fastx \
 --outSAMattributes All \
 --alignEndsType Local \
 --outFilterMultimapNmax 3 \
 --outFilterMismatchNmax 1 \
 --alignIntronMax 10000 \
 --outFilterMatchNmin 20
# --outWigType wiggle --outWigNorm RPM

pigz *mate1

#################################################
## sort the bam file
samtools sort -@ $threads $output"Aligned.out.bam" -o $output"Aligned.sortedByCoord.out.bam"
samtools index -@ $threads $output"Aligned.sortedByCoord.out.bam"
rm $output"Aligned.out.bam"

done
```


#### 3.2.5 Estimating gene expression levels with either `RSEM` or `featureCounts`

Both `RSEM` and `featureCounts` can be employed to quantify gene expression levels. For the purpose of this analysis, we will utilize `RSEM` as a representative tool.

1. Estimating transcript abundance using Ribo-seq data

```bash
$ cd
$ mkdir -p ./sce/4.ribo-seq/4.quantification/
$ cd ./sce/4.ribo-seq/4.quantification/

#################################################
# quantify the isoforms expression
for bam in ../3.star/*Aligned.toTranscriptome.out.bam
do

rsem-calculate-expression -p 12 \
 --no-bam-output --alignments \
 -q $bam ../../1.reference/rsem-index/rsem `\basename $bam Aligned.toTranscriptome.out.bam`

done
```

2. Integrating Ribo-seq quantification values for all samples

```bash
#################################################
# merge the gene expression
merge_rsem -c expected_count -l *.genes.results -o gene.expected_count.txt &>> merge_rsem.log
merge_rsem -c TPM -l *.genes.results -o gene.TPM.txt &>> merge_rsem.log
merge_rsem -c FPKM -l *.genes.results -o gene.FPKM.txt &>> merge_rsem.log

#################################################
# merge the isoforms expression
merge_rsem -c expected_count -l *.isoforms.results -o isoforms.expected_count.txt &>> merge_rsem.log
merge_rsem -c TPM -l *.isoforms.results -o isoforms.TPM.txt &>> merge_rsem.log
merge_rsem -c FPKM -l *.isoforms.results -o isoforms.FPKM.txt &>> merge_rsem.log
```


## 4. Perform RNA-seq data analysis of `GSE67387` with `RiboParser`

### 4.0 Prepare the directory to store the results

```bash
$ cd
$ mkdir -p ./3.rna-seq/5.riboparser/
$ cd ./3.rna-seq/5.riboparser/
$ mkdir 01.qc 02.digestion 03.offset 04.density 05.merge \
 06.periodicity 07.metaplot 08.coverage 09.correlation 10.shuffle
```

### 4.1 Quality check of sequencing data

1. Quality check of RNA-seq data

```bash
$ cd ./01.qc

#################################################
# check the ribo-seq quality
for bam in ../../3.star/*Aligned.toTranscriptome.out.bam
do
prefix_name=$(basename $bam Aligned.toTranscriptome.out.bam)

rpf_Check -b $bam -s --thread 10 \
 -t ../../../1.reference/norm/gene.norm.txt \
 -o $prefix_name &>> $prefix_name".log"

done
```

2. Integrating RNA-seq quality check results for all samples
```bash

#################################################
# merge the rna-seq quality results
merge_length -l *length_distribution.txt -o sce
merge_saturation -l *gene_saturation.txt -o sce

cd ..
```

### 4.2 Enzymatic bias in NGS library preparation
1. Bias in restriction enzyme digestion and ligation in sequencing data

```bash
$ cd ./02.digestion/

#################################################
# check the reads digestion
for bam in ../01.qc/*.bam
do
prefix_name=$(basename $bam .bam)

rpf_Digest -b $bam -m 25 -M 50 --scale \
 -s ../../../1.reference/norm/gene.norm.rna.fa \
 -t ../../../1.reference/norm/gene.norm.txt \
 -o $prefix_name &>> $prefix_name".log"

done
```

2. Integrating reads digestion results for all samples

```bash
#################################################
# merge the rpf digestion
merge_digestion -l ./02.digestion/*pwm.txt -o sce

cd ..
```

### 4.3 Use `RiboParser` to create the `offset` table

1. Create the `offset` table for RNA-seq

`Offset` prediction is unnecessary for RNA-seq analysis. A constant `offset` of 12 can be assigned to all entries in the table.

```bash
$ cd ./03.offset/

#################################################
# set the offset table
for bam in ../01.qc/*.bam
do

prefix_name=$(basename $bam .bam)
rna_Offset -m 27 -M 50 -e 12 -o $prefix_name &>> $prefix_name".log"

done
```

### 4.4 Convert the `BAM` file to reads density

Transform the read counts in a `BAM` file into density values and save them in a `TXT` format file.

1. Transform RNA-seq data output

```bash
$ cd ./04.density/

#################################################
# convert the reads to density
for bam in ../01.qc/*.bam
do
prefix_name=$(basename $bam .bam)

rna_Density -b $bam -m 27 -M 33 -l --thread 10 \
 -p ../03.offset/$prefix_name"_offset.txt" \
 -s ../../../1.reference/norm/gene.norm.rna.fa \
 -t ../../../1.reference/norm/gene.norm.txt \
 -o $prefix_name &>> $prefix_name".log"

done

cd ..
```


### 4.5 Integrating density file for all samples

Data from various batches and samples can be integrated for unified analysis. If not integrated, individual analysis for each sample is required, increasing the number of operational steps.

1. Integrating RNA-seq density results for all samples
```bash
$ cd ./05.merge/

#################################################
# create the samples file: RNA.file.list
merge_dst_list -l ../04.density/*_rna.txt -o RNA.file.list

cat RNA.file.list

Name File  Type
wt_rna_YPD1 /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944924_rna.txt RNA
wt_rna_YPD2 /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944925_rna.txt RNA
wt_rna_YPD3 /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944926_rna.txt RNA
ncs2d_rna_YPD1  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944927rna.txt RNA
ncs2d_rna_YPD2  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944928_rna.txt RNA
ncs2d_rna_YPD3  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944929_rna.txt RNA
elp6d_rna_YPD1  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944930_rna.txt RNA
elp6d_rna_YPD2  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944931_rna.txt RNA
elp6d_rna_YPD3  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944932_rna.txt RNA
ncs2d_elp6d_rna_YPD1  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944933_rna.txt RNA
ncs2d_elp6d_rna_YPD2  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944934_rna.txt RNA
ncs2d_elp6d_rna_YPD3  /home/sce/3.rna-seq/5.riboparser/04.density/SRR1944935_rna.txt RNA

#################################################
# merge all the RNA-seq files
rpf_Merge -l RNA.file.list -o RNA &>> RNA.log

cd ..
```

### 4.6 Calculate tri-nucleotide periodicity

1. Check the tri-nucleotide periodicity of RNA-seq
```bash
$ cd ./06.periodicity/

#################################################
# check the periodicity
rpf_Periodicity \
 -r ../05.merge/RNA_merged.txt \
 -m 30 --tis 0 --tts 0 -o RNA &>> RNA.log
```

### 4.7 `Meta-gene` analysis

Investigation of reads density in the vicinity of start and stop codons using `meta-gene` analysis.

1. `Meta-gene` analysis of RNA-seq

```bash
$ cd ./07.metaplot/

#################################################
# metagene analysis
rpf_Metaplot \
 -t ../../../1.reference/norm/gene.norm.txt \
 -r ../05.merge/RNA_merged.txt \
 -m 50 --mode bar -o RNA &>> RNA.log
```

### 4.8 Gene coverage

Examine the distribution of read density along the gene body.

1. Check gene density of RNA-seq

```bash
$ cd ./08.coverage/

#################################################
# check the reads density along with the gene body
rpf_Coverage \
 -t ../../../1.reference/norm/gene.norm.txt \
 -r ../05.merge/RNA_merged.txt \
 -m 50 --outlier \
 -b 10,100,10 \
 -n --heat \
 -o RNA &>> RNA.log

rpf_Percent \
 -t ../../../1.reference/norm/gene.norm.txt \
 -r ../05.merge/RNA_merged.txt \
 -n -m 50 \
 -f 0 \
 -o RNA &>> RNA.log
```

### 4.9 Check the repeatability of samples

1. Check the repeatability of RNA-seq

```bash
$ cd ./09.correlation/

#################################################
# calculate the samples replication of RNA-seq
rpf_Corr \
 -r ../05.merge/RNA_merged.txt \
 -o RNA &>> RNA.log
```


## 5. Perform Ribo-seq data analysis of `GSE67387` with `RiboParser`

### 5.0 Prepare the directory to store the results

```bash
$ cd
$ mkdir -p ./4.ribo-seq/5.riboparser/
$ cd ./4.ribo-seq/5.riboparser/
$ mkdir 01.qc 02.digestion 03.offset 04.density 05.merge \
 06.periodicity 07.metaplot 08.coverage 09.correlation 10.quantification \
 11.pausing_score 12.codon_occupancy 13.codon_decoding_time 14.codon_selection_time \
 15.coefficient_of_variation 16.meta_codon 17.shuffle 18.retrieve
```

### 5.1 Quality check of sequencing data

To ensure the reliability of downstream analyses, rigorous quality control (QC) of ribosome profiling data must include systematic evaluation of peak detection and gene detection rate.

`Peak Detection`:\
Calculate the dominant RPF length (expected range: 26–34 nt) and its proportion within the total reads (QC pass threshold: ≥70% in expected range).

`Gene Detection Rate`:\
Sampling from 5% to 95% of the data, the number of covered genes and their expression levels are calculated. When the slope of the gene count curve approaches 0, it typically indicates that the sequencing data is near saturation.

1. Explanation of `rpf_Check`

```bash
$ rpf_Check -h

Check the RPFs mapping condition.

Step1: Checking the input Arguments.

usage: rpf_Check [-h] -t TRANSCRIPT -b BAM -o OUTPUT [--thread THREAD] [-g {0,1}] [-a {star,hisat2,bowtie2}] [-r] [-l] [-s]

This script is used to summary the BAM condition.

options:
  -h, --help            show this help message and exit.
  --thread THREAD       the number of threads (default: 1). Suitable for large bam files > 1G.
                        It will take a lot of memory.
  -g {0,1}              filter the number of reads mapped loci (default: 0). [0]: all reads will be
                        used; [1]: reads with unique mapped loci will be used
  -a {star,hisat2,bowtie2}
                        screen the reads mapped loci from BAM file generate with different reads alignment methods (default: star).
  -r                    reads aligned to negative strand will also be counted. (default: False).
  -l                    only keep the longest transcripts (default: False).
  -s                    whether to calculate RPF saturation. (default: False). This step will take 
                        a lot of time and memory.

Required arguments:
  -t TRANSCRIPT         the input file name of gene annotation.
  -b BAM                the input file name of bam.
  -o OUTPUT             the prefix of output file.
```

2. Quality check of ribo-seq data

```bash
$ cd ./01.qc/

#################################################
# check the ribo-seq quality
for bam in ../../3.star/*Aligned.toTranscriptome.out.bam
do
prefix_name=$(basename $bam Aligned.toTranscriptome.out.bam)

rpf_Check -b $bam -s --thread 10 \
 -t ../../../1.reference/norm/gene.norm.txt \
 -o $prefix_name &>> $prefix_name".log"

done
```

3.  Results of `rpf_Check`

```bash
# results of sample SRR1944912
SRR1944912.bam # Filtered and sorted BAM file
SRR1944912.bam.bai # Index file for the BAM file
SRR1944912_gene_saturation.pdf # Barplot for gene saturation analysis
SRR1944912_gene_saturation.png # Barplot for gene saturation analysis
SRR1944912_gene_saturation.txt # Statistical results of gene saturation analysis
SRR1944912_length_distribution.pdf # Line graph of reads length distribution
SRR1944912_length_distribution.png # Line graph of reads length distribution
SRR1944912_length_distribution.txt # Statistical results of reads length distribution
SRR1944912.log # Log file of program execution
SRR1944912_reads_saturation.pdf # Boxplot for reads saturation analysis
SRR1944912_reads_saturation.png # Boxplot for reads saturation analysis
SRR1944912_reads_saturation.txt # Statistical results of reads saturation analysis
```


4. Explanation of `merge_length` and `merge_saturation`

```bash
Step1: Checking the input Arguments.
usage: merge_length [-h] -l LIST [LIST ...] -o OUTPUT

This script is used to merge the length distribution files

options:
  -h, --help            show this help message and exit

Required arguments:
  -l LIST [LIST ...], --list LIST [LIST ...]
                        List for length distribution files (e.g., '*length_distribution.txt').
  -o OUTPUT             prefix of output file name.
```

5. Integrating Ribo-seq quality check results for all samples

```bash
#################################################
# merge the ribo-seq quality results
merge_length -l *length_distribution.txt -o RIBO
merge_saturation -l *gene_saturation.txt -o RIBO

cd ..
```


### 5.2 Enzymatic bias in NGS library preparation

Ribonuclease digestion and ligation steps in Ribo-seq can introduce biases that affect the representation of RNA fragments.

`Ribonuclease Digestion Bias`:\
Ribonucleases may preferentially cleave certain RNA sequences or structures, leading to uneven fragmentation and skewed representation of RNA species in the data.

`Ligation Bias`:\
Ligation efficiency can vary by sequence context, RNA fragment length, or secondary structure.   This may result in overrepresentation of certain fragments due to more favorable ligation sites.

To reduce these biases, it's important to analyze digestion and ligation efficiency, optimize protocols, and include controls to ensure data accuracy and reliability.

1. Explanation of `rpf_Digest`

```bash
$ rpf_Digest -h

Step1: Checking the input Arguments.

usage: rpf_Digest [-h] -t TRANSCRIPT -s SEQUENCE -b BAM -o OUTPUT [-l] [--scale] [-m MIN] [-M MAX]

This script is used to Detect the digestion sites.

options:
  -h, --help     show this help message and exit
  -l             only retain the transcript with longest CDS of each gene (default: False).Recommended : True
  --scale        scale the motif matrix (default: False).
  -m MIN         the minimum reads length to keep (default: 20 nt).
  -M MAX         the maximum reads length to keep (default: 100 nt).

Required arguments:
  -t TRANSCRIPT  the name of input transcript file in TXT format.
  -s SEQUENCE    the name of input transcript sequence file in FA format.
  -b BAM         the name of mapping file in BAM format.
  -o OUTPUT      the name of output file. (prefix + _digestion_sites.txt)
```

2. Bias in restriction enzyme digestion and ligation in sequencing data

```bash
$ cd ./02.digestion/

#################################################
# check the reads digestion
for bam in ../01.qc/*.bam
do
prefix_name=$(basename $bam .bam)

rpf_Digest -b $bam -m 27 -M 33 --scale \
 -s ../../../1.reference/norm/gene.norm.rna.fa \
 -t ../../../1.reference/norm/gene.norm.txt \
 -o $prefix_name &>> $prefix_name".log"

done
```

3. Results of `rpf_Digest`

```bash
# results of SRR1944912
SRR1944912_3end_counts.txt # Statistical values of different bases at the 3'-end of reads
SRR1944912_3end_pwm.txt # PWM matrix of bases at the 3'-end of reads
SRR1944912_3end_seqlogo2.pdf # Seqlogo of bases at the 3'-end of reads
SRR1944912_5end_counts.txt # Statistical values of different bases at the 5'-end of reads
SRR1944912_5end_pwm.txt # PWM matrix of bases at the 5'-end of reads
SRR1944912_5end_seqlogo2.pdf # Seqlogo of bases at the 5'-end of reads
SRR1944912_digestion_sites.txt # Open reading frame aligned at the ends of reads
SRR1944912.log # Log file of program execution
SRR1944912_scaled_digestion_sites_plot.pdf # Heatmap of open reading frame aligned at the ends of reads
```


4. Explanation of `merge_digestion`

```bash
$ merge_digestion -h

Step1: Checking the input Arguments.
usage: merge_digestion [-h] -l LIST [LIST ...] -o OUTPUT

This script is used to merge the reads digestion files

options:
  -h, --help            show this help message and exit

Required arguments:
  -l LIST [LIST ...], --list LIST [LIST ...]
                        List for digestion files (e.g., '*_5end_pwm.txt').
  -o OUTPUT             prefix of output file name.
```

5. Integrating reads digestion results for all samples
```bash
#################################################
# merge the rpf digestion
merge_digestion -l *pwm.txt -o RIBO

cd ..
```


### 5.3 Use `RiboParser` to predict the optimal offset

The ability of Ribo-seq to analyze translation at codon resolution relies on accurately identifying the codons located at ribosomal `A-`, `P-`, and `E-`sites for each RPF.
The offset represents the distance of the 5’ end of the RPF to the first nucleotide of the `P-site` codon. 
Two commonly used methods are the ribosome structure-based model (`RSBM`) and the start/stop-based model (`SSCBM`).

1. Explanation of `rpf_Offset`

```bash
$ rpf_Offset -h

Step1: Checking the input Arguments.

usage: rpf_Offset [-h] -t TRANSCRIPT -b BAM -o OUTPUT [--mode {SSCBM,RSBM}] [-a {both,tis,tts}] [-l] [-m MIN] [-M MAX] [-p EXP_PEAK] [-s SHIFT] [--silence] [-d]

This script is used to detect the P-site offset.

options:
  -h, --help           show this help message and exit
  --mode {SSCBM,RSBM}  specify the mode of offset detect [SSCBM, RSBM]. (default: SSCBM).
  -a {both,tis,tts}    specify the alignment of reads for offset detect [both, tis, tts]. (default: both).
  -l                   only retain the transcript with longest CDS of each gene (default: False).
  -m MIN               the minimum reads length to keep (default: 27 nt).
  -M MAX               the maximum reads length to keep (default: 33 nt).
  -p EXP_PEAK          expected RPFs length fitted to ribosome structure [~30 nt] (default: 30 nt).
  -s SHIFT             psite shift for different RPFs length. (default: 2 nt).
  --silence            discard the warning information. (default: True).
  -d                   output the details of offset (default: False).

Required arguments:
  -t TRANSCRIPT        the name of input transcript filein TXT format.
  -b BAM               the name of mapping file in BAM format.
  -o OUTPUT            the prefix of output file. (prefix + _offset.txt)
```

2. Predict the optimal offset of Ribo-seq

```bash
$ cd ../03.offset/

#################################################
# predict the offset table
for bam in ../01.qc/*.bam
do
pre