STAR 使用手册

Mapping RNA-seq Reads with STAR

基础用法：Mapping RNA-seq reads to the reference genome

• 通过Illumina 或者Ion Torrent测序后得到千万条相对较短的序列(reads)片段，这些片段是来自于RNA分子。
• 最基础的用法是将这些reads比对到参考基因组上，然后进行各种基础的下游分析，包括：
  • 计算转录/基因表达量
  • 差异基因表达
  • 发现新的剪切链接和转录本
  • 基因组间的信号可视化
• STAR对计算机的需求：
  • Unix、linux或Mac OS X操作系统
  • 10 × 内存，如果参考基因组大小为3GB，那么即内存最小需求为30GB，建议32GB以上的内存。
  • 储存空间>100GB

1. 软件准备与安装：略

注释基因组的准备：

• 注释文件GTF/GFF下载并解压(当然也可以在ensemble网站下载其他版本的GTF文件)

wget ftp://ftp.ensembl.org/pub/release-79/gtf/homo_sapiens/Homo_sapiens.GRCh38.79.gtf.gz
gunzip Homo_sapiens.GRCh38.79.gtf.gz

• 准备测序原始数据：略

2. 命令行(没有建立索引，input为压缩格式)

STAR \
--runThreadN 12 --genomeDir ~/star/genome/ \
--sjdbGTFfile ~/star/Homo_sapiens.GRCh38.79.gtf --sjdbOverhang 100 \
--readFilesIn ~/star/ENCFF001RFH.fastq.gz ~/star/ENCFF001RFG.fastq.gz \
--readFilesCommand zcat

• 注意：

  • 当fastq文件未压缩时，删除--readFilesCommand zcat选项

  • 可以不是用注释文件进行比对，删除--sjdbGTFfile~/star/Homo_sapiens.GRCh38.79.gtf --sjdbOverhang 100 options命令行（不推荐），在缺少注释文件时，可以使用2-pass方式进行mapping。

  • 对于双段测序数据 --sjdbGTFfile命令行后键入空格，然后输入read-1和read-2的位置，单端输入一个。

3. 比对过程中的状态信息

ar 31 01:34:01 ..... Started STAR run
Mar 31 01:34:49 ..... Finished GTF processing
Mar 31 01:37:10 ..... Finished inserting 1st pass junctions into genome
Mar 31 01:37:10 ..... Started mapping
Mar 31 01:54:53 ..... Finished successfull

4. 日志输出

过程日志将被记录并保存到Log.progress.out文件中，包含了每分钟显示reads被处理的数量，和一些统计信息：

cat Log.progress.out


Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped
Unmapped Unmapped Unmapped
M/hr number length unique length MMrate multi multi + MM short other
Mar 31 01:38:12 299.7 5161748 202 92.2% 201.0 0.3% 6.0% 0.1% 0.0% 1.7% 0.0%
Mar 31 01:39:12 356.2 12069587 202 92.2% 200.9 0.3% 6.0% 0.1% 0.0% 1.7% 0.0%
Mar 31 01:40:13 347.7 17674136 202 92.2% 200.9 0.3% 6.0% 0.1% 0.0% 1.7% 0.0%
......
Mar 31 01:54:38 356.7 103831542 202 92.1% 200.9 0.3% 6.0% 0.1% 0.0% 1.8% 0.0%
ALL DONE!

5. 程序运行完毕

程序运行完毕后将出现以下文件

Aligned.out.sam        # 比对后SAM文件，最重要
Log.final.out          # 运行mapping的summary
Log.out                # 包含各种运行时间的信息，通过这些信息可以进行Debug
Log.progress.out
SJ.out.tab             # 包含了高可信度的splice junction信息

可选的方案1：建立基因组索引

1. 准备工作

• 硬件需求（如前）

• 软件安装（如前）

• Input files下载

  wget ftp://ftp.ensembl.org/pub/release-79/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
  

2. 运行命令行

/STAR \ 
--runThreadN 12              # 线程数
--runMode genomeGenerate     # 运行模式
--genomeDir ./ \             # index输出
--genomeFastaFiles ./Homo_sapiens.GRCh38.dna.primary_assembly.fa  \     # 参考基因组位置
--sjdbGTFfile ./Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf \     # 注释文件位置 
--sjdbOverhang 75                                                       # 测序长度-1

3. 状态信息输出

Mar 30 23:28:52 ..... Started STAR run
Mar 30 23:28:52 ... Starting to generate Genome files
Mar 30 23:29:56 ... starting to sort Suffix Array. This may take a long time...
Mar 30 23:30:17 ... sorting Suffix Array chunks and saving them to disk...
Mar 30 23:46:51 ... loading chunks from disk, packing SA...
Mar 30 23:48:57 ... writing Suffix Array to disk ...
Mar 30 23:49:20 ... Finished generating suffix array
Mar 30 23:49:20 ... starting to generate Suffix Array index...
Mar 31 00:09:28 ... writing SAindex to disk
Mar 31 00:09:29 ..... Finished successfully

4. 输出文件（现在版本有15个文件输出）

chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
Genome
genomeParameters.txt
Log.out
SA
SAindex

可选方案2：Mapping RNA-seq reads with 2-pass procedure

STAR进行2次mapping，第一次发现新的剪切并将其插入基因组索引，第二次会在mapping一次注释文件（包括在第一次发现的新的剪切），从而发现新剪切的灵敏度。特别是在缺少注释文件情况下（没有某些物种的注释信息），这种方式特别推荐。

1.准备工作

与基础用法需求一致

2. 命令行

STAR \
--runThreadN 12 \
--genomeDir ~/star/genome/ \
--sjdbGTFfile ~/star/Homo_sapiens.GRCh38.79.gtf \
--sjdbOverhang 100 \
--readFilesIn ~/star/ENCFF001RFH.fastq.gz ~/star/ENCFF001RFG.fastq.gz 
--readFilesCommand zcat \
--twopassMode Basic            # 这行命令激活2-pass

• 注意：在缺少注释文件，移除--sjdbGTFfile~/star/Homo_sapiens.GRCh38.79.gtf -- sjdbOverhang 100命令行即可。

可选方案3： Mapping reads and generating unsorted and

coordinate-sorted BAM files

BAM与SAM文件信息一致，但更小

1. 准备工作

• 硬件需求（如前）

• 软件安装（如前）

• Input files下载

2. 命令行

STAR \
--runThreadN 12 
--genomeDir ~/star/genome/ \
--sjdbGTFfile ~/star/Homo_sapiens.GRCh38.79.gtf 
--sjdbOverhang 100 \
--readFilesIn ~/star/ENCFF001RFH.fastq.gz ~/star/ENCFF001RFG.fastq.gz 
--readFilesCommand zcat \
--outSAMtype BAM SortedByCoordinate Unsorted  # 输出unsorted和sorted BAM文件

• 注意：--outSAMtype命令行激活Unsorted或者SortedByCoordinate，可选：
  • --outSAMtype BAM Unsorted
  • --outSAMtype BAM SortedByCoordinate

3. 输出文件

输出文件与基础用法相比:多了Aligned.out.sam或者Aligned.sortedByCoord.out.bam

可选方案4: Generating signal files for visualization on genome browsers for stranded RNA-seq data

整个基因组的RNAseq“信号”计算为交叉（即映射到）每个基因组位置（核苷酸）的读数的数量。这些信号可以通过基因组浏览器如：UCSC基因浏览器 (http://genome.ucsc.edu/) or IGV 浏览器 (https://www.broadinstitute.org/igv/)进行可视化。

1. 准备工作

与基础用法一致

• 从方案3的命令行产生信号文件Aligned.sortedByCoord.out.bam（如果不需要不排序的bam文件，仅使用--outSAMtype BAM SortedByCoordinate选项即可）

2.命令行（产生signal 文件）：

STAR
--runMode inputAlignmentsFromBAM \
--inputBAMfile Aligned.sortedByCoord.out.bam \
--outWigType bedGraph \
--outWigStrand Stranded

3. 输出文件

Signal.UniqueMultiple.str1.out.bg         
Signal.UniqueMultiple.str2.out.bg
Signal.Unique.str1.out.bg
Signal.Unique.str2.out.bg

注意：格式均为BedGraph格式

*.Unique.          # 仅包括来自唯一映射读取的信号 
*.UniqueMultiple.  # 包括来自唯一和多映射读取的信号 

信号由基因组两条链分别生成：*.str1. 和*.str2.，它们对应于不同的链，具体取决于文库制备协议。 为了 其中第一个读数位于 RNA 分子的相反链上的协议（例如 Illumina 链 Tru-Seq），*.str1. 对应于 (-) 链，*.str2.对应于 (+) 链。

可选方案5:Generating signal files for visualization on genome browsers for un-stranded RNA-seq data

1. 准备工作

与可选方案4相同

2. 命令行

STAR \
--runMode inputAlignmentsFromBAM \
--inputBAMfile Aligned.sortedByCoord.out.bam \
--outWigType bedGraph \
--outWigStrand Untranded 

3. 输出文件

Signal.UniqueMultiple.str1.out.bg
Signal.Unique.str1.out.bg

与方案4相比，相当没有区分正负链了

可选方案6：Mapping RNA-seq reads and generating chimeric alignments to detect fusion transcripts and circular RNA

标准 STAR 输出仅包括基因组的线性比对，即它不包括嵌合（融合）比对，其中包括跨嵌合连接的读取和配对末端读取与非一致（嵌合）配对对齐。 该协议描述了如何产生嵌合比对并将它们输出到单独的文件中。

1. 准备工作

与基础用法一致

2. 命令行

STAR \
--runThreadN 12 
--genomeDir ~/star/genome/ \
--sjdbGTFfile ~/star/Homo_sapiens.GRCh38.79.gtf \
--sjdbOverhang 100 \
--readFilesIn ~/star/ENCFF001RFH.fastq.gz ~/star/ENCFF001RFG.fastq.gz \
--readFilesCommand zcat \
--chimSegmentMin 20       # 激活嵌合体输出。N是嵌合片段中每一个的最小允许长度

3. 输出文件

与方案4-5抑制，多了两个嵌合体文件

Chimeric.out.sam                 # sam格式
Chimeric.out.junction

可选方案7：Mapping RNA-seq reads, generating output in transcriptomic coordinates and using RSEM to quantify expression of transcripts and genes

量化 RNA 转录本和基因是 RNA-seq 数据分析中最重要的任务之一。 RSEM (Li and Dewey, 2011) 是一种流行的软件包，能够使用 RNA-seq 数据量化注释基因和转录本。 在此协议中，STAR 输出转录组坐标中的基因组比对，然后由 RSEM 用于产生转录/基因量化。

1.准备工作

如前

安装RSEM包

2.命令行

STAR \
--runThreadN 12 
--genomeDir ~/star/genome/ \
--sjdbGTFfile ~/star/Homo_sapiens.GRCh38.79.gtf \
--sjdbOverhang 100 \
--readFilesIn ~/star/ENCFF001RFH.fastq.gz ~/star/ENCFF001RFG.fastq.gz \
--readFilesCommand zcat \
--quantMode TranscriptomeSAM           # 转录本定量

3.输出文件

较4-5多了Aligned.toTranscriptome.out.bam文件，注意此文件是带注释的转录序列，而基因组 SAM/BAM 文件是基因组序列（染色体）

4.准备 RSEM参考文件

rsem-prepare-reference --gtf ~/star/Homo_sapiens.GRCh38.79.gtf \
~/star/genome/Homo_sapiens.GRCh38.dna.primary_assembly.fa ./ref

5. 运行RSEM对STAR 转录组BAM文件进行定量

rsem-calculate-expression --bam --no-bam-output -p 12 --paired-end --forwardprob 0 \
~/star/alt_rsem/Aligned.toTranscriptome.out.bam ~/star/rsem_ref/ref ~/star/
alt_rsem/Quant \
>& ~/star/alt_rsem/rsem.log

6. RSEM的输出文件

包括Quant.isoforms.results和 Quant.genes.results文件，转录本和基因表达水平的两个文件

可选方案8：Mapping RNA-seq reads and running Cufflinks to

assemble and quantify transcripts for stranded RNA-seq data

可选方案9：Mapping RNA-seq reads and running Cufflinks to

assemble and quantify transcripts for un-stranded RNA-seq data

理解

对于STAR 比对后的结果分析：

会出现一个Log.final.out文件，里面可以查看各种信息，其中最为重要的信息一般为：Uniquely mapped reads %”或者mapping rate，好的文库比对率超过90%或80%以上，当比对率<50%时可能出现的问题包括

• 去除核小体RNA不够充分。在RNAseq建库时，提取总RNA后需要去除rRNA，方法包括：
  • ribo-depletion techniques
    • 一种是通过DNA探针或者RNA探针杂交捕获rRNA再用磁珠吸附去除的方法，这种方法称为Ribo-Zero-seq
    • 另一种是利用特异性核酸酶处理的方法提取RNA，称为RNA酶消化法
  • Poly-A+ selection
    • 通过对多腺苷聚集酸化polyA+的RNA进行富集并筛选
  • poly-dT reverse transcription priming

由于rRNA包含许多的同源相似序列，因此RNAseq的reads将同一个基因上，导致multi-mapping reads% >15%的情形。

• 低测序质量：可以通过Mismatch rate per base, %,Deletion rate per base, Insertion rate per bas进行测序错误率进行评估。当然，这些指标也可能代表某些基因型变体。一般Illumina测序mismatch error rate应小于5%，indel error rates小于0.05%，高错误率暗示较低的测序质量。另外Average mapped length与相应Average input read length显著降低也暗示低的测序质量。

建议通过绘制 FASTQ 文件中“质量分数”的分布来评估测序质量。 测序质量差可能会导致未映射读数的数量增加。

• 外源性 RNA/DNA 污染。不同物种的RNA/DNA杂合子不能在参考基因组上进行mapping，因此将会增加unmapped的比例。如果测序质量较好，然而，% of reads unmapped: too short”或% of reads unmapped: other（>15%）可能表明存在外源性 RNA/DNA 污染 。

建议将几个未映射的reads BLAST 到 NCBI 序列数据库，以识别可能的污染源。

• 计算处理问题。通常，计算处理问题会导致运行失败并且不会生成映射统计输出（请参阅故障排除部分）。 但是，在某些情况下，映射作业将成功完成，从而产生极低的映射率 (<5%)。 一些典型的错误包括 (i) 使用错误的物种基因组索引目录； (ii) 使用与 read-1 和 read-2 相同的文件