macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

転写産物レベルで正確なリードカウントを行う RSEM

2021 1/9  タイトル修正

2021 1/15 コマンドと説明追記

2021 4/27 ベンチマーク論文追加

 

 RNA-Seqは転写産物の量を測定する方法に革命を起こしている。RNA-Seqデータからのトランスクリプト定量における重要な課題は、複数の遺伝子やアイソフォームにマップされたリードの取り扱いである。この問題は、配列決定されたゲノムがない場合のde novoトランスクリプトームアセンブリを用いた定量化において特に重要であり、どのトランスクリプトが同じ遺伝子のアイソフォームであるかを決定することは困難である。第二の重要な問題は、RNA-Seq実験のデザインであり、リード数、リードの長さ、リードがcDNA断片の片方または両方の末端から来るかどうかという点である。
 本研究では、シングルエンドまたはペアエンドのRNA-Seqデータから遺伝子とアイソフォームのアバンダンスを定量化するためのユーザーフレンドリーなソフトウェアパッケージであるRSEMを紹介する。RSEMは、アバンダンス推定値、95%信頼区間、可視化ファイルを出力し、RNA-Seqデータのシミュレーションも可能である。他の既存のツールとは対照的に、このソフトウェアはリファレンスゲノムを必要としない。したがって、de novoトランスクリプトームアセンブラと組み合わせることで、RSEMはゲノム配列のない種の正確なトランスクリプト定量を可能にする。シミュレーションおよび実データセットにおいて、RSEMはリファレンスゲノムに依存した定量法と比較して優れた性能を有している。また、RSEMが曖昧にマッピングされたリードを効果的に使用できることを利用して、ゲノムレベルの正確なアバンダンス推定値は、ショートシングルエンドリードを大量に使用した場合に最もよく得られることを示した。一方、単一遺伝子内のアイソフォームの相対的な頻度の推定は、各遺伝子の可能なスプライス形態の数に応じて、ペアエンドリードを使用することによって改善される可能性がある。
 RSEMは、RNA-Seqデータから転写産物の量を定量するための正確で使いやすいソフトウェアツールである。基準となるゲノムの存在に依存しないため、特にde novo転写産物アセンブリ定量化に有用である。また、現在では比較的高価なRNA-Seqを用いた定量実験をコスト効率よく設計するための貴重な指針となっている。

 

 

インストール

Github

#bioconda (link)
conda create -n rsem -y python=3.8
conda activate rsem
conda install -c bioconda rsem -y

> rsem-prepare-reference -h 

NAME

    rsem-prepare-reference - Prepare transcript references for RSEM and

    optionally build BOWTIE/BOWTIE2/STAR/HISAT2(transcriptome) indices.

 

SYNOPSIS

     rsem-prepare-reference [options] reference_fasta_file(s) reference_name

 

ARGUMENTS

    reference_fasta_file(s)

        Either a comma-separated list of Multi-FASTA formatted files OR a

        directory name. If a directory name is specified, RSEM will read all

        files with suffix ".fa" or ".fasta" in this directory. The files

        should contain either the sequences of transcripts or an entire

        genome, depending on whether the '--gtf' option is used.

 

    reference name

        The name of the reference used. RSEM will generate several

        reference-related files that are prefixed by this name. This name can

        contain path information (e.g. '/ref/mm9').

 

OPTIONS

    --gtf <file>

        If this option is on, RSEM assumes that 'reference_fasta_file(s)'

        contains the sequence of a genome, and will extract transcript

        reference sequences using the gene annotations specified in <file>,

        which should be in GTF format.

 

        If this and '--gff3' options are off, RSEM will assume

        'reference_fasta_file(s)' contains the reference transcripts. In this

        case, RSEM assumes that name of each sequence in the Multi-FASTA files

        is its transcript_id.

 

        (Default: off)

 

    --gff3 <file>

        The annotation file is in GFF3 format instead of GTF format. RSEM will

        first convert it to GTF format with the file name

        'reference_name.gtf'. Please make sure that 'reference_name.gtf' does

        not exist. (Default: off)

 

    --gff3-RNA-patterns <pattern>

        <pattern> is a comma-separated list of transcript categories, e.g.

        "mRNA,rRNA". Only transcripts that match the <pattern> will be

        extracted. (Default: "mRNA")

 

    --gff3-genes-as-transcripts

        This option is designed for untypical organisms, such as viruses,

        whose GFF3 files only contain genes. RSEM will assume each gene as a

        unique transcript when it converts the GFF3 file into GTF format.

 

    --trusted-sources <sources>

        <sources> is a comma-separated list of trusted sources, e.g.

        "ENSEMBL,HAVANA". Only transcripts coming from these sources will be

        extracted. If this option is off, all sources are accepted. (Default:

        off)

 

    --transcript-to-gene-map <file>

        Use information from <file> to map from transcript (isoform) ids to

        gene ids. Each line of <file> should be of the form:

 

        gene_id transcript_id

 

        with the two fields separated by a tab character.

 

        If you are using a GTF file for the "UCSC Genes" gene set from the

        UCSC Genome Browser, then the "knownIsoforms.txt" file (obtained from

        the "Downloads" section of the UCSC Genome Browser site) is of this

        format.

 

        If this option is off, then the mapping of isoforms to genes depends

        on whether the '--gtf' option is specified. If '--gtf' is specified,

        then RSEM uses the "gene_id" and "transcript_id" attributes in the GTF

        file. Otherwise, RSEM assumes that each sequence in the reference

        sequence files is a separate gene.

 

        (Default: off)

 

    --allele-to-gene-map <file>

        Use information from <file> to provide gene_id and transcript_id

        information for each allele-specific transcript. Each line of <file>

        should be of the form:

 

        gene_id transcript_id allele_id

 

        with the fields separated by a tab character.

 

        This option is designed for quantifying allele-specific expression. It

        is only valid if '--gtf' option is not specified. allele_id should be

        the sequence names presented in the Multi-FASTA-formatted files.

 

        (Default: off)

 

    --polyA

        Add poly(A) tails to the end of all reference isoforms. The length of

        poly(A) tail added is specified by '--polyA-length' option. STAR

        aligner users may not want to use this option. (Default: do not add

        poly(A) tail to any of the isoforms)

 

    --polyA-length <int>

        The length of the poly(A) tails to be added. (Default: 125)

 

    --no-polyA-subset <file>

        Only meaningful if '--polyA' is specified. Do not add poly(A) tails to

        those transcripts listed in <file>. <file> is a file containing a list

        of transcript_ids. (Default: off)

 

    --bowtie

        Build Bowtie indices. (Default: off)

 

    --bowtie-path <path>

        The path to the Bowtie executables. (Default: the path to Bowtie

        executables is assumed to be in the user's PATH environment variable)

 

    --bowtie2

        Build Bowtie 2 indices. (Default: off)

 

    --bowtie2-path <path>

        The path to the Bowtie 2 executables. (Default: the path to Bowtie 2

        executables is assumed to be in the user's PATH environment variable)

 

    --star

        Build STAR indices. (Default: off)

 

    --star-path <path>

        The path to STAR's executable. (Default: the path to STAR executable

        is assumed to be in user's PATH environment variable)

 

    --star-sjdboverhang <int>

        Length of the genomic sequence around annotated junction. It is only

        used for STAR to build splice junctions database and not needed for

        Bowtie or Bowtie2. It will be passed as the --sjdbOverhang option to

        STAR. According to STAR's manual, its ideal value is

        max(ReadLength)-1, e.g. for 2x101 paired-end reads, the ideal value is

        101-1=100. In most cases, the default value of 100 will work as well

        as the ideal value. (Default: 100)

 

    --hisat2-hca

        Build HISAT2 indices on the transcriptome according to Human Cell

        Atlas (HCA) SMART-Seq2 pipeline. (Default: off)

 

    --hisat2-path <path>

        The path to the HISAT2 executables. (Default: the path to HISAT2

        executables is assumed to be in the user's PATH environment variable)

 

    -p/--num-threads <int>

        Number of threads to use for building STAR's genome indices. (Default:

        1)

 

    -q/--quiet

        Suppress the output of logging information. (Default: off)

 

    -h/--help

        Show help information.

 

PRIOR-ENHANCED RSEM OPTIONS

    --prep-pRSEM

        A Boolean indicating whether to prepare reference files for pRSEM,

        including building Bowtie indices for a genome and selecting training

        set isoforms. The index files will be used for aligning ChIP-seq reads

        in prior-enhanced RSEM and the training set isoforms will be used for

        learning prior. A path to Bowtie executables and a mappability file in

        bigWig format are required when this option is on. Currently, Bowtie2

        is not supported for prior-enhanced RSEM. (Default: off)

 

    --mappability-bigwig-file <string>

        Full path to a whole-genome mappability file in bigWig format. This

        file is required for running prior-enhanced RSEM. It is used for

        selecting a training set of isoforms for prior-learning. This file can

        be either downloaded from UCSC Genome Browser or generated by GEM

        (Derrien et al., 2012, PLoS One). (Default: "")

 

DESCRIPTION

    This program extracts/preprocesses the reference sequences for RSEM and

    prior-enhanced RSEM. It can optionally build Bowtie indices (with

    '--bowtie' option) and/or Bowtie 2 indices (with '--bowtie2' option) using

    their default parameters. It can also optionally build STAR indices (with

    '--star' option) using parameters from ENCODE3's STAR-RSEM pipeline. For

    prior-enhanced RSEM, it can build Bowtie genomic indices and select

    training set isoforms (with options '--prep-pRSEM' and

    '--mappability-bigwig-file <string>'). If an alternative aligner is to be

    used, indices for that particular aligner can be built from either

    'reference_name.idx.fa' or 'reference_name.n2g.idx.fa' (see OUTPUT for

    details). This program is used in conjunction with the

    'rsem-calculate-expression' program.

 

OUTPUT

    This program will generate 'reference_name.grp', 'reference_name.ti',

    'reference_name.transcripts.fa', 'reference_name.seq',

    'reference_name.chrlist' (if '--gtf' is on), 'reference_name.idx.fa',

    'reference_name.n2g.idx.fa', optional Bowtie/Bowtie 2 index files, and

    optional STAR index files.

 

    'reference_name.grp', 'reference_name.ti', 'reference_name.seq', and

    'reference_name.chrlist' are used by RSEM internally.

 

    'reference_name.transcripts.fa' contains the extracted reference

    transcripts in Multi-FASTA format. Poly(A) tails are not added and it may

    contain lower case bases in its sequences if the corresponding genomic

    regions are soft-masked.

 

    'reference_name.idx.fa' and 'reference_name.n2g.idx.fa' are used by

    aligners to build their own indices. In these two files, all sequence

    bases are converted into upper case. In addition, poly(A) tails are added

    if '--polyA' option is set. The only difference between

    'reference_name.idx.fa' and 'reference_name.n2g.idx.fa' is that

    'reference_name.n2g.idx.fa' in addition converts all 'N' characters to 'G'

    characters. This conversion is in particular desired for aligners (e.g.

    Bowtie) that do not allow reads to overlap with 'N' characters in the

    reference sequences. Otherwise, 'reference_name.idx.fa' should be used to

    build the aligner's index files. RSEM uses 'reference_name.idx.fa' to

    build Bowtie 2 indices and 'reference_name.n2g.idx.fa' to build Bowtie

    indices. For visualizing the transcript-coordinate-based BAM files

    generated by RSEM in IGV, 'reference_name.idx.fa' should be imported as a

    "genome" (see Visualization section in README.md for details).

 

    If the whole genome is indexed for prior-enhanced RSEM, all the index

    files will be generated with prefix as 'reference_name_prsem'. Selected

    isoforms for training set are listed in the file

    'reference_name_prsem.training_tr_crd'

 

EXAMPLES

    1) Suppose we have mouse RNA-Seq data and want to use the UCSC mm9 version

    of the mouse genome. We have downloaded the UCSC Genes transcript

    annotations in GTF format (as mm9.gtf) using the Table Browser and the

    knownIsoforms.txt file for mm9 from the UCSC Downloads. We also have all

    chromosome files for mm9 in the directory '/data/mm9'. We want to put the

    generated reference files under '/ref' with name 'mouse_0'. We do not add

    any poly(A) tails. Please note that GTF files generated from UCSC's Table

    Browser do not contain isoform-gene relationship information. For the UCSC

    Genes annotation, this information can be obtained from the

    knownIsoforms.txt file. Suppose we want to build Bowtie indices and Bowtie

    executables are found in '/sw/bowtie'.

 

    There are two ways to write the command:

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --bowtie \

                            --bowtie-path /sw/bowtie \                  

                            /data/mm9/chr1.fa,/data/mm9/chr2.fa,...,/data/mm9/chrM.fa \

                            /ref/mouse_0

 

    OR

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --bowtie \

                            --bowtie-path /sw/bowtie \

                            /data/mm9 \

                            /ref/mouse_0

 

    2) Suppose we also want to build Bowtie 2 indices in the above example and

    Bowtie 2 executables are found in '/sw/bowtie2', the command will be:

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --bowtie \

                            --bowtie-path /sw/bowtie \

                            --bowtie2 \

                            --bowtie2-path /sw/bowtie2 \

                            /data/mm9 \

                            /ref/mouse_0

 

    3) Suppose we want to build STAR indices in the above example and save

    index files under '/ref' with name 'mouse_0'. Assuming STAR executable is

    '/sw/STAR', the command will be:

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --star \

                            --star-path /sw/STAR \

                            -p 8 \

                            /data/mm9/chr1.fa,/data/mm9/chr2.fa,...,/data/mm9/chrM.fa \

                            /ref/mouse_0

 

    OR

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --star \

                            --star-path /sw/STAR \

                            -p 8 \

                            /data/mm9

                            /ref/mouse_0

 

    STAR genome index files will be saved under '/ref/'.

 

    4) Suppose we want to prepare references for prior-enhanced RSEM in the

    above example. In this scenario, both STAR and Bowtie are required to

    build genomic indices - STAR for RNA-seq reads and Bowtie for ChIP-seq

    reads. Assuming their executables are under '/sw/STAR' and '/sw/Bowtie',

    respectively. Also, assuming the mappability file for mouse genome is

    '/data/mm9.bigWig'. The command will be:

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --star \

                            --star-path /sw/STAR \

                            -p 8 \

                            --prep-pRSEM \

                            --bowtie-path /sw/Bowtie \

                            --mappability-bigwig-file /data/mm9.bigWig \

                            /data/mm9/chr1.fa,/data/mm9/chr2.fa,...,/data/mm9/chrM.fa \

                            /ref/mouse_0

 

    OR

 

     rsem-prepare-reference --gtf mm9.gtf \

                            --transcript-to-gene-map knownIsoforms.txt \

                            --star \

                            --star-path /sw/STAR \

                            -p 8 \

                            --prep-pRSEM \

                            --bowtie-path /sw/Bowtie \

                            --mappability-bigwig-file /data/mm9.bigWig \

                            /data/mm9

                            /ref/mouse_0

 

    Both STAR and Bowtie's index files will be saved under '/ref/'. Bowtie

    files will have name prefix 'mouse_0_prsem'

 

    5) Suppose we only have transcripts from EST tags stored in 'mm9.fasta'

    and isoform-gene information stored in 'mapping.txt'. We want to add 125bp

    long poly(A) tails to all transcripts. The reference_name is set as

    'mouse_125'. In addition, we do not want to build Bowtie/Bowtie 2 indices,

    and will use an alternative aligner to align reads against either

    'mouse_125.idx.fa' or 'mouse_125.idx.n2g.fa':

 

     rsem-prepare-reference --transcript-to-gene-map mapping.txt \

                            --polyA

                            mm9.fasta \

                            mouse_125

 

rsem-calculate-expression -h

NAME

    rsem-calculate-expression - Estimate gene and isoform expression from

    RNA-Seq data.

 

SYNOPSIS

     rsem-calculate-expression [options] upstream_read_file(s) reference_name sample_name 

     rsem-calculate-expression [options] --paired-end upstream_read_file(s) downstream_read_file(s) reference_name sample_name 

     rsem-calculate-expression [options] --alignments [--paired-end] input reference_name sample_name

 

ARGUMENTS

    upstream_read_files(s)

        Comma-separated list of files containing single-end reads or upstream

        reads for paired-end data. By default, these files are assumed to be

        in FASTQ format. If the --no-qualities option is specified, then FASTA

        format is expected.

 

    downstream_read_file(s)

        Comma-separated list of files containing downstream reads which are

        paired with the upstream reads. By default, these files are assumed to

        be in FASTQ format. If the --no-qualities option is specified, then

        FASTA format is expected.

 

    input

        SAM/BAM/CRAM formatted input file. If "-" is specified for the

        filename, the input is instead assumed to come from standard input.

        RSEM requires all alignments of the same read group together. For

        paired-end reads, RSEM also requires the two mates of any alignment be

        adjacent. In addition, RSEM does not allow the SEQ and QUAL fields to

        be empty. See Description section for how to make input file obey

        RSEM's requirements.

 

    reference_name

        The name of the reference used. The user must have run

        'rsem-prepare-reference' with this reference_name before running this

        program.

 

    sample_name

        The name of the sample analyzed. All output files are prefixed by this

        name (e.g., sample_name.genes.results)

 

BASIC OPTIONS

    --paired-end

        Input reads are paired-end reads. (Default: off)

 

    --no-qualities

        Input reads do not contain quality scores. (Default: off)

 

    --strandedness <none|forward|reverse>

        This option defines the strandedness of the RNA-Seq reads. It

        recognizes three values: 'none', 'forward', and 'reverse'. 'none'

        refers to non-strand-specific protocols. 'forward' means all

        (upstream) reads are derived from the forward strand. 'reverse' means

        all (upstream) reads are derived from the reverse strand. If

        'forward'/'reverse' is set, the '--norc'/'--nofw' Bowtie/Bowtie 2

        option will also be enabled to avoid aligning reads to the opposite

        strand. For Illumina TruSeq Stranded protocols, please use 'reverse'.

        (Default: 'none')

 

    -p/--num-threads <int>

        Number of threads to use. Both Bowtie/Bowtie2, expression estimation

        and 'samtools sort' will use this many threads. (Default: 1)

 

    --alignments

        Input file contains alignments in SAM/BAM/CRAM format. The exact file

        format will be determined automatically. (Default: off)

 

    --fai <file>

        If the header section of input alignment file does not contain

        reference sequence information, this option should be turned on.

        <file> is a FAI format file containing each reference sequence's name

        and length. Please refer to the SAM official website for the details

        of FAI format. (Default: off)

 

    --bowtie2

        Use Bowtie 2 instead of Bowtie to align reads. Since currently RSEM

        does not handle indel, local and discordant alignments, the Bowtie2

        parameters are set in a way to avoid those alignments. In particular,

        we use options '--sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1

        --score-min L,0,-0.1' by default. The last parameter of '--score-min',

        '-0.1', is the negative of maximum mismatch rate. This rate can be set

        by option '--bowtie2-mismatch-rate'. If reads are paired-end, we

        additionally use options '--no-mixed' and '--no-discordant'. (Default:

        off)

 

    --star

        Use STAR to align reads. Alignment parameters are from ENCODE3's

        STAR-RSEM pipeline. To save computational time and memory resources,

        STAR's Output BAM file is unsorted. It is stored in RSEM's temporary

        directory with name as 'sample_name.bam'. Each STAR job will have its

        own private copy of the genome in memory. (Default: off)

 

    --hisat2-hca

        Use HISAT2 to align reads to the transcriptome according to Human Cell

        Atlast SMART-Seq2 pipeline. In particular, we use HISAT parameters "-k

        10 --secondary --rg-id=$sampleToken --rg SM:$sampleToken --rg

        LB:$sampleToken --rg PL:ILLUMINA --rg PU:$sampleToken --new-summary

        --summary-file $sampleName.log --met-file $sampleName.hisat2.met.txt

        --met 5 --mp 1,1 --np 1 --score-min L,0,-0.1 --rdg 99999999,99999999

        --rfg 99999999,99999999 --no-spliced-alignment --no-softclip --seed

        12345". If inputs are paired-end reads, we additionally use parameters

        "--no-mixed --no-discordant". (Default: off)

 

    --append-names

        If gene_name/transcript_name is available, append it to the end of

        gene_id/transcript_id (separated by '_') in files

        'sample_name.isoforms.results' and 'sample_name.genes.results'.

        (Default: off)

 

    --seed <uint32>

        Set the seed for the random number generators used in calculating

        posterior mean estimates and credibility intervals. The seed must be a

        non-negative 32 bit integer. (Default: off)

 

    --single-cell-prior

        By default, RSEM uses Dirichlet(1) as the prior to calculate posterior

        mean estimates and credibility intervals. However, much less genes are

        expressed in single cell RNA-Seq data. Thus, if you want to compute

        posterior mean estimates and/or credibility intervals and you have

        single-cell RNA-Seq data, you are recommended to turn on this option.

        Then RSEM will use Dirichlet(0.1) as the prior which encourage the

        sparsity of the expression levels. (Default: off)

 

    --calc-pme

        Run RSEM's collapsed Gibbs sampler to calculate posterior mean

        estimates. (Default: off)

 

    --calc-ci

        Calculate 95% credibility intervals and posterior mean estimates. The

        credibility level can be changed by setting '--ci-credibility-level'.

        (Default: off)

 

    -q/--quiet

        Suppress the output of logging information. (Default: off)

 

    -h/--help

        Show help information.

 

    --version

        Show version information.

 

OUTPUT OPTIONS

    --sort-bam-by-read-name

        Sort BAM file aligned under transcript coordidate by read name.

        Setting this option on will produce deterministic maximum likelihood

        estimations from independent runs. Note that sorting will take long

        time and lots of memory. (Default: off)

 

    --no-bam-output

        Do not output any BAM file. (Default: off)

 

    --sampling-for-bam

        When RSEM generates a BAM file, instead of outputting all alignments a

        read has with their posterior probabilities, one alignment is sampled

        according to the posterior probabilities. The sampling procedure

        includes the alignment to the "noise" transcript, which does not

        appear in the BAM file. Only the sampled alignment has a weight of 1.

        All other alignments have weight 0. If the "noise" transcript is

        sampled, all alignments appeared in the BAM file should have weight 0.

        (Default: off)

 

    --output-genome-bam

        Generate a BAM file, 'sample_name.genome.bam', with alignments mapped

        to genomic coordinates and annotated with their posterior

        probabilities. In addition, RSEM will call samtools (included in RSEM

        package) to sort and index the bam file.

        'sample_name.genome.sorted.bam' and

        'sample_name.genome.sorted.bam.bai' will be generated. (Default: off)

 

    --sort-bam-by-coordinate

        Sort RSEM generated transcript and genome BAM files by coordinates and

        build associated indices. (Default: off)

 

    --sort-bam-memory-per-thread <string>

        Set the maximum memory per thread that can be used by 'samtools sort'.

        <string> represents the memory and accepts suffices 'K/M/G'. RSEM will

        pass <string> to the '-m' option of 'samtools sort'. Note that the

        default used here is different from the default used by samtools.

        (Default: 1G)

 

ALIGNER OPTIONS

    --seed-length <int>

        Seed length used by the read aligner. Providing the correct value is

        important for RSEM. If RSEM runs Bowtie, it uses this value for

        Bowtie's seed length parameter. Any read with its or at least one of

        its mates' (for paired-end reads) length less than this value will be

        ignored. If the references are not added poly(A) tails, the minimum

        allowed value is 5, otherwise, the minimum allowed value is 25. Note

        that this script will only check if the value >= 5 and give a warning

        message if the value < 25 but >= 5. (Default: 25)

 

    --phred33-quals

        Input quality scores are encoded as Phred+33. This option is used by

        Bowtie, Bowtie 2 and HISAT2. (Default: on)

 

    --phred64-quals

        Input quality scores are encoded as Phred+64 (default for GA Pipeline

        ver. >= 1.3). This option is used by Bowtie, Bowtie 2 and HISAT2.

        (Default: off)

 

    --solexa-quals

        Input quality scores are solexa encoded (from GA Pipeline ver. < 1.3).

        This option is used by Bowtie, Bowtie 2 and HISAT2. (Default: off)

 

    --bowtie-path <path>

        The path to the Bowtie executables. (Default: the path to the Bowtie

        executables is assumed to be in the user's PATH environment variable)

 

    --bowtie-n <int>

        (Bowtie parameter) max # of mismatches in the seed. (Range: 0-3,

        Default: 2)

 

    --bowtie-e <int>

        (Bowtie parameter) max sum of mismatch quality scores across the

        alignment. (Default: 99999999)

 

    --bowtie-m <int>

        (Bowtie parameter) suppress all alignments for a read if > <int> valid

        alignments exist. (Default: 200)

 

    --bowtie-chunkmbs <int>

        (Bowtie parameter) memory allocated for best first alignment

        calculation (Default: 0 - use Bowtie's default)

 

    --bowtie2-path <path>

        (Bowtie 2 parameter) The path to the Bowtie 2 executables. (Default:

        the path to the Bowtie 2 executables is assumed to be in the user's

        PATH environment variable)

 

    --bowtie2-mismatch-rate <double>

        (Bowtie 2 parameter) The maximum mismatch rate allowed. (Default: 0.1)

 

    --bowtie2-k <int>

        (Bowtie 2 parameter) Find up to <int> alignments per read. (Default:

        200)

 

    --bowtie2-sensitivity-level <string>

        (Bowtie 2 parameter) Set Bowtie 2's preset options in --end-to-end

        mode. This option controls how hard Bowtie 2 tries to find alignments.

        <string> must be one of "very_fast", "fast", "sensitive" and

        "very_sensitive". The four candidates correspond to Bowtie 2's

        "--very-fast", "--fast", "--sensitive" and "--very-sensitive" options.

        (Default: "sensitive" - use Bowtie 2's default)

 

    --star-path <path>

        The path to STAR's executable. (Default: the path to STAR executable

        is assumed to be in user's PATH environment variable)

 

    --star-gzipped-read-file

        (STAR parameter) Input read file(s) is compressed by gzip. (Default:

        off)

 

    --star-bzipped-read-file

        (STAR parameter) Input read file(s) is compressed by bzip2. (Default:

        off)

 

    --star-output-genome-bam

        (STAR parameter) Save the BAM file from STAR alignment under genomic

        coordinate to 'sample_name.STAR.genome.bam'. This file is NOT sorted

        by genomic coordinate. In this file, according to STAR's manual,

        'paired ends of an alignment are always adjacent, and multiple

        alignments of a read are adjacent as well'. (Default: off)

 

    --hisat2-path <path>

        The path to HISAT2's executable. (Default: the path to HISAT2

        executable is assumed to be in user's PATH environment variable)

 

ADVANCED OPTIONS

    --tag <string>

        The name of the optional field used in the SAM input for identifying a

        read with too many valid alignments. The field should have the format

        <tagName>:i:<value>, where a <value> bigger than 0 indicates a read

        with too many alignments. (Default: "")

 

    --fragment-length-min <int>

        Minimum read/insert length allowed. This is also the value for the

        Bowtie/Bowtie2 -I option. (Default: 1)

 

    --fragment-length-max <int>

        Maximum read/insert length allowed. This is also the value for the

        Bowtie/Bowtie 2 -X option. (Default: 1000)

 

    --fragment-length-mean <double>

        (single-end data only) The mean of the fragment length distribution,

        which is assumed to be a Gaussian. (Default: -1, which disables use of

        the fragment length distribution)

 

    --fragment-length-sd <double>

        (single-end data only) The standard deviation of the fragment length

        distribution, which is assumed to be a Gaussian. (Default: 0, which

        assumes that all fragments are of the same length, given by the

        rounded value of --fragment-length-mean)

 

    --estimate-rspd

        Set this option if you want to estimate the read start position

        distribution (RSPD) from data. Otherwise, RSEM will use a uniform

        RSPD. (Default: off)

 

    --num-rspd-bins <int>

        Number of bins in the RSPD. Only relevant when '--estimate-rspd' is

        specified. Use of the default setting is recommended. (Default: 20)

 

    --gibbs-burnin <int>

        The number of burn-in rounds for RSEM's Gibbs sampler. Each round

        passes over the entire data set once. If RSEM can use multiple

        threads, multiple Gibbs samplers will start at the same time and all

        samplers share the same burn-in number. (Default: 200)

 

    --gibbs-number-of-samples <int>

        The total number of count vectors RSEM will collect from its Gibbs

        samplers. (Default: 1000)

 

    --gibbs-sampling-gap <int>

        The number of rounds between two succinct count vectors RSEM collects.

        If the count vector after round N is collected, the count vector after

        round N + <int> will also be collected. (Default: 1)

 

    --ci-credibility-level <double>

        The credibility level for credibility intervals. (Default: 0.95)

 

    --ci-memory <int>

        Maximum size (in memory, MB) of the auxiliary buffer used for

        computing credibility intervals (CI). (Default: 1024)

 

    --ci-number-of-samples-per-count-vector <int>

        The number of read generating probability vectors sampled per sampled

        count vector. The crebility intervals are calculated by first sampling

        P(C | D) and then sampling P(Theta | C) for each sampled count vector.

        This option controls how many Theta vectors are sampled per sampled

        count vector. (Default: 50)

 

    --keep-intermediate-files

        Keep temporary files generated by RSEM. RSEM creates a temporary

        directory, 'sample_name.temp', into which it puts all intermediate

        output files. If this directory already exists, RSEM overwrites all

        files generated by previous RSEM runs inside of it. By default, after

        RSEM finishes, the temporary directory is deleted. Set this option to

        prevent the deletion of this directory and the intermediate files

        inside of it. (Default: off)

 

    --temporary-folder <string>

        Set where to put the temporary files generated by RSEM. If the folder

        specified does not exist, RSEM will try to create it. (Default:

        sample_name.temp)

 

    --time

        Output time consumed by each step of RSEM to 'sample_name.time'.

        (Default: off)

 

PRIOR-ENHANCED RSEM OPTIONS

    --run-pRSEM

        Running prior-enhanced RSEM (pRSEM). Prior parameters, i.e. isoform's

        initial pseudo-count for RSEM's Gibbs sampling, will be learned from

        input RNA-seq data and an external data set. When pRSEM needs and only

        needs ChIP-seq peak information to partition isoforms (e.g. in pRSEM's

        default partition model), either ChIP-seq peak file (with the

        '--chipseq-peak-file' option) or ChIP-seq FASTQ files for target and

        input and the path for Bowtie executables are required (with the

        '--chipseq-target-read-files <string>', '--chipseq-control-read-files

        <string>', and '--bowtie-path <path> options), otherwise, ChIP-seq

        FASTQ files for target and control and the path to Bowtie executables

        are required. (Default: off)

 

    --chipseq-peak-file <string>

        Full path to a ChIP-seq peak file in ENCODE's narrowPeak, i.e. BED6+4,

        format. This file is used when running prior-enhanced RSEM in the

        default two-partition model. It partitions isoforms by whether they

        have ChIP-seq overlapping with their transcription start site region

        or not. Each partition will have its own prior parameter learned from

        a training set. This file can be either gzipped or ungzipped.

        (Default: "")

 

    --chipseq-target-read-files <string>

        Comma-separated full path of FASTQ read file(s) for ChIP-seq target.

        This option is used when running prior-enhanced RSEM. It provides

        information to calculate ChIP-seq peaks and signals. The file(s) can

        be either ungzipped or gzipped with a suffix '.gz' or '.gzip'. The

        options '--bowtie-path <path>' and '--chipseq-control-read-files

        <string>' must be defined when this option is specified. (Default: "")

 

    --chipseq-control-read-files <string>

        Comma-separated full path of FASTQ read file(s) for ChIP-seq conrol.

        This option is used when running prior-enhanced RSEM. It provides

        information to call ChIP-seq peaks. The file(s) can be either

        ungzipped or gzipped with a suffix '.gz' or '.gzip'. The options

        '--bowtie-path <path>' and '--chipseq-target-read-files <string>' must

        be defined when this option is specified. (Default: "")

 

    --chipseq-read-files-multi-targets <string>

        Comma-separated full path of FASTQ read files for multiple ChIP-seq

        targets. This option is used when running prior-enhanced RSEM, where

        prior is learned from multiple complementary data sets. It provides

        information to calculate ChIP-seq signals. All files can be either

        ungzipped or gzipped with a suffix '.gz' or '.gzip'. When this option

        is specified, the option '--bowtie-path <path>' must be defined and

        the option '--partition-model <string>' will be set to 'cmb_lgt'

        automatically. (Default: "")

 

    --chipseq-bed-files-multi-targets <string>

        Comma-separated full path of BED files for multiple ChIP-seq targets.

        This option is used when running prior-enhanced RSEM, where prior is

        learned from multiple complementary data sets. It provides information

        of ChIP-seq signals and must have at least the first six BED columns.

        All files can be either ungzipped or gzipped with a suffix '.gz' or

        '.gzip'. When this option is specified, the option '--partition-model

        <string>' will be set to 'cmb_lgt' automatically. (Default: "")

 

    --cap-stacked-chipseq-reads

        Keep a maximum number of ChIP-seq reads that aligned to the same

        genomic interval. This option is used when running prior-enhanced

        RSEM, where prior is learned from multiple complementary data sets.

        This option is only in use when either

        '--chipseq-read-files-multi-targets <string>' or

        '--chipseq-bed-files-multi-targets <string>' is specified. (Default:

        off)

 

    --n-max-stacked-chipseq-reads <int>

        The maximum number of stacked ChIP-seq reads to keep. This option is

        used when running prior-enhanced RSEM, where prior is learned from

        multiple complementary data sets. This option is only in use when the

        option '--cap-stacked-chipseq-reads' is set. (Default: 5)

 

    --partition-model <string>

        A keyword to specify the partition model used by prior-enhanced RSEM.

        It must be one of the following keywords:

 

        - pk

          Partitioned by whether an isoform has a ChIP-seq peak overlapping

          with its transcription start site (TSS) region. The TSS region is

          defined as [TSS-500bp, TSS+500bp]. For simplicity, we refer this

          type of peak as 'TSS peak' when explaining other keywords.

 

        - pk_lgtnopk

          First partitioned by TSS peak. Then, for isoforms in the 'no TSS

          peak' set, a logistic model is employed to further classify them

          into two partitions.

 

        - lm3, lm4, lm5, or lm6

          Based on their ChIP-seq signals, isoforms are classified into 3, 4,

          5, or 6 partitions by a linear regression model.

 

        - nopk_lm2pk, nopk_lm3pk, nopk_lm4pk, or

        nopk_lm5pk

          First partitioned by TSS peak. Then, for isoforms in the 'with TSS

          peak' set, a linear regression model is employed to further classify

          them into 2, 3, 4, or 5 partitions.

 

        - pk_lm2nopk, pk_lm3nopk, pk_lm4nopk, or

        pk_lm5nopk

          First partitioned by TSS peak. Then, for isoforms in the 'no TSS

          peak' set, a linear regression model is employed to further classify

          them into 2, 3, 4, or 5 partitions.

 

        - cmb_lgt

          Using a logistic regression to combine TSS signals from multiple

          complementary data sets and partition training set isoform into

          'expressed' and 'not expressed'. This partition model is only in use

          when either '--chipseq-read-files-multi-targets <string>' or

          '--chipseq-bed-files-multi-targets <string> is specified.

 

        Parameters for all the above models are learned from a training set.

        For detailed explanations, please see prior-enhanced RSEM's paper.

        (Default: 'pk')

 

DEPRECATED OPTIONS

        The options in this section are deprecated. They are here only for

        compatibility reasons and may be removed in future releases.

 

    --sam

        Inputs are alignments in SAM format. (Default: off)

 

    --bam

        Inputs are alignments in BAM format. (Default: off)

 

    --strand-specific

        Equivalent to '--strandedness forward'. (Default: off)

 

    --forward-prob <double>

        Probability of generating a read from the forward strand of a

        transcript. Set to 1 for a strand-specific protocol where all

        (upstream) reads are derived from the forward strand, 0 for a

        strand-specific protocol where all (upstream) read are derived from

        the reverse strand, or 0.5 for a non-strand-specific protocol.

        (Default: off)

 

DESCRIPTION

    In its default mode, this program aligns input reads against a reference

    transcriptome with Bowtie and calculates expression values using the

    alignments. RSEM assumes the data are single-end reads with quality

    scores, unless the '--paired-end' or '--no-qualities' options are

    specified. Alternatively, users can use STAR to align reads using the

    '--star' option. RSEM has provided options in 'rsem-prepare-reference' to

    prepare STAR's genome indices. Users may use an alternative aligner by

    specifying '--alignments', and providing an alignment file in SAM/BAM/CRAM

    format. However, users should make sure that they align against the

    indices generated by 'rsem-prepare-reference' and the alignment file

    satisfies the requirements mentioned in ARGUMENTS section.

 

    One simple way to make the alignment file satisfying RSEM's requirements

    is to use the 'convert-sam-for-rsem' script. This script accepts

    SAM/BAM/CRAM files as input and outputs a BAM file. For example, type the

    following command to convert a SAM file, 'input.sam', to a ready-for-use

    BAM file, 'input_for_rsem.bam':

 

      convert-sam-for-rsem input.sam input_for_rsem

 

    For details, please refer to 'convert-sam-for-rsem's documentation page.

 

NOTES

    1. Users must run 'rsem-prepare-reference' with the appropriate reference

    before using this program.

 

    2. For single-end data, it is strongly recommended that the user provide

    the fragment length distribution parameters (--fragment-length-mean and

    --fragment-length-sd). For paired-end data, RSEM will automatically learn

    a fragment length distribution from the data.

 

    3. Some aligner parameters have default values different from their

    original settings.

 

    4. With the '--calc-pme' option, posterior mean estimates will be

    calculated in addition to maximum likelihood estimates.

 

    5. With the '--calc-ci' option, 95% credibility intervals and posterior

    mean estimates will be calculated in addition to maximum likelihood

    estimates.

 

    6. The temporary directory and all intermediate files will be removed when

    RSEM finishes unless '--keep-intermediate-files' is specified.

 

    With the '--run-pRSEM' option and associated options (see section

    'PRIOR-ENHANCED RSEM OPTIONS' above for details), prior-enhanced RSEM will

    be running. Prior parameters will be learned from supplied external data

    set(s) and assigned as initial pseudo-counts for isoforms in the

    corresponding partition for Gibbs sampling.

 

OUTPUT

    sample_name.isoforms.results

        File containing isoform level expression estimates. The first line

        contains column names separated by the tab character. The format of

        each line in the rest of this file is:

 

        transcript_id gene_id length effective_length expected_count TPM FPKM

        IsoPct [posterior_mean_count posterior_standard_deviation_of_count

        pme_TPM pme_FPKM IsoPct_from_pme_TPM TPM_ci_lower_bound

        TPM_ci_upper_bound TPM_coefficient_of_quartile_variation

        FPKM_ci_lower_bound FPKM_ci_upper_bound

        FPKM_coefficient_of_quartile_variation]

 

        Fields are separated by the tab character. Fields within "" are

        optional. They will not be presented if neither '--calc-pme' nor

        '--calc-ci' is set.

 

        'transcript_id' is the transcript name of this transcript. 'gene_id'

        is the gene name of the gene which this transcript belongs to (denote

        this gene as its parent gene). If no gene information is provided,

        'gene_id' and 'transcript_id' are the same.

 

        'length' is this transcript's sequence length (poly(A) tail is not

        counted). 'effective_length' counts only the positions that can

        generate a valid fragment. If no poly(A) tail is added,

        'effective_length' is equal to transcript length - mean fragment

        length + 1. If one transcript's effective length is less than 1, this

        transcript's both effective length and abundance estimates are set to

        0.

 

        'expected_count' is the sum of the posterior probability of each read

        comes from this transcript over all reads. Because 1) each read

        aligning to this transcript has a probability of being generated from

        background noise; 2) RSEM may filter some alignable low quality reads,

        the sum of expected counts for all transcript are generally less than

        the total number of reads aligned.

 

        'TPM' stands for Transcripts Per Million. It is a relative measure of

        transcript abundance. The sum of all transcripts' TPM is 1 million.

        'FPKM' stands for Fragments Per Kilobase of transcript per Million

        mapped reads. It is another relative measure of transcript abundance.

        If we define l_bar be the mean transcript length in a sample, which

        can be calculated as

 

        l_bar = \sum_i TPM_i / 10^6 * effective_length_i (i goes through every

        transcript),

 

        the following equation is hold:

 

        FPKM_i = 10^3 / l_bar * TPM_i.

 

        We can see that the sum of FPKM is not a constant across samples.

 

        'IsoPct' stands for isoform percentage. It is the percentage of this

        transcript's abandunce over its parent gene's abandunce. If its parent

        gene has only one isoform or the gene information is not provided,

        this field will be set to 100.

 

        'posterior_mean_count', 'pme_TPM', 'pme_FPKM' are posterior mean

        estimates calculated by RSEM's Gibbs sampler.

        'posterior_standard_deviation_of_count' is the posterior standard

        deviation of counts. 'IsoPct_from_pme_TPM' is the isoform percentage

        calculated from 'pme_TPM' values.

 

        'TPM_ci_lower_bound', 'TPM_ci_upper_bound', 'FPKM_ci_lower_bound' and

        'FPKM_ci_upper_bound' are lower(l) and upper(u) bounds of 95%

        credibility intervals for TPM and FPKM values. The bounds are

        inclusive (i.e. [l, u]).

 

        'TPM_coefficient_of_quartile_variation' and

        'FPKM_coefficient_of_quartile_variation' are coefficients of quartile

        variation (CQV) for TPM and FPKM values. CQV is a robust way of

        measuring the ratio between the standard deviation and the mean. It is

        defined as

 

        CQV := (Q3 - Q1) / (Q3 + Q1),

 

        where Q1 and Q3 are the first and third quartiles.

 

    sample_name.genes.results

        File containing gene level expression estimates. The first line

        contains column names separated by the tab character. The format of

        each line in the rest of this file is:

 

        gene_id transcript_id(s) length effective_length expected_count TPM

        FPKM [posterior_mean_count posterior_standard_deviation_of_count

        pme_TPM pme_FPKM TPM_ci_lower_bound TPM_ci_upper_bound

        TPM_coefficient_of_quartile_variation FPKM_ci_lower_bound

        FPKM_ci_upper_bound FPKM_coefficient_of_quartile_variation]

 

        Fields are separated by the tab character. Fields within "" are

        optional. They will not be presented if neither '--calc-pme' nor

        '--calc-ci' is set.

 

        'transcript_id(s)' is a comma-separated list of transcript_ids

        belonging to this gene. If no gene information is provided, 'gene_id'

        and 'transcript_id(s)' are identical (the 'transcript_id').

 

        A gene's 'length' and 'effective_length' are defined as the weighted

        average of its transcripts' lengths and effective lengths (weighted by

        'IsoPct'). A gene's abundance estimates are just the sum of its

        transcripts' abundance estimates.

 

    sample_name.alleles.results

        Only generated when the RSEM references are built with allele-specific

        transcripts.

 

        This file contains allele level expression estimates for

        allele-specific expression calculation. The first line contains column

        names separated by the tab character. The format of each line in the

        rest of this file is:

 

        allele_id transcript_id gene_id length effective_length expected_count

        TPM FPKM AlleleIsoPct AlleleGenePct [posterior_mean_count

        posterior_standard_deviation_of_count pme_TPM pme_FPKM

        AlleleIsoPct_from_pme_TPM AlleleGenePct_from_pme_TPM

        TPM_ci_lower_bound TPM_ci_upper_bound

        TPM_coefficient_of_quartile_variation FPKM_ci_lower_bound

        FPKM_ci_upper_bound FPKM_coefficient_of_quartile_variation]

 

        Fields are separated by the tab character. Fields within "[]" are

        optional. They will not be presented if neither '--calc-pme' nor

        '--calc-ci' is set.

 

        'allele_id' is the allele-specific name of this allele-specific

        transcript.

 

        'AlleleIsoPct' stands for allele-specific percentage on isoform level.

        It is the percentage of this allele-specific transcript's abundance

        over its parent transcript's abundance. If its parent transcript has

        only one allele variant form, this field will be set to 100.

 

        'AlleleGenePct' stands for allele-specific percentage on gene level.

        It is the percentage of this allele-specific transcript's abundance

        over its parent gene's abundance.

 

        'AlleleIsoPct_from_pme_TPM' and 'AlleleGenePct_from_pme_TPM' have

        similar meanings. They are calculated based on posterior mean

        estimates.

 

        Please note that if this file is present, the fields 'length' and

        'effective_length' in 'sample_name.isoforms.results' should be

        interpreted similarly as the corresponding definitions in

        'sample_name.genes.results'.

 

    sample_name.transcript.bam

        Only generated when --no-bam-output is not specified.

 

        'sample_name.transcript.bam' is a BAM-formatted file of read

        alignments in transcript coordinates. The MAPQ field of each alignment

        is set to min(100, floor(-10 * log10(1.0 - w) + 0.5)), where w is the

        posterior probability of that alignment being the true mapping of a

        read. In addition, RSEM pads a new tag ZW:f:value, where value is a

        single precision floating number representing the posterior

        probability. Because this file contains all alignment lines produced

        by bowtie or user-specified aligners, it can also be used as a

        replacement of the aligner generated BAM/SAM file.

 

    sample_name.transcript.sorted.bam and

    sample_name.transcript.sorted.bam.bai

        Only generated when --no-bam-output is not specified and

        --sort-bam-by-coordinate is specified.

 

        'sample_name.transcript.sorted.bam' and

        'sample_name.transcript.sorted.bam.bai' are the sorted BAM file and

        indices generated by samtools (included in RSEM package).

 

    sample_name.genome.bam

        Only generated when --no-bam-output is not specified and

        --output-genome-bam is specified.

 

        'sample_name.genome.bam' is a BAM-formatted file of read alignments in

        genomic coordinates. Alignments of reads that have identical genomic

        coordinates (i.e., alignments to different isoforms that share the

        same genomic region) are collapsed into one alignment. The MAPQ field

        of each alignment is set to min(100, floor(-10 * log10(1.0 - w) +

        0.5)), where w is the posterior probability of that alignment being

        the true mapping of a read. In addition, RSEM pads a new tag

        ZW:f:value, where value is a single precision floating number

        representing the posterior probability. If an alignment is spliced, a

        XS:A:value tag is also added, where value is either '+' or '-'

        indicating the strand of the transcript it aligns to.

 

    sample_name.genome.sorted.bam and

    sample_name.genome.sorted.bam.bai

        Only generated when --no-bam-output is not specified, and

        --sort-bam-by-coordinate and --output-genome-bam are specified.

 

        'sample_name.genome.sorted.bam' and

        'sample_name.genome.sorted.bam.bai' are the sorted BAM file and

        indices generated by samtools (included in RSEM package).

 

    sample_name.time

        Only generated when --time is specified.

 

        It contains time (in seconds) consumed by aligning reads, estimating

        expression levels and calculating credibility intervals.

 

    sample_name.log

        Only generated when --alignments is not specified.

 

        It captures alignment statistics outputted from the user-specified

        aligner.

 

    sample_name.stat

        This is a folder instead of a file. All model related statistics are

        stored in this folder. Use 'rsem-plot-model' can generate plots using

        this folder.

 

        'sample_name.stat/sample_name.cnt' contains alignment statistics. The

        format and meanings of each field are described in

        'cnt_file_description.txt' under RSEM directory.

 

        'sample_name.stat/sample_name.model' stores RNA-Seq model parameters

        learned from the data. The format and meanings of each filed of this

        file are described in 'model_file_description.txt' under RSEM

        directory.

 

        The following four output files will be generated only by

        prior-enhanced RSEM

 

        - 'sample_name.stat/sample_name_prsem.all_tr_features'

          It stores isofrom features for deriving and assigning pRSEM prior.

          The first line is a header and the rest is one isoform per line. The

          description for each column is:

 

          * trid: transcript ID from input annotation

 

          * geneid: gene ID from input anntation

 

          * chrom: isoform's chromosome name

 

          * strand: isoform's strand name

 

          * start: isoform's end with the lowest genomic loci

 

          * end: isoform's end with the highest genomic loci

 

          * tss_mpp: average mappability of [TSS-500bp, TSS+500bp], where TSS

            is isoform's transcription start site, i.e. 5'-end

 

          * body_mpp: average mappability of (TSS+500bp, TES-500bp), where TES

            is isoform's transcription end site, i.e. 3'-end

 

          * tes_mpp: average mappability of [TES-500bp, TES+500bp]

 

          * pme_count: isoform's fragment or read count from RSEM's posterior

            mean estimates

 

          * tss: isoform's TSS loci

 

          * tss_pk: equal to 1 if isoform's [TSS-500bp, TSS+500bp] region

            overlaps with a RNA Pol II peak; 0 otherwise

 

          * is_training: equal to 1 if isoform is in the training set where

            Pol II prior is learned; 0 otherwise

 

        - 'sample_name.stat/sample_name_prsem.all_tr_prior'

          It stores prior parameters for every isoform. This file does not

          have a header. Each line contains a prior parameter and an isoform's

          transcript ID delimited by ` # `.

 

        - 'sample_name.stat/sample_name_uniform_prior_1.isoforms.results'

          RSEM's posterior mean estimates on the isoform level with an initial

          pseudo-count of one for every isoform. It is in the same format as

          the 'sample_name.isoforms.results'.

 

        - 'sample_name.stat/sample_name_uniform_prior_1.genes.results'

          RSEM's posterior mean estimates on the gene level with an initial

          pseudo-count of one for every isoform. It is in the same format as

          the 'sample_name.genes.results'.

 

        When learning prior from multiple external data sets in prior-enhanced

        RSEM, two additional output files will be generated.

 

        - 'sample_name.stat/sample_name.pval_LL'

          It stores a p-value and a log-likelihood. The p-value indicates

          whether the combination of multiple complementary data sets is

          informative for RNA-seq quantification. The log-likelihood shows how

          well pRSEM's Dirichlet-multinomial model fits the read counts of

          partitioned training set isoforms.

 

        - 'sample_name.stat/sample_name.lgt_mdl.RData'

          It stores an R object named 'glmmdl', which is a logistic regression

          model on the training set isoforms and multiple external data sets.

 

        In addition, extra columns will be added to

        'sample_name.stat/all_tr_features'

 

        * is_expr: equal to 1 if isoform has an abundance >= 1 TPM and a

          non-zero read count from RSEM's posterior mean estimates; 0

          otherwise

 

        * "$external_data_set_basename": log10 of external data's signal at

          [TSS-500, TSS+500]. Signal is the number of reads aligned within

          that interval and normalized to RPKM by read depth and interval

          length. It will be set to -4 if no read aligned to that interval.

 

          There are multiple columns like this one, where each represents an

          external data set.

 

        * prd_expr_prob: predicted probability from logistic regression model

          on whether this isoform is expressed or not. A probability higher

          than 0.5 is considered as expressed

 

        * partition: group index, to which this isoforms is partitioned

 

        * prior: prior parameter for this isoform

 

EXAMPLES

    Assume the path to the bowtie executables is in the user's PATH

    environment variable. Reference files are under '/ref' with name

    'mouse_125'.

 

    1) '/data/mmliver.fq', single-end reads with quality scores. Quality

    scores are encoded as for 'GA pipeline version >= 1.3'. We want to use 8

    threads and generate a genome BAM file. In addition, we want to append

    gene/transcript names to the result files:

 

     rsem-calculate-expression --phred64-quals \

                               -p 8 \

                               --append-names \

                               --output-genome-bam \

                               /data/mmliver.fq \

                               /ref/mouse_125 \

                               mmliver_single_quals

 

    2) '/data/mmliver_1.fq' and '/data/mmliver_2.fq', stranded paired-end

    reads with quality scores. Suppose the library is prepared using TruSeq

    Stranded Kit, which means the first mate should map to the reverse strand.

    Quality scores are in SANGER format. We want to use 8 threads and do not

    generate a genome BAM file:

 

     rsem-calculate-expression -p 8 \

                               --paired-end \

                               --strandedness reverse \

                               /data/mmliver_1.fq \

                               /data/mmliver_2.fq \

                               /ref/mouse_125 \

                               mmliver_paired_end_quals

 

    3) '/data/mmliver.fa', single-end reads without quality scores. We want to

    use 8 threads:

 

     rsem-calculate-expression -p 8 \

                               --no-qualities \

                               /data/mmliver.fa \

                               /ref/mouse_125 \

                               mmliver_single_without_quals

 

    4) Data are the same as 1). This time we assume the bowtie executables are

    under '/sw/bowtie'. We want to take a fragment length distribution into

    consideration. We set the fragment length mean to 150 and the standard

    deviation to 35. In addition to a BAM file, we also want to generate

    credibility intervals. We allow RSEM to use 1GB of memory for CI

    calculation:

 

     rsem-calculate-expression --bowtie-path /sw/bowtie \

                               --phred64-quals \

                               --fragment-length-mean 150.0 \

                               --fragment-length-sd 35.0 \

                               -p 8 \

                               --output-genome-bam \

                               --calc-ci \

                               --ci-memory 1024 \

                               /data/mmliver.fq \

                               /ref/mouse_125 \

                               mmliver_single_quals

 

    5) '/data/mmliver_paired_end_quals.bam', BAM-formatted alignments for

    paired-end reads with quality scores. We want to use 8 threads:

 

     rsem-calculate-expression --paired-end \

                               --alignments \

                               -p 8 \

                               /data/mmliver_paired_end_quals.bam \

                               /ref/mouse_125 \

                               mmliver_paired_end_quals

 

    6) '/data/mmliver_1.fq.gz' and '/data/mmliver_2.fq.gz', paired-end reads

    with quality scores and read files are compressed by gzip. We want to use

    STAR to aligned reads and assume STAR executable is '/sw/STAR'. Suppose we

    want to use 8 threads and do not generate a genome BAM file:

 

     rsem-calculate-expression --paired-end \

                               --star \

                               --star-path /sw/STAR \

                               --gzipped-read-file \

                               --paired-end \

                               -p 8 \

                               /data/mmliver_1.fq.gz \

                               /data/mmliver_2.fq.gz \

                               /ref/mouse_125 \

                               mmliver_paired_end_quals

 

    7) In the above example, suppose we want to run prior-enhanced RSEM

    instead. Assuming we want to learn priors from a ChIP-seq peak file

    '/data/mmlive.narrowPeak.gz':

 

     rsem-calculate-expression --star \

                               --star-path /sw/STAR \

                               --gzipped-read-file \

                               --paired-end \

                               --calc-pme \

                               --run-pRSEM \

                               --chipseq-peak-file /data/mmliver.narrowPeak.gz \

                               -p 8 \

                               /data/mmliver_1.fq.gz \

                               /data/mmliver_2.fq.gz \

                               /ref/mouse_125 \

                               mmliver_paired_end_quals

 

    8) Similar to the example in 7), suppose we want to use the partition

    model 'pk_lm2nopk' (partitioning isoforms by Pol II TSS peak first and

    then partitioning 'no TSS peak' isoforms into two bins by a linear

    regression model), and we want to partition isoforms by RNA Pol II's

    ChIP-seq read files '/data/mmliver_PolIIRep1.fq.gz' and

    '/data/mmliver_PolIIRep2.fq.gz', and the control ChIP-seq read files

    '/data/mmliver_ChIPseqCtrl.fq.gz'. Also, assuming Bowtie's executables are

    under '/sw/bowtie/':

 

     rsem-calculate-expression --star \

                               --star-path /sw/STAR \

                               --gzipped-read-file \

                               --paired-end \

                               --calc-pme \

                               --run-pRSEM \

                               --chipseq-target-read-files /data/mmliver_PolIIRep1.fq.gz,/data/mmliver_PolIIRep2.fq.gz \

                               --chipseq-control-read-files /data/mmliver_ChIPseqCtrl.fq.gz \

                               --partition-model pk_lm2nopk \

                               --bowtie-path /sw/bowtie \

                               -p 8 \

                               /data/mmliver_1.fq.gz \

                               /data/mmliver_2.fq.gz \

                               /ref/mouse_125 \

                               mmliver_paired_end_quals

 

    9) Similar to the example in 8), suppose we want to derive prior from four

    histone modification ChIP-seq read data sets: '/data/H3K27Ac.fastq.gz',

    '/data/H3K4me1.fastq.gz', '/data/H3K4me2.fastq.gz', and

    '/data/H3K4me3.fastq.gz'. Also, assuming Bowtie's executables are under

    '/sw/bowtie/':

 

     rsem-calculate-expression --star \

                               --star-path /sw/STAR \

                               --gzipped-read-file \

                               --paired-end \

                               --calc-pme \

                               --run-pRSEM \

                               --partition-model cmb_lgt \

                               --chipseq-read-files-multi-targets /data/H3K27Ac.fastq.gz,/data/H3K4me1.fastq.gz,/data/H3K4me2.fastq.gz,/data/H3K4me3.fastq.gz \

                               --bowtie-path /sw/bowtie \

                               -p 8 \

                               /data/mmliver_1.fq.gz \

                               /data/mmliver_2.fq.gz \

                               /ref/mouse_125 \

                               mmliver_paired_end_quals

 

rsem-generate-data-matrix

$ rsem-generate-data-matrix

Usage: rsem-generate-data-matrix sampleA.[alleles/genes/isoforms].results sampleB.[alleles/genes/isoforms].results ... > output_name.matrix

All result files should have the same file type. The 'expected_count' columns of every result file are extracted to form the data matrix.

他にも多くのコマンドがある。

> rsem-

$ rsem-

rsem-bam2readdepth                                     rsem-extract-reference-transcripts                     rsem-get-unique                                        rsem-preref                                            rsem-scan-for-paired-end-reads

rsem-bam2wig                                           rsem-for-ebseq-calculate-clustering-info               rsem-gff3-to-gtf                                       rsem-refseq-extract-primary-assembly                   rsem-simulate-reads

rsem-build-read-index                                  rsem-for-ebseq-find-DE                                 rsem-parse-alignments                                  rsem-run-ebseq                                         rsem-synthesis-reference-transcripts

rsem-calculate-credibility-intervals                   rsem-for-ebseq-generate-ngvector-from-clustering-info  rsem-plot-model                                        rsem-run-em                                            rsem-tbam2gbam

rsem-calculate-expression                              rsem-gen-transcript-plots                              rsem-plot-transcript-wiggles                           rsem-run-gibbs                                         

rsem-control-fdr                                       rsem-generate-data-matrix                              rsem-prepare-reference                                 rsem-sam-validator                                     

 

 

実行方法

1、indexing

GTF ファイルとゲノムのfasta、最後にindex名を指定する(1と同じにする)。ランの過程でbowtie2/star/hisat2のindexも作成される( この例ではxxx.index.〜)。--gff3を使えばGTFの代わりにGFF3のアノテーションを与えることもできる。

#bowtie2
rsem-prepare-reference --gtf genome.gtf --bowtie2 --bowtie2-path <path>/<to>/<your>/<bowtie2-path> -p 20 genome.fa bowtie2_index

#star
rsem-prepare-reference --gtf genome.gtf --star --star-path <path>/<to>/<your>/<star-path> -p 20 genome.fa star_index

#hisat2
rsem-prepare-reference --gtf genome.gtf --hisat2-hca --hisat2-path <path>/<to>/<your>/<hisat2-path> -p 20 genome.fa hisat2_index
  • --gtf   If this option is on, RSEM assumes that 'reference_fasta_file(s)'  contains the sequence of a genome, and will extract transcript reference sequences using the gene annotations specified in <file>,  which should be in GTF format.
    If this and '--gff3' options are off, RSEM will assume  'reference_fasta_file(s)' contains the reference transcripts. In this   case, RSEM assumes that name of each sequence in the Multi-FASTA files  is its transcript_id. (Default: off)
  • --gff3   The annotation file is in GFF3 format instead of GTF format. RSEM will first convert it to GTF format with the file name 'reference_name.gtf'. Please make sure that 'reference_name.gtf' does not exist. (Default: off)
  • --bowtie2   Use Bowtie 2 instead of Bowtie to align reads. Since currently RSEM
    does not handle indel, local and discordant alignments, the Bowtie2 parameters are set in a way to avoid those alignments. In particular, we use options '- sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1' by default. The last parameter of '--score-min', '-0.1', is the negative of maximum mismatch rate. This rate can be set by option '--bowtie2-mismatch-rate'. If reads are paired-end, we additionally use options '--no-mixed' and '--no-discordant'. (Default: off)
  • --star   Use STAR to align reads. Alignment parameters are from ENCODE3's
    STAR-RSEM pipeline. To save computational time and memory resources,
    STAR's Output BAM file is unsorted. It is stored in RSEM's temporary
    directory with name as 'sample_name.bam'. Each STAR job will have
    its own private copy of the genome in memory. (Default: off)
  • --hisat2-hca   Use HISAT2 to align reads to the transcriptome according to Human. Cell Atlast SMART-Seq2 pipeline. In particular, we use HISAT
    parameters "-k 10 --secondary --rg-id=$sampleToken --rg
    SM:$sampleToken --rg LB:$sampleToken --rg PL:ILLUMINA --rg
    PU:$sampleToken --new-summary --summary-file $sampleName.log
    --met-file $sampleName.hisat2.met.txt --met 5 --mp 1,1 --np 1
    --score-min L,0,-0.1 --rdg 99999999,99999999 --rfg 99999999,99999999
    --no-spliced-alignment --no-softclip --seed 12345". If inputs are
    paired-end reads, we additionally use parameters "--no-mixed
    --no-discordant". (Default: off)
  • --transcript-to-gene-map    Use information from <file> to map from transcript (isoform) ids to  gene ids. Each line of <file> should be of the form:  gene_id transcript_id with the two fields separated by a tab character. If you are using a GTF file for the "UCSC Genes" gene set from the UCSC Genome Browser, then the "knownIsoforms.txt" file (obtained from the "Downloads" section of the UCSC Genome Browser site) is of this format.  If this option is off, then the mapping of isoforms to genes depends on whether the '--gtf' option is specified. If '--gtf' is specified,  then RSEM uses the "gene_id" and "transcript_id" attributes in the GTF  file. Otherwise, RSEM assumes that each sequence in the reference sequence files is a separate gene.  (Default: off)
  • --prep-pRSEM    A Boolean indicating whether to prepare reference files for pRSEM,   including building Bowtie indices for a genome and selecting training set isoforms. The index files will be used for aligning ChIP-seq reads   in prior-enhanced RSEM and the training set isoforms will be used for   learning prior. A path to Bowtie executables and a mappability file in  bigWig format are required when this option is on. Currently, Bowtie2  is not supported for prior-enhanced RSEM. (Default: off)
  • --mappability-bigwig-file   Full path to a whole-genome mappability file in bigWig format. This   file is required for running prior-enhanced RSEM. It is used for  selecting a training set of isoforms for prior-learning. This file can be either downloaded from UCSC Genome Browser or generated by GEM (Derrien et al., 2012, PLoS One). (Default: "")

インデックスxxx.index~が出力される。注意点として、STARを使う場合はメモリが多い環境で実行すること(32GB程度では足りないことが多い)。メモリが少ない環境でSTAR indexを作成すると、genomeParameterファイルが出力されない事がある。このファイルがないと次のステップでエラーを起こす。

--transcript-to-gene-map”についてはこちらを参照。

How to get --transcript-to-gene-map <file> in RSEM?

 

2、rsem-calculateを実行。オプションの後にペアエンドfastq、index名(bowtie2の場合、step1のコマンドでbowtie2.indexと名前をつけているので”bowtie2.index”と指定)、出力ファイル名を指定する。gzipped fastqは使えないようなので解凍して指定する。

#bowtie2
rsem-calculate-expression --paired-end -p 20 --bowtie2 --bowtie2-path <path>/<to>/<your>/<bowtie2-path> \
sample1_R1.fq sample1_R2.fq bowtie2_index sample1

#star
rsem-calculate-expression --paired-end -p 20 --star --star-path <path>/<to>/<your>/<star-path> \
sample1_R1.fq sample1_R2.fq star_index sample1

#hisat2
rsem-calculate-expression --paired-end -p 20 --hisat2-hca --hisat2-path <path>/<to>/<your>/<hisat2-path> \
sample1_R1.fq sample1_R2.fq hisat2_index sample1

sample1.genes.resultsとsample1.isofomrs.resultsその他のファイルが出力される。

bamが最後の方で生じ、ディスクに保存される。このリファレンスはゲノムでは無く転写産物リファレンスである。転写産物をリファレンスにしてマッピングしているので、IGV等のビューアで見る際には間違えないようにしたい。STARならstar_index.transcripts.faになる(別の見方をすれば、転写産物に当てているので、bowtieではスプリットアラインメントはどうなっているのかとか、hisat2 等のスプリットアラインメントの距離の最大値等は気にしなくて良い)。ただbamはソートされていないので、ビューアで見るなら、samtools sortしてindexのbam.baiも作る必要はある。

 

3、複数の結果を統合 (expected countの統合*1)

#gene level
rsem-generate-data-matrix sample*genes.results

#transcripts(isoform) level
rsem-generate-data-matrix sample*isoforms.results

 expected read count(RSEMの確率的リードカウント値で正規化されていない生の値)の表が出力される。先頭行の名前だけ修正してiDEP(修正)などにロードすれば、すぐに結果を得ることができる。その時はまずPCAなどを行なって、妥当なグループになっているかどうか確認すること。

 

 

 追記

マッピングを別途実行し、その結果のbamファイルRSEMに提供することもできる。その場合はアラインメントオプション、特にリピートのマッピング設定に注意する。

例えばSTARの場合

1、indexing and mapping

#index
STAR --runMode genomeGenerate --genomeDir STAR_index --genomeFastaFiles genome.fa --sjdbGTFfile annotation.gtf

#mapping
STAR --runThreadN 12 --genomeDir STAR_index --readFilesIn pair_R1.fq.gz pair_R2.fq.gz --readFilesCommand zcat --quantMode TranscriptomeSAM --outSAMtype BAM SortedByCoordinate --genomeLoad NoSharedMemory --outFilterMultimapNmax 1 --outFileNamePrefix sample1

2、read count

bamファイルを指定する際は--alignmentsのフラグを立てる。

rsem-calculate-expression --alignments --paired-end -p 20 --bowtie2 input.bam bowtie2_index sample1
  • --alignments   Input file contains alignments in SAM/BAM/CRAM format. The exact file format will be determined automatically. (Default: off)

 

 

 

引用

RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome

Bo Li, Colin N Dewey

BMC Bioinformatics. 2011 Aug 4 

 

参考


 関連


*1 genes.results やisofomrs.resultsにはカウント値として確率的カウント値、TPMとFPKMの正規化されたカウント値の3つプリントされている。rsem-generate-data-matrixやrsem-generate-data-matrixは確率的カウント値を統合する。切り替えるオプションはないので、TPMやFPKMを取り出したいならその列だけ抜き出して横に繋ぐ。

#TPM
cut -f 6 sample1.genes.results > sample1
cut -f 6 sample2.genes.results > sample2
cut -f 6 sample3.genes.results > sample3
paste sample1 sample2 sample3 > TPM

#FPKM
cut -f 7 sample1.genes.results > sample1
cut -f 7 sample2.genes.results > sample2
cut -f 7 sample3.genes.results > sample3
paste sample1 sample2 sample3 > FPKM

 

追記

ベンチマークペーパー

2つの独立したデータセットを用いて、7つの競合パイプラインを評価した。性能は全般的に低く、2つの手法が明らかに劣っており、RSEMが残りの手法をわずかに上回っていた。


関連