LncExpDB

1 RNA-seq data processing pipeline

This pipeline is used to analyze RNA sequencing data obtained from different biological contexts with a reference genome and annotation. It takes sample summary and runinfo files derived from NCBI SRAas inputs, performs quality control and alignment steps, and finally produces the gene/transcript reads count matrices.

Environment requirement

1. SRAtoolkit (ver. 2.9.6.1)
2. Trimmomatic (ver. 0.39)
3. HISAT2 (ver. 2.2.1)
4. STAR
5. RSeQC (ver. 2.6.4)
6. FeatureCounts
7. Kallisto
8. StringTie (ver. 2.0.6)

2 Data and Methods

2.1 Download sra files

The sra files are downloaded with the following command:

wget -c -t 0 download_url -P output_path

The -c parameter enables resume capability. -t 0 sets the retry attempts to infinite.-Pspecifies that the downloaded file will be saved in the output_path directory.

2.2 Convert sra files into fastq files

The sra files are converted into fastq files by fasterq-dump with the following command:

fasterq-dump -p -e 32 --split-3 sra_file -O output_dir

-p enables multi-threading to speed up the data extraction process. -e 32 specifies the use of 32 threads for processing the data.--split-3 splits the sra file into R1 and R2 fastq files for paired-end sequencing data or single fastq file for single-end sequencing data. -O specifies the output directory where the extracted fastq files will be saved.

2.3 Merge multiple fastq files

The runs from the same sample are merged with the following command:
For paired-end sequencing data

cat sample/*_1.fastq >sample/sample_1.fastq

cat sample/*_2.fastq >sample/sample_2.fastq

For single-end sequencing data

 cat sample/*.fastq >sample/sample.fastq

2.4 Infer library strandedness from fastq files

Library preparation protocols for RNA-Seq can be stranded or unstranded. The alignment and quantification tools depend on the strandedness to improve the specificity of alignment and reads allocation, separately. Manually curating this information from GEO or SRA is time-consuming. Fortunately, the strandedness can be auto inferred by infer_experiment.py script of RSeQC package with the following command:

python2 infer_experiment.py -i sample_Alignment-unsorted.sort.bam -r hg38BED > infer_experiment_log
            -i sample_Alignment-unsorted.sort.bam

specifies the input BAM file.
Here, this file is generated from randomly selecting a sample to align by HISAT2. -r provides a BED file containing genomic regions of all genes. infer_experiment_log contains the fraction of reads mapped to forward strand and reverse strand. This fraction is used to infer the strandedness.

2.5 Filter low-quality reads and sequencing adapters with Trimmomatic

Trimmomatic is used to filter out low-quality reads and sequencing adapters with the following command:
For paired-end sequencing data

java -jar Trimmomatic PE -threads threads -phred33 sample_1.fastq sample_2.fastq sample_paired_3.fastq sample_unpaired_3.fastq sample_paired_4.fastq sample_unpaired_4.fastq ILLUMINACLIP:%s:2:30:10:2:TRUE HEADCROP:3 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:15

For single-end sequencing data

java -jar Trimmomatic SE -threads threads -phred33 sample.fastq sample_paired.fastq ILLUMINACLIP:%s:2:30:10:2:TRUE HEADCROP:3 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:15

SE and PE indicate that the input is single-end and paired-end sequencing data, separately. -threads specifies the number of threads. -phred33 indicates that the quality scores in the fastq files are encoded using the Phred+33 format. sample_1.fastq and sample_2.fastq represent the forward and reverse reads of paired-end sequencing data. sample_paired_3.fastq and sample_paired_4.fastqare the paired reads that pass the filtering. sample_unpaired_3.fastq and sample_unpaired_4.fastqrecord the unpaired reads after filtering. ILLUMINACLIP cuts adapter and other illumina-specific sequences from the read.HEADCROP:3 indicates that the first 3 bases from the start of the read will be removed. LEADING:3 will remove low-quality bases from the start of the read, with a quality threshold of 3. Meanwhile, TRAILING:3 removes low-quality bases from the end of each sequencewith a quality threshold of 3.SLIDINGWINDOW:4:15 uses a sliding window approach to trim the sequence based on quality, with a window size of 4 bases. If the average quality within the window falls below 15, the sequence is trimmed from the end of that window.MINLEN:15 represents thatreads shorter than 15 bases after trimming will be discards.

2.6 Transcript quantification with Kallisto

Kallisto, a fast, accurate and alignment-free program, has been widely used for transcript quantification. Here, Kallisto is used to quantify transcripts with the following command:
For paired-end sequencing data

Kallisto quant -i kallisto_index -t threadsinferred_strandednesssample_paired_3.fastq sample_paired_4.fastq -o output_dir

For single-end sequencing data

Kallisto quant -i kallisto_index --single -l 150 -s 10 -t threadsinferred_strandednesssample_paired.fastq -o output_dir

-ispecifies the index file to be used by Kallisto, where kallisto_index is the pre-built transcriptome index. It can be constructed with the following command:

kallisto index -i kallisto_index lncbook.transcript.fa

lncbook.transcript.fa

contains all transcript sequences.

inferred_strandednessis

inferred from RSeQC.–singleindicates that the input sequencing data is single-end rather than paired-end. -l 150 specifies the average fragment length of the sequencing reads.-s 10 specifies the standard deviation of the fragment length. -t specifies the number of threads.

2.7 Align reads to genome with STAR

STAR is used to align reads to reference genome with the following command:
For paired-end sequencing data

STAR --runThreadN 4 --genomeDir star_reference --out FilterType BySJout --out FilterMultimapNmax 20 --out FilterMismatchNmax 999 --outFilterMismatchNoverLmax 0.04 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 --alignSJoverhangMin 8 --align SJDBoverhangMin 1 --sjdbScore 1 --out SAMtype BAM SortedByCoordinate --out BAMsortingThreadN 16 –readFilesIn sample_paired_3.fastq sample_paired_4.fastq --out FileNamePrefix sample

For single-end sequencing data

STAR --runThreadN 4 --genomeDir star_reference --out FilterTypeBySJout --out FilterMultimapNmax 20 --out FilterMismatchNmax 999 --outFilterMismatchNoverLmax 0.04 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 --alignSJoverhangMin 8 --align SJDBoverhangMin 1 --sjdbScore 1 --out SAMtype BAM SortedByCoordinate --out BAMsortingThreadN 16 --readFilesInsample_paired.fastq --out FileNamePrefixsample

This parameter combination was used by ENCODE project to quantify abundance of deep-sequencing samples. Specifically, --runThreadNspecifies the number of threads. --genomeDirspecifies the directory containing the STAR genome index files.--out FilterTypeBySJout specifies the default filtering mode where alignments are filtered based on splice junctions.--out FilterMultimapNmax 20 allows up to 20 multiple mappings. --out FilterMismatch Nmaxindicates the max allowed mismatches in the alignment. --out FilterMismatchNoverLmax0.04 limits the maximum number of mismatches to 4% of the read length. --align IntronMinspecifies the minimum intron length to be considered for alignment. --align IntronMaxspecifies the maximum intron length to be considered for alignment. --align MatesGapMaxspecifies the maximum gap size between paired-end reads. --align SJoverhangMinspecifies the minimum overhang length of the splice junctions to be considered for alignment. --sjdb Scoresets the score for each splice junction. A value of 1 is the minimum score threshold for splice junctions, used to filter out low-quality splice junctions and enhance the identification of high-quality ones.--out SAMtype BAM SortedByCoordinate specifies the output file format and sorting. BAM indicates that the output will be in BAM format, and SortedByCoordinate means that the alignments will be sorted by coordinate in the output file.--out BAMsorting ThreadN specifies the number of threads to be used for sorting the BAM file. –readFilesIn specifies the input fastq files. –outFileNamePrefix sets the prefix for the output files.

2.8 Gene quantification

In LncExpDB 1.0, gene expression was quantified using featureCounts.

featureCounts is used to quantify the reads on each gene with the following command:
For paired-end sequencing data

featureCounts -T threads -s inferred_strandedness -p -t exon -g gene_id -a lncbook_gtf -o sample_gene_quant.txt sample.Aligned.sortedByCoord.out.bam

For single-end sequencing data

featureCounts -T threads -s inferred_strandedness-t exon -g gene_id -a lncbook_gtf -o sample_gene_quant.txt sample.Aligned.sortedByCoord.out.bam -T

specifies the number of threads. -s defines the strand specificity of the reads. -p indicates that the input data consists of paired-end reads. -t exon specifies that the counting should be performed on exon features. -g gene_id specifies that the counts should be summarized by gene_id, using the gene identifier field from the annotation file.-a provides the annotation file. -o sets the output file name.

In LncExpDB 2.0, gene expression was quantified using RSEM for greater accuracy.

RSEM is used to quantify gene and isoform expression from RNA-seq data with the following command:
For paired-end sequencing data

frsem-calculate-expression --paired-end --bam --estimate-rspd --strandedness none -p 4 sampleAligned.sortedByCoord.out.bam rsem_reference/sample sample_RSEM

For single-end sequencing data

rsem-calculate-expression --single-cell-stranded --bam --estimate-rspd --strandedness none -p 4 sampleAligned.sortedByCoord.out.bam rsem_reference/sample sample_RSEM

specifies the number of threads. -s defines the strand specificity of the reads. -p indicates that the input data consists of paired-end reads. -t exon specifies that the counting should be performed on exon features. -g gene_id specifies that the counts should be summarized by gene_id, using the gene identifier field from the annotation file.-a provides the annotation file. -o sets the output file name.
-p sets the number of threads.
--strandedness sets read strand specificity (none/forward/reverse). --paired-end for paired-end reads; omit for single-end. --bam input is BAM instead of FASTQ.
--estimate-rspd estimates read start position distribution. rsem_reference/sample is the reference index. sample_RSEM sets the output file prefix.

3. Pipleline Usage

The pipeline can be performed with the following command:

python RNA-seq.py -a sra_result.csv -b sra_runinfo.csv -o output_path -t threads

The sra_result.csv and sra_runinfo.csv can be obtained from NCBI SRA. The sra_result.csv contains the basic sample information such as the accession, organism name, sequencing instrument, library strategy, file size and etc. The example that includes the complete columns is shown below:
The sra_runinfo.csv contains the basic information for each sequencing run. The basic information contains the total base number, average read length, file size, download link and etc. The example that includes the complete columns is shown below:
The output_path specifies the directory where all output files will be stored. The -t parameter specifies the number of threads used for Trimmomatic, Kallisto and featureCounts.

4. Output

The structure of the output folder is as follows:

The projectID of each dataset isusedas the name of the output directory. For each sample under this project, the read count matrices produced by featureCounts(sample_gene_quant.txt) and Kallisto (abundance.tsv) are stored in subfolders named with the sample IDs. The output logs are stored in the corresponding logs folder.

5. Abundance normalization pipeline

Normalization of RNA-Seq data has been proven essential to ensure accurate inferences and replication of biological findings. The major goal of normalization is to remove biology-unrelevant variances introducing by sequencing platform, library protocol, sequencing depth, gene/transcript length, or GC content. While comparing the expression difference from different datasets, it is necessary to remove the variances introducing by sequenicng platform and library protocol. Here, we only compare the expression difference between samples from the same dataset, ensuring the consistence of sequencing platform and library protocol. TMM (trimmed mean of M values), an effective between-sample normalization method to remove variances under the assumption that the most genes show little variation across different samples [A scaling normalization method for differential expression analysis of RNA-seq data], has been proven effective in identifying differentially expressed genes [Selecting between-sample RNA-Seq normalization methods from the perspective of their assumptions]. We perform gene/transcript abundance normalization with the calc NormFactors function in edgeR package. Then convert normalized abundance matrices into CPM, FPKM and TPM matrices.