转录组分析 | 使用STAR进行比对

前几期，小编已经教大家完成了RNA-seq数据的质控，下面就要正式开始转录组分析啦！

通过二代测序我们可以获得150bp左右的reads，如果想要知道reads是从哪个转录本上测出来的，就需要将reads比对到参考基因组上。比对的算法很复杂，但简单理解就是看reads与基因组上哪个区域一致。

常用的比对工具有Tophat2、Hisat2和STAR。不同的工具有各自的优势，目前比较流行的工具是Hisat2和STAR，它俩的比对速度都比较快，STAR的uniquely mapping reads比例较高，对于我们做多倍体物种分析的人来说，STAR的优势非常大，所以小编以STAR为例教大家进行reads比对。

## 下载 STAR
wget -c https://github.com/alexdobin/STAR/archive/2.7.3a.tar.gz
## 解压 STAR
tar -xvzf 2.7.3a.tar.gz
## 运行 STAR
./STAR-2.7.3a/bin/Linux_x86_64/STAR

在进行reads比对前，我们需要先构建基因组索引。

## 构建基因组序列
STAR --runThreadN 6 
    --runMode genomeGenerate 
    --genomeDir index_dir 
    --genomeFastaFiles genome.fasta 
    --sjdbGTFfile genome.gtf 
    --sjdbOverhang 149

--runThreadN：线程数。

--runMode genomeGenerate：构建基因组索引。

--genomeDir：索引目录。（index_dir一定要是存在的文件夹，需提前建好）

--genomeFastaFiles：基因组文件。

--sjdbGTFfile：基因组注释文件。

--sjdbOverhang：reads长度减1。

索引构建完成后，就可以看到index_dir中生成了以下文件：

有了索引后，我们就可以进行reads比对了。

## 进行 reads 比对
STAR --twopassMode Basic 
    --quantMode TranscriptomeSAM GeneCounts 
    --runThreadN 6 
    --genomeDir index_dir 
    --alignIntronMin 20 
    --alignIntronMax 50000 
    --outSAMtype BAM SortedByCoordinate 
    --sjdbOverhang 149 
    --outSAMattrRGline ID:sample SM:sample PL:ILLUMINA 
    --outFilterMismatchNmax 2 
    --outSJfilterReads Unique 
    --outSAMmultNmax 1 
    --outFileNamePrefix out_prefix 
    --outSAMmapqUnique 60 
    --readFilesCommand gunzip -c 
    --readFilesIn seq1.fq.gz seq2.fq.gz

--twopassMode Basic：使用two-pass模式进行reads比对。简单来说就是先按索引进行第一次比对，而后把第一次比对发现的新剪切位点信息加入到索引中进行第二次比对。

--quantMode TranscriptomeSAM GeneCounts：将reads比对至转录本序列。

--runThreadN：线程数。

--genomeDir：索引目录。

--alignIntronMin：最短的内含子长度。（根据GTF文件计算）

--alignIntronMax：最长的内含子长度。（根据GTF文件计算）

--outSAMtype BAM SortedByCoordinate：输出BAM文件并进行排序。

--sjdbOverhang：reads长度减1。

--outSAMattrRGline：ID代表样本ID，SM代表样本名称，PL为测序平台。在使用GATK进行SNP Calling时同一SM的样本可以合并在一起。

--outFilterMismatchNmax：比对时允许的最大错配数。

--outSJfilterReads Unique：对于跨越剪切位点的reads（junction reads)，只考虑跨越唯一剪切位点的reads。

--outSAMmultNmax：每条reads输出比对结果的数量。

--outFileNamePrefix：输出文件前缀。

--outSAMmapqUnique 60：将uniquely mapping reads的MAPQ值调整为60，满足下游使用GATK进行分析的需要。

--readFilesCommand：对FASTQ文件进行操作。

--readFilesIn：输入FASTQ文件的路径。

比对完成后，我们可以看到输出目录下有以下文件：

我们可以使用samtools查看生成的BAM文件。

## 查看 BAM 文件
samtools view CK-1_Aligned.sortedByCoord.out.bam |head -n 5

可以看到，以"Aligned.sortedByCoord.out.bam"为后缀的BAM文件中，reads比对到的位置是基因组位置。

以"Aligned.toTranscriptome.out.bam"为后缀的BAM文件中，reads比对到的位置是转录本位置。

"Log.final.out"里记录了许多比对情况的统计信息。

STAR的参数非常多，大家在实际应用过程中可以参考它的Manual。

虽然STAR的uniquely mapping reads比例比较高，但运行时所需的内存非常大，大家在使用时一定要注意提供足够大的内存。

参考资料：

https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf