Overview

RNA-Binding Site Seeker (rbsSeeker) is a unified software tool designed to identify significant RNA-binding protein (RBP) binding sites, including (peaks and individual cross-linking sites) from a variety of CLIP-seq datasets such as HITS-CLIP, PAR-CLIP, iCLIP and eCLIP. It is also compatible with miCLIP data analysis.

Leveraging Poisson and Hypergeometric modeling, rbsSeeker efficiently detects highly reliable peaks and cross-linking sites, including those characterized by deletions, truncations, and mutations.

• System requirements
• Run time
• Installation
• Input
• Output
• Basic usage
• Run rbsSeeker on testing dataset
• Acknowledgements
• Contact

System requirements

rbsSeeker has been tested on Linux systems with the following specifications:

• RAM: 64GB
• CPU: 16+ cores

Run time

Written in C and C++, rbsSeeker is highly efficient for analyzing CLIP-seq data. Most tasks are completed within 10 minutes during testing.

Installation

• Installing rbsSeeker on a Linux server is straightforward. Use the following commands:

  # Assume the installation directory is /username/software
  
  cd /username/software
  
  git clone https://github.com/kerenzhou062/rbsSeeker.git
  
  cd ./rbsSeeker
  
  make clean
  make
  
  export PATH=$PATH:/username/software/rbsSeeker/bin

Input

rbsSeeker requires the following input files:

• Genome sequence in FASTA format: Specify using the --fa argument. Chromosome names in this file must match those in the input bam file.

• FAI index file: Specify using the--fai argument. Generate this file using samtools:

  # # Example: Generate an FAI index for hg38.fa
  
  samtools faidx hg38.fa

• BAM file: Specify using the --bam argument. This file should contain sequence alignment data generated by tools such as STAR, HISAT2 and Tophat.

Note: rbsSeeker only accepts BAM files from single-end reads for peak/site calling.

Output

rbsSeeker may have 8 following output files in bed format (0-base) depends on the input arguments and your dataset.

Output file	Description
*_Peak.bed	identified peak-calling results
*_Mutation.bed	identified mutation sites (CIMS sites)
*_TC.bed or *_GA.bed	identified T-to-C mutation sites (when --cvs TC is set) or G-to-A mutation sites (when --cvs GA is set).
*_Truncation.bed	identified truncation sites (CITS)
*_Deletion.bed	identified deletion sites
*_Insertion.bed	identified insertion sites
*_End.bed	identified read termination sites

Here's the description of columns in the outputs:

Column name	Description
chrom	chromosome name
chromStart	Start genomic coordiate of the event (e.g. peak, deletion, truncation or mutation) (0-base)
chromEnd	End genomic coordinate of the event
name	Unique event ID
score	RPM of peak/site height
strand	Genomic sense (+) or antisense (-) strand of the event
extendSeq	Sequence extended ±10 bp from the site or peak center
motifPos	Start position of the motif in the extendSeq column (0-based). -1 indicates no motif found.
type	Type of peak or site
log10(p-value)	Log10-transformed p-value
log10(q-value)	Log10-transformed q-value
readNum	This column represents the number of reads that support a specific binding event.
height	This value represents the total number of reads covering a specific nucleotide position, regardless of whether they contain a variation event or not.
heightRpm	This is the height normalized to Reads Per Million (RPM) mapped reads.
mfold	This column represents the fold-enrichment of a binding signal compared to a background model.
ratio	this is the proportion of reads covering a site that also carry the specific variation.

Example

Below is an example demonstrating how to use rbsSeeker to identify N6-methyladenosine (m6A) sites at single-base resolution from miCLIP data.

• Preprocessing Ensure that raw miCLIP reads are properly processed (e.g., adapters trimmed, PCR duplicates removed, barcodes eliminated) and aligned to the appropriate genome (e.g., hg38). Use the resulting BAM file (miCLIP.sorted.bam) for downstream analysis.

• Running rbsSeeker

  # Example: Call m6A sites with -t 129600000 for human transcriptome
  # This step typically finishes within 10 minutes, depending on dataset size
  rbsSeeker -T CT -L 20 -t 129600000 -n 1 -H 3 -d 1 -p 0.05 -q 0.1 \
    -o ./output -P miCLIP --fa hg38.fa --fai hg38.fa.fai --treat miCLIP.sorted.bam > miCLIP.rbsSeeker.log

• Output files from rbsSeeker results
  

  • miCLIP_CT.bed
  • miCLIP_End.bed
  • miCLIP_Insertion.bed
  • miCLIP_Mutation.bed
  • miCLIP_Peak.bed
  • miCLIP_Truncation.bed

• Identify potential m6A sites (DRACH motif) (supposed the miCLIP was performed with Abcam antibody)

  • Substitutions at m6A site from miCLIP_CT.bed

    #substitutions at m6A site
    awk 'BEGIN{FS="\t";OFS="\t";}
    {
      if (FNR >1) {
        if ($8 == 8) {
          seq = substr($7, 9, 5);
          if (seq ~ /[AGT][AG]AC[ACT]/) {
            $4="mut|"seq"|"FNR;
            print $1, $2, $3, $4, $5, $6;
          }
        }
      }
    }' miCLIP_Mutation.bed > miCLIP.mut.bed

  • C->T mucations at +1 position of m6A site from miCLIP_Mutation.bed

    #C->T mucations at +1 position of m6A site
    awk 'BEGIN{FS="\t";OFS="\t";}
    {
      if (FNR >1) {
        if ($8 == 7) {
          if ($6 == "+") {
            $2 = $2 - 1;
            $3 = $2 + 1;
          }else{
            $3 = $3 + 1;
            $2 = $3 - 1;
          }
          seq = substr($7, 8, 5);
          if (seq ~ /[AGT][AG]AC[ACT]/) {
            $4="CT|"seq"|"FNR;
            print $1, $2, $3, $4, $5, $6;
          }
        }
      }
    }' miCLIP_CT.bed > miCLIP.CT.bed

  • Truncations at +2 position of m6A site from miCLIP_Truncation.bed

    #truncations at +2 position of m6A site
    awk 'BEGIN{FS="\t";OFS="\t";}
    {
      if (FNR >1) {
        if ($8 == 6) {
          if ($6 == "+") {
            $2 = $2 - 2;
            $3 = $2 + 1;
          }else{
            $3 = $3 + 2;
            $2 = $3 - 1;
          }
          seq = substr($7, 7, 5);
          if (seq ~ /[AGT][AG]AC[ACT]/) {
            $4="trunc|"seq"|"FNR;
            print $1, $2, $3, $4, $5, $6;
          }
        }
      }
    }' miCLIP_Truncation.bed > miCLIP.trunc.bed

• Pool identified sites together and get the final result miCLIP.merge.bed (bedtools is required)

  cat miCLIP.mut.bed miCLIP.CT.bed miCLIP.trunc.bed | awk 'BEGIN{FS="\t";OFS="\t";}
    {
      split($4,arr,"|");
      seq=arr[2];
      key=$1"\t"$2"\t"$3"\t"seq"\t"$6;
      readNumArr[key] += $5;
    }
    END{
      for (key in readNumArr) {
        split(key, arr,"\t");
        print arr[1], arr[2], arr[3], arr[4], readNumArr[key], arr[5];
      }
    }' | sort -k1,1 -k2,2n | bedtools intersect -a stdin \
    -b miCLIP.mut.bed miCLIP.CT.bed miCLIP.trunc.bed \
    -names mut CT trunc -s -wa -wb | \
    awk 'BEGIN{FS="\t";OFS="\t";}
    {
      key=$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6;
      if(key in hashArrA){
        hashArrA[key] += 1;
        hashArrB[key] = hashArrB[key]","$7;
      }else{
        hashArrA[key] = 1;
        hashArrB[key] = $7;
      }
    }
    END{
      for (key in hashArrA){
        print key, hashArrB[key], hashArrA[key];
      }
    }' | sort -k1,1 -k2,2n | \
    awk 'BEGIN{FS="\t";OFS="\t";}{$4=$4"|"FNR;print}'> miCLIP.merge.bed

Basic Usage

The available options of rbsSeeker are as follow:

Usage: rbsSeeker [options] --fa <genome file> --fai <genome fai> --treat <mapped alignments>
[options]
-v/--verbose               : verbose information
-V/--version               : rbsSeeker version
-h/--help                  : help informations
-R/--PCR                   : remove pcr duplictions[default is not removed]
-e/--rm                    : remove the muations in start or end sites[default is not removed]
-N/--norm                  : normalized the reads to locus[default is not normalized]
-k/--skip                  : skip reads spanning the intron[default is not skip]
-K/--rnafold               : RNAfold for your sequences[default is not]
--fa <string>              : genome file with FASTA format
--fai <string>             : genome fai file with FAI format
--treat <string>           : alignments treatment file with BAM format
--control <string>         : alignments control file with BAM format
-o/--outdir <string>       : output dir
-P/--prefix <string>       : prefix for output files
-t/--transcriptome <int>   : transcriptome size[e.g. in human, default=129600000]
-T/--cvs <string>          : conversion string[e.g. TC in PAR-CLIP, CT in miCLIP]
-c/--min-peak-len <int>    : minimum length for a peak [default>=10]
-C/--max-peak-len <int>    : maximum length for a peak [default<=1000000]
-i/--min-read-len <int>    : minimum read length [default>=18]
-a/--max-read-len <int>    : maximum read length [default<=1000000]
-n/--min-read-num <double> : minimum number of reads for calling a peak [Default=1]
-L/--max-locus-num <int>   : maximum locus number of reads for mapping to genome [Default=5]
-H/--min-height <double>   : minimum read height for calling a site [Default=5]
-r/--rpm <double>          : minimum rpm height for calling a site [Default=1]
-d/--min-var <double>      : minimum read number for calling a variational site [Default=1]
-p/--pval <double>         : minimum p value for calling a site [Default<=0.05]
-q/--qval <double>         : minimum q value for calling a site [Default<=0.05]
-M/--motif<string>         : search motif in the identified sequences  [default=NULL]
-s/--min-ratio<double>     : minimum ratio for variation [default>=0]
-S/--max-ratio<double>     : maximum ratio for variation [default<=1.0]
-m/--mfold<double>         : minimum fold-change for variation[default>=2]
-w/--window<int>           : window length for calculating the lamda of each peak [default=500]

Run rbsSeeker on testing dataset

Please refer to this guide for instructions on running rbsSeeker on a testing dataset.

Acknowledgements

We extend our gratitude to all contributors of public codes and libraries (e.g., BamTools) utilized by rbsSeeker.

Contact

• Jian-Hua Yang yangjh7@mail.sysu.edu.cn, RNA Information Center, School of Life Sciences, Sun Yat-Sen University
  
• Keren Zhou kzhou@stjude.org, Department of Pathology, St. Jude Children’s Research Hospital, Memphis, TN, USA