比对软件STAR的使用
在之前的学习和练习里，比对这一步我使用过bowtie2（DNA比对）和hisat2（RNA-seq比对），现在学习另一个很火的软件：STAR。STAR能够发现非典型拼接和嵌合（融合）转录本，并能够比对全长RNA序列（关于转录组比对STAR软件使用）。
STAR的官方使用说明书（PDF版）：https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf
因为我用的是服务器，所以这里就不贴如何在自己的电脑上安装STAR了。STAR的安装方法见：https://github.com/alexdobin/STAR
这里推荐一个网站，这里介绍了几种不同比对软件的使用方法（bowtie2, bwa, STAR, Tophat四种），非常的简单粗暴：here。因为这里我使用的是STAR，所以只对该网址的STAR部分的使用说明进行介绍。你也可以使用STAR的官方说明书，里面有更详细的参数选择。
在HOMER网站上，对几种比对软件进行的说明，可以看到STAR的速度是最快的，当然你得保证你有这么多的内存来运行：
bowtie : fast, works well
bowtie2 : fast, can perform local alignments too
BWA - Fast, allows indels, commonly used for genome/exome resequencing
Subread - Very fast, (also does splice alignment)
STAR - Extremely fast (also does splice alignment, requires at least 30 Gb memory)
虽然STAR的比对速度非常的快，但是这是建立在你有足够的内存基础上的，比如你如果使用小鼠或者人类的基因组来比对，那么至少需要30G内存来支持STAR的运行，所以一般的笔记本电脑都无法达到这个要求。（比如我的电脑只有8个G，基本啥啥都干不了）
Step 1 - Build a genome index构建索引
进行比对之前，需要对基因组构建索引。STAR软件可以使用hg19基因组的索引文件，所以你如果习惯使用hg19，那么就可以跳过这一步。这一步的代码格式是：
STAR --runMode genomeGenerate --runThreadN <# cpus> --genomeDir <genome output directory> --genomeFastaFiles <input Genome FASTA file>
其中：
--runThreadN是指你要用几个cpu来运行；
--genomeDir构建索引输出文件的目录；
--genomeFastaFiles你的基因组fasta文件所在的目录
比如说，我需要用GRCh38基因组来构建索引，我的代码（请忽略下面从srun到--pty之间的代码，从STAR开始，前面是我在服务器里运行代码时加的一些信息）：
$ srun --mem=128G --cpus-per-task=20 --time=8:00:00 --pty STAR \
--runMode genomeGenerate \
--runThreadN 40 \
--genomeDir ./ \ #生成文件存到当前文件夹
--genomeFastaFiles /gpfs/share/apps/iGenomes/Homo_sapiens/NCBI/GRCh38/Sequence/WholeGenomeFasta/genome.fa \
--sjdbGTFfile /gpfs/share/apps/iGenomes/Homo_sapiens/NCBI/GRCh38/Annotation/Genes.gencode/genes.gtf
#上面我加上了另一个参数--sjdbGTFfile，指示的是基因组的注释文件位置，你也可以不加这一个，但是在后续比对的时候需要指定基因组注释文件，否则会报错
Jun 28 09:42:40 ..... started STAR run
Jun 28 09:42:40 ... starting to generate Genome files
Jun 28 09:44:06 ... starting to sort Suffix Array. This may take a long time...
Jun 28 09:44:25 ... sorting Suffix Array chunks and saving them to disk...
Jun 28 09:56:13 ... loading chunks from disk, packing SA...
Jun 28 09:57:33 ... finished generating suffix array
Jun 28 09:57:33 ... generating Suffix Array index
Jun 28 10:01:16 ... completed Suffix Array index
Jun 28 10:01:16 ..... processing annotations GTF
Jun 28 10:01:53 ..... inserting junctions into the genome indices
Jun 28 10:05:15 ... writing Genome to disk ...
Jun 28 10:05:16 ... writing Suffix Array to disk ...
Jun 28 10:05:26 ... writing SAindex to disk
Jun 28 10:05:27 ..... finished successfully
#你可以看到我调用了服务器里128G的内存，线程数用了40，仍然运行了20分钟。线程数越大，越占内存。
运行后，生成了如下的文件：
$ ll
total 29037186
-rw------- 1 fy04 fy04 1207 Jun 28 09:43 chrLength.txt
-rw------- 1 fy04 fy04 4500 Jun 28 09:43 chrNameLength.txt
-rw------- 1 fy04 fy04 3293 Jun 28 09:43 chrName.txt
-rw------- 1 fy04 fy04 2143 Jun 28 09:43 chrStart.txt
-rw------- 1 fy04 fy04 41738007 Jun 28 10:01 exonGeTrInfo.tab
-rw------- 1 fy04 fy04 16839015 Jun 28 10:01 exonInfo.tab
-rw------- 1 fy04 fy04 2526461 Jun 28 10:01 geneInfo.tab
-rw------- 1 fy04 fy04 3208569972 Jun 28 10:05 Genome
-rw------- 1 fy04 fy04 842 Jun 28 10:05 genomeParameters.txt
-rw------- 1 fy04 fy04 68326103 Jun 28 10:05 Log.out
-rw------- 1 fy04 fy04 24786270724 Jun 28 10:05 SA
-rw------- 1 fy04 fy04 1565873619 Jun 28 10:05 SAindex
-rw------- 1 fy04 fy04 10194666 Jun 28 10:01 sjdbInfo.txt
-rw------- 1 fy04 fy04 11054143 Jun 28 10:01 sjdbList.fromGTF.out.tab
-rw------- 1 fy04 fy04 8999291 Jun 28 10:01 sjdbList.out.tab
-rw------- 1 fy04 fy04 13332206 Jun 28 10:01 transcriptInfo.tab
以上是针对小鼠或者人的基因组的构建索引的方法，如果你想对比较小的基因组构建索引，比如：酵母，需要加一个参数：--genomeSAindexNbases 6。
而我在网上也查询了其他的一些参考文章，有的在这个参数上添加的是10。（STAR recommended parameters for S cerevisiae (yeast) genome）
比如我构建的酵母Saccharomyces cerevisiae的索引代码：
$ srun --mem=32G --cpus-per-task=20 --time=2:00:00 --pty STAR \
--runMode genomeGenerate \
--runThreadN 40 \
--genomeDir ./ \
--genomeFastaFiles ./genome_rename.fa \
--genomeSAindexNbases 10 \
--sjdbGTFfile ./Saccharomyces_cerevisiae.R64-1-1.100.gtf
这里需要注意的是，如果你构建索引的时候报错，比如说提示你的染色体名称在fasta和gtf文件里不一致：
我第一次试的时候就出现了这种报错，发现我的fasta文件里染色体都是chr开头的：
$ head genome.fa
>chrI
CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACC
CACACACACACATCCTAACACTACCCTAACACAGCCCTAATCTAACCCTG
GCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTAC
CCTGTCCCATTCAACCATACCACTCCGAACCACCATCCATCCCTCTACTT
ACTACCACTCACCCACCGTTACCCTCCAATTACCCATATCCAACCCACTG
CCACTTACCCTACCATTACCCTACCATCCACCATGACCTACTCACCATAC
TGTTCTTCTACCCACCATATTGAAACGCTAACAAATGATCGTAAATAACA
CACACGTGCTTACCCTACCACTTTATACCACCACCACATGCCATACTCAC
CCTCACTTGTATACTGATTTTACGTACGCACACGGATGCTACAGTATATA
然而我的gtf文件里并没有chr，而是只有编号I, II, III这种。所以要把fasta文件里的修改一下：
$ sed 's/chr//g' genome.fa > genome_rename.fa
$ head genome_rename.fa
>I
CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACC
CACACACACACATCCTAACACTACCCTAACACAGCCCTAATCTAACCCTG
GCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTAC
CCTGTCCCATTCAACCATACCACTCCGAACCACCATCCATCCCTCTACTT
ACTACCACTCACCCACCGTTACCCTCCAATTACCCATATCCAACCCACTG
CCACTTACCCTACCATTACCCTACCATCCACCATGACCTACTCACCATAC
TGTTCTTCTACCCACCATATTGAAACGCTAACAAATGATCGTAAATAACA
CACACGTGCTTACCCTACCACTTTATACCACCACCACATGCCATACTCAC
CCTCACTTGTATACTGATTTTACGTACGCACACGGATGCTACAGTATATA
然后就可以愉快的构建索引了，酵母的索引文件非常快，刷刷的就完事了~
Step 2 - Align RNA-Seq Reads to the genome with STAR比对
首先要明确的是，STAR软件的参数实在是太多了，那么经常用的就那么几个。
下面的代码里，其中--twopassMode参数凭个人爱好是否使用，这个参数是非常非常耗时的，加上这个参数以后的运行时间是正常的两倍。而且非常消耗内存（star的twopassMode问题）。但是比对结果更为精准。
以下是我的代码：
$ srun --mem=128G --cpus-per-task=20 --time=2:00:00 --pty STAR \
--runThreadN 40 \ #线程数
--runMode alignReads \#比对模式
--readFilesCommand zcat \#说明你的fastq文件是压缩形式的，就是.gz结尾的，不加的话会报错
--quantMode TranscriptomeSAM GeneCounts \ #将reads比对至转录本序列
--twopassMode Basic \#先按索引进行第一次比对，而后把第一次比对发现的新剪切位点信息加入到索引中进行第二次比对。这个参数可以保证更精准的比对情况，但是费时也费内存。
--outSAMtype BAM Unsorted \ #输出BAM文件，不进行排序。如果不加这一行，只输出SAM文件。
--outSAMunmapped None \
--genomeDir /gpfs/home/fangy04/downloads/STAR_index/GRCh38/ \ #索引文件目录
--readFilesIn /gpfs/home/fangy04/downloads/SRR8112732_1.fastq.gz /gpfs/home/fangy04/downloads/SRR8112732_2.fastq.gz \ #两个fastq文件目录
--outFileNamePrefix DRB_TT_seq_SRR8112732 #输出文件前缀
然后运行的时候会出现进度显示：
Jun 28 16:44:01 ..... started STAR run
Jun 28 16:44:01 ..... loading genome
Jun 28 16:44:15 ..... started 1st pass mapping
Jun 28 16:51:41 ..... finished 1st pass mapping
Jun 28 16:51:42 ..... inserting junctions into the genome indices
Jun 28 16:54:06 ..... started mapping
Jun 28 17:06:27 ..... finished mapping
Jun 28 17:06:52 ..... finished successfully
我用了服务器里128G只比对了一对fastq文件，用时20分钟，嗯，大家可以体验一下~
STAR比对后会输出好几个文件。上面的代码里，如果你没有加--outSAMtype BAM Unsorted这一个参数的设置，只会输出SAM。在这些文件里，最重要的就是“*.Aligned.out.sam”，你可以使用samtools对其进行bam文件的转换以及排序。但是如果你设置了输出文件为bam，那么直接用samtools进行排序就好了。
输出的文件：
9216920116 Jun 28 17:06 DRB_TT_seq_SRR8112732Aligned.out.bam #这个文件是最重要的，用来后续进行remove duplicates和sort
1166235552 Jun 28 17:06 DRB_TT_seq_SRR8112732Aligned.toTranscriptome.out.bam #这个文件是那些比对到转录本上的reads组成的bam文件
2034 Jun 28 17:06 DRB_TT_seq_SRR8112732Log.final.out
20188 Jun 28 17:06 DRB_TT_seq_SRR8112732Log.out
2571 Jun 28 17:06 DRB_TT_seq_SRR8112732Log.progress.out
1585521 Jun 28 17:06 DRB_TT_seq_SRR8112732ReadsPerGene.out.tab
6732305 Jun 28 17:06 DRB_TT_seq_SRR8112732SJ.out.tab #剪切的信息
8192 Jun 28 16:51 DRB_TT_seq_SRR8112732_STARgenome
8192 Jun 28 16:51 DRB_TT_seq_SRR8112732_STARpass1
如果你想查看比对率之类的信息，可以查看*.Log.final.out文件：
$ cat DRB_TT_seq_SRR8112732Log.final.out
Started job on | Jun 28 16:44:01
Started mapping on | Jun 28 16:54:06
Finished on | Jun 28 17:06:52
Mapping speed, Million of reads per hour | 369.36
Number of input reads | 78590839 #fastq文件的信息
Average input read length | 150 #read长度
UNIQUE READS:
Uniquely mapped reads number | 58096443 #唯一比对上的reads数量
Uniquely mapped reads % | 73.92%
Average mapped length | 144.37
Number of splices: Total | 1880228
Number of splices: Annotated (sjdb) | 1554548 #剪切数
Number of splices: GT/AG | 1513028
Number of splices: GC/AG | 50306
Number of splices: AT/AC | 4877
Number of splices: Non-canonical | 312017 #非典型剪切数
Mismatch rate per base, % | 0.73%
Deletion rate per base | 0.01%
Deletion average length | 1.61
Insertion rate per base | 0.01%
Insertion average length | 1.63
MULTI-MAPPING READS: #多重比对数
Number of reads mapped to multiple loci | 4439909
% of reads mapped to multiple loci | 5.65%
Number of reads mapped to too many loci | 405437
% of reads mapped to too many loci | 0.52%
UNMAPPED READS: #未比对上的reads
Number of reads unmapped: too many mismatches | 0
% of reads unmapped: too many mismatches | 0.00%
Number of reads unmapped: too short | 15355117
% of reads unmapped: too short | 19.54%
Number of reads unmapped: other | 293933
% of reads unmapped: other | 0.37%
CHIMERIC READS: #嵌合的reads数
Number of chimeric reads | 0
% of chimeric reads | 0.00%
而输出文件里的*SJ.out.tab文件，则是splice junctions的一些信息，其中需要注意的是：对于junction的位置信息，STAR则是按照intron的起始和终止位置来定，而其他的一些软件则是按照exon的位置来决定的（参考:比对软件STAR的简单使用）：
$ head DRB_TT_seq_SRR8112732SJ.out.tab
chr1 10078 10178 2 2 0 0 1 26
chr1 10084 10178 2 2 0 0 1 32
chr1 10090 10178 2 2 0 0 1 32
chr1 10102 10178 2 2 1 0 4 37
chr1 10113 10183 2 2 1 0 1 37
chr1 10137 10282 2 2 1 0 1 30
chr1 10137 10288 2 2 1 1 1 37
chr1 10137 10294 2 2 1 0 1 30
chr1 10162 10397 2 2 1 0 1 13
chr1 10179 10225 2 2 1 3 2 27
#第一列：染色体名称
#第二列：内含子起始位点
#第三列：内含子结束位点
#第四列：链。1是正链，2是负链。
#第五列：内含子的motif: 0: 非典型; 1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5:AT/AC, 6: GT/AT
#第六列：0: 未注释的； 1：注释的
#第七列：横跨这个junction上的唯一比对上的reads数
#第八列：横跨这个junction上的muoti-mapping的reads数
#第九列：maximum spliced alignment overhang
以上是个简单的学习笔记，如果想了解更多的更详细的参数设置，可以阅读STAR的官方说明书：https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf（这是最新的版本说明书，是2020年6月更新的）
