star: 比对RNA-seq到参考基因组获取剪切事件

# 安装

## \## 获取star源代码
# Get latest STAR source from releases
wget https://github.com/alexdobin/STAR/archive/2.7.10b.tar.gz
tar -xzf 2.7.10b.tar.gz
cd STAR-2.7.10b

# Alternatively, get STAR source using git
# git clone https://github.com/alexdobin/STAR.git

## 编译

# Compile
cd STAR/source
make STAR

## 或者使用conda

  • 上面从源码编译的方法太慢,或访问github网络不稳定情况下,可以使用conda;
conda install -c bioconda star

## 检查安装是否成功

STAR
Usage: STAR  [options]... --genomeDir /path/to/genome/index/   --readFilesIn R1.fq R2.fq
Spliced Transcripts Alignment to a Reference (c) Alexander Dobin, 2009-2022

STAR version=2.7.10b
STAR compilation time,server,dir=2022-11-01T09:53:26-04:00 :/home/dobin/data/STAR/STARcode/STAR.master/source
For more details see:
<https://github.com/alexdobin/STAR>
<https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf>

To list all parameters, run STAR --help

# star 基本流程

## 建立基因组索引

## 比对reads到基因组

  • 基于构建的比对索引,将reads(fasta和fastq格式)比对到基因组;
  • 输出文件可以有:alignments (SAM/BAM), mapping summary statistics, splice junctions, unmapped reads, signal (wiggle) tracks etc

# 产生基因组索引

## 基本命令:

STAR --runThreadN 30 --runMode genomeGenerate \
 --genomeDir /dataBase/hg38StarIndex \
 --genomeFastaFiles /dataBase/GRCh38.p13.genome.fa \
 --sjdbGTFfile /dataBase/gencode.v42.annotation.gtf \
 --sjdbOverhang ReadLength-1

---runThreadN: 线程
--runMode:genomeGenerate 为产生基因组索引
--genomeDir:index文件输出位置
--genomeFastaFiles /path/to/genome/fasta1 /path/to/genome/fasta2 ... 基因组fa文件
--sjdbGTFfile:基因组注释文件
--sjdbOverhang ReadLength-1 #已经注释的剪切位点两端的序列长度,一般是read 长度-1,当read 长度是变化的时候,就用最大的read 长度-1

  • 最终生成的基因组文件包括binary genome sequence, suffix arrays, text chromosome names/lengths, splice junctions coordinates, and transcripts/genes information.

## 基因组序列文件和注释文件的选择

### 基因组序列文件: chromosomes/scaffolds/patches 应该怎么选择?

  • star 建议包含主要的染色体(e.g., for human chr1-22,chrX,chrY,chrM,) 和un-place, un-localized scaffolds. 通常,un-place, un-localized scaffolds仅让基因组长度增加了几个MegaBase,然而,大量reads比对到这些scaffolds 上核糖体RNA(rRNA)repeats区域。如果基因组中不包括这些scaffolds,这些reads 将会被视为unmapped,或者错位匹配到基因组其他位置。正常情况下,基因组不应该有patches 和alternative haplotypes 序列。
  • 基因组序列获取方式:
    • ENSEMBL: files marked with .dna.primary.assembly, such as: ftp://ftp.ensembl. org/pub/release-77/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_ assembly.fa.gz
    • GENCODE: files marked with PRI (primary). Strongly recommended for mouse and human: http://www.gencodegenes.org/.

### 基因组注释文件:

  • star 建议包含注释文件应该包含最全的注释。但是要注意GTF与fasta中染色体的名字要一致。UCSC 采用chr1, chr2, ... , ENSEMBL 采用1, 2, ... 命名, the ENSEMBL 和UCSC的 FASTA 和GTF 文件如果要混着用,就需要考虑重命名。
  • 一般情况下,使用--sjdbGTFfile GTF文件,如果使用GFF文件的话,需要设置--sjdbGTFtagExonParentTranscript Parent, 默认情况为--sjdbGTFtagExonParentTranscript transcript id

### 使用剪切注释数据

  • --sjdbFileChrStartEnd /path/to/sjdbFile.txt可以添加剪切注释数据
    文件中包含信息:
Chr \tab Start \tab End \tab Strand=+/-/.

Start 和End是内含子第一个和最后一个碱基(1-based chromosome coordinates);STAR可以从sjdbFile.txt和GTF文件中获取剪切事件。

### 较小的基因组

对于较小的基因组,参数--genomeSAindexNbases必须按比例缩小,为min(14,log2(GenomeLength)/2-1)。例如,对于1megaBase基基因组,使用9,对于100000基基因组,使用7。

### 基因组有许多序列片段

一个基因组有>5,000 chro- somes/scaffolds时,可能需要减少--genomeChrBinNb位以减少RAM消耗。建议使用以下缩放比例:--genomeChrBinNbits = min(18,log2[max(GenomeLength/NumberOfReferences, ReadLength)])。3 gigaBase基因组,有100000 条chromosomes/scaffolds,,可以设置为15。

# 比对

## 基本命令:

STAR --runThreadN 20 --genomeDir /dataBase/hg38StarIndex \
 --runMode alignReads \
 --quantMode TranscriptomeSAM GeneCounts \
 --readFilesIn  R1.fastq.gz R2.fastq.gz \
 --outSAMtype BAM SortedByCoordinate \
 --readFilesCommand gunzip -c \
 --outFileNamePrefix test

---runThreadN: 线程
--genomeDir:基因组index位置
--runMode:alignReads 比对基因组,与前面的genomeGenerate 不同;
--quantMode:定量类型;TranscriptomeSAM:输出转录组的比对情况到一个单独的文件;GeneCounts:每个基因的reads数目统计
--readFilesIn: 测序数据,可以是fasta和fastq格式的
--outSAMtype: 输出文件格式
--readFilesCommand: 测序文件是压缩文件的话,解压命令
--outFileNamePrefix: 输出文件前缀

# 比对结果

testLog.out: 程序运行日志
testLog.progress.out: 程序运行中实时统计
testLog.final.out: 比对结束之后的,比对情况统计
testAligned.sortedByCoord.out.bam:排序之后的比对结果    
testSJ.out.tab:剪切事件
testAligned.toTranscriptome.out.bam:转录组比对结果  
testReadsPerGene.out.tab: 每个基因read计数情况
  • testLog.out: 包含程序运行的详细过程
  • testLog.progress.out: 程序运行中实时统计, 处理的reads和mapped上的reads;每分钟更新一次
Time    Speed        Read     Read   Mapped   Mapped   Mapped   Mapped Unmapped Unmapped Unmapped Unmapped
                    M/hr      number   length   unique   length   MMrate    multi   multi+       MM    short    other
Jan 01 15:02:35    426.6     8295517       79    35.8%     75.8     0.3%     7.5%     1.9%     0.0%    54.6%     0.2%
ALL DONE!
  • TIA1Log.final.out:比对结束之后的,比对情况统计;STAR将双端测序的一对reads记为一个read,samtools flagstat/idxstats是将一对reads分开计算成两个read。这个文件中统计的read就是一对reads。并且具有唯一匹配和多个匹配的read都会分别统计。
Started job on |       Jan 01 15:00:50
                             Started mapping on |       Jan 01 15:01:25
                                    Finished on |       Jan 01 15:02:35
       Mapping speed, Million of reads per hour |       426.63

                          Number of input reads |       8295517
                      Average input read length |       79
                                    UNIQUE READS:
                   Uniquely mapped reads number |       2971995
                        Uniquely mapped reads % |       35.83%
                          Average mapped length |       75.76
                       Number of splices: Total |       81528
            Number of splices: Annotated (sjdb) |       35799
                       Number of splices: GT/AG |       41246
                       Number of splices: GC/AG |       1103
                       Number of splices: AT/AC |       277
               Number of splices: Non-canonical |       38902
                      Mismatch rate per base, % |       0.30%
                         Deletion rate per base |       0.03%
                        Deletion average length |       1.67
                        Insertion rate per base |       0.00%
                       Insertion average length |       1.29
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       618387
             % of reads mapped to multiple loci |       7.45%
        Number of reads mapped to too many loci |       154144
             % of reads mapped to too many loci |       1.86%
                                  UNMAPPED READS:
  Number of reads unmapped: too many mismatches |       0
       % of reads unmapped: too many mismatches |       0.00%
            Number of reads unmapped: too short |       4530903
                 % of reads unmapped: too short |       54.62%
                Number of reads unmapped: other |       20088
                     % of reads unmapped: other |       0.24%
                                  CHIMERIC READS:
                       Number of chimeric reads |       0
                            % of chimeric reads |       0.00%
  • testSJ.out.tab 剪切事件统计,star将剪切时间的第一个和最后一个碱基定义为内含子的第一个和最后一个碱基,值得注意的是,其他许多软件使用外显子碱基位置表示剪切事件。
chr1    15463   15623   1       1       0       0       1       21
第1列:染色体号
第2列:内含子第一个碱基 (1-based)
第3列:内含子最后一个碱基 (1-based)
第4列:链, (0: 为定义, 1: +, 2: -)
第5列:内含子基序:0:非规范;1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5: AT/AC, 6: GT/AT
第6列:0:未注释的,1:在剪切数据库中有注释的,2,在剪切数据库中有注释的,也在GTF文件中存在;
第7列:剪切事件上只有唯一匹配的reads
第8列:剪切事件上具有多次匹配的reads
第9列:maximum spliced alignment overhang
  • testReadsPerGene.out.tab: 每个基因read计数情况
N_unmapped      4705135 4705135 4705135
N_multimapping  618387  618387  618387
N_noFeature     2046544 2864860 2107430
N_ambiguous     78709   2231    32042
ENSG00000290825.1       0       0       0
ENSG00000223972.6       0       0       0
ENSG00000227232.5       1       0       1
column 1: gene ID
column 2: counts for unstranded RNA-seq
18
column 3: counts for the 1st read strand aligned with RNA (htseq-count option -s yes)
column 4: counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)

#原文

Github: Spliced Transcripts Alignment to a Reference
STARmanual

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 159,015评论 4 362
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,262评论 1 292
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 108,727评论 0 243
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,986评论 0 205
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,363评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,610评论 1 219
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,871评论 2 312
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,582评论 0 198
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,297评论 1 242
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,551评论 2 246
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,053评论 1 260
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,385评论 2 253
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,035评论 3 236
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,079评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,841评论 0 195
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,648评论 2 274
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,550评论 2 270

推荐阅读更多精彩内容