一、比对原理

背景

1、短读长、非连续性的转录本结构、持续增长的通量是转录组实验的现状
2、高的比对错误率、低的比对速度、读长的限制（很多软件要求读长<200nt,并且不适合三代测序的比对读长）、比对的偏差、对于剪切行为的识别不敏感（很多基于DNA这样连续序列的比对软件只能通过提供剪切序列数据库或者只部分比对的策略完成剪接体的识别，这样会遗漏很多新型剪接体）是当前很多比对软件的弊端
3、目前对比对算法的要求为：1）可以在本就存在错配、插入或缺失的参考基因组中完成精确匹配 2）部分序列存在剪切行为，并融合形成剪切的RNA，对这部分序列能精确的匹配

原理

1、Seed search

1、核心思想为按顺序查找maximal mappable prefix (MMP???)
2、MMP定义：给定一个read R,read 上的位点为i,参考序列为G,那MMP(R,i,G)就定义为R从i位置开始与G实现匹配的最长子串
3、Seed search的过程为：假设存在一个没有错配的存在一个剪切位点的read，从read的第一个碱基开始检索MMP，直到无法匹配为止，定义为MMP1即seed1，之后继续重新匹配，定义为MMP2即seed2，直到将read匹配完。
4、这种策略可以精确找到经典的剪切位点，并可以找出非经典的剪切位点
5、同时该策略对存在错配或indel的read也可以进行比对

seed search

2、Clustering、stitching、scores

1、先选择一些锚，然后锚周围的seed进行clustering
2、同一个cluster的seed进行stitching，同时可以检测到嵌合体
3、最终对stitching的结果进行打分（依据匹配、错配、插入、缺失...），得分最高的stitching被选为该read的最优比对

二、比对方案

参考帮助文档：https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf
https://github.com/alexdobin/STAR

1、建索引

--runThreadN NumberOfThreads
--runMode genomeGenerate
--genomeDir /path/to/genomeDir
--genomeFastaFiles /path/to/genome/fasta1 /path/to/genome/fasta2 ..
--sjdbGTFfile /path/to/annotations.gtf
--sjdbOverhang ReadLength-1
#runThreadN：使用的线程数目
#runMode：未看懂，选择genomeGenerate即可
#genomeDir：用来输出索引的文件夹名字，需要有写入权限，最好是个空文件夹
#genomeFastaFiles:参考基因组名
#sjdbGTFfile：gtf文件名，可选项，STAR可以从中提取剪切位点，可以提高比对准确性，推荐！
#sjdbOverhang：如果read读长是100nt，那么这个参数值为100-1=99
#最终可以获得五个文件，binary genome sequence, suffix arrays, text chromosome names/lengths,
#splice junctions coordinates, and transcripts/genes information
注意：
1、使用的fasta和gtf要配套，因为不同网站上fasta和gtf的注释信息格式略有不同，不配套容易出错
2、如果参考基因组很小
For small genomes, the parameter --genomeSAindexNbases must to be scaled down, with a typical
value of min(14, log2(GenomeLength)/2 - 1). For example, for 1 megaBase genome, this is equal
to 9, for 100 kiloBase genome, this is equal to 7

2、比对

--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2 ]
--outFileNamePrefix /path/to/output/dir/prefix
--outSAMtype BAM SortedByCoordinate
--quantMode GeneCounts
#genomeDir:保存索引的目录
#readFilesIn：fastq所在的文件夹
#outFileNamePrefix：输出文件的路径及前缀
#outSAMtype：输出的bam文件的类型，可选值为BAM Unsorted（输出为排序的bam），BAM SortedByCoordinate（排过序的bam，与samtolls sort结果类似），BAM Unsorted SortedByCoordinate（两种文件都输出）
#quantMode：是否输出基因定量结果，默认为-（不输出），可选项为TranscriptomeSAM（output SAM/BAM alignments to transcriptome into a separate file）和GeneCounts（count reads per gene）或者TranscriptomeSAM GeneCounts（两者都输出），注意在只用quantMode参数的时候一定要保证提供了gtf文件，即--sjdbGTFfile是有值的！
注意：
1、star支持fasta和fastq格式的文件
2、如果read文件是压缩形式的，如果是*.gz, --readFilesCommand zcat或者 --readFilesCommand gunzip -c，如果是bzip2格式的， --readFilesCommand bunzip2 -c
3、对于readfilesin参数，支持多个测序文件一起比对
如果是单端测序的多测序文件一起比对，则每个文件需要用","隔开，如--readFilesIn sample1.fq,sample2.fq,sample3.fq     注意中间没有空格！！！
如果是双端测序的多测序文件比对，格式为S1read1,S2read1,S3read1 S1read2,S2read2,S3read2
即read1和read2内部用逗号隔开，read1和read2之间用空格隔开！
还有一些更专业的参数可以参考帮助文档！

3、输出

1、Log file
Log.out（一些进程日志，方便之后debug）
Log.progress.out（进程中的统计工作，每一分钟记录一次，包括当前比对率等信息）
Log.final.out（最终的比对结果统计文件，对之后的质控很有用！）
2、SAM文件（Aligned.out.sam）
3、Unsorted或sorted-by-coordinate的BAM文件
通过调整--outSAMtype 可以得到不同形式的bam文件Aligned.out.bam或Aligned.sortedByCoord.out.bam两类文件
4、splice junction
SJ.out.tab文件中包含了可信度很高的剪切位点信息，具体每一列的意义参考帮助文档