歡迎關(guān)注微信公眾號(hào)“生信小王子”！

前幾期，小編已經(jīng)教大家完成了RNA-seq數(shù)據(jù)的質(zhì)控，下面就要正式開始轉(zhuǎn)錄組分析啦！

通過(guò)二代測(cè)序我們可以獲得150bp左右的reads，如果想要知道reads是從哪個(gè)轉(zhuǎn)錄本上測(cè)出來(lái)的，就需要將reads比對(duì)到參考基因組上。比對(duì)的算法很復(fù)雜，但簡(jiǎn)單理解就是看reads與基因組上哪個(gè)區(qū)域一致。

常用的比對(duì)工具有Tophat2、Hisat2和STAR。不同的工具有各自的優(yōu)勢(shì)，目前比較流行的工具是Hisat2和STAR，它倆的比對(duì)速度都比較快，STAR的uniquely mapping reads比例較高，對(duì)于我們做多倍體物種分析的人來(lái)說(shuō)，STAR的優(yōu)勢(shì)非常大，所以小編以STAR為例教大家進(jìn)行reads比對(duì)。

## 下載 STAR

wget -c https://github.com/alexdobin/STAR/archive/2.7.3a.tar.gz

## 解壓 STAR

tar -xvzf 2.7.3a.tar.g

z## 運(yùn)行 STAR

./STAR-2.7.3a/bin/Linux_x86_64/STAR

在進(jìn)行reads比對(duì)前，我們需要先構(gòu)建基因組索引。

##?構(gòu)建基因組索引

STAR --runThreadN 6 --runMode genomeGenerate --genomeDir index_dir --genomeFastaFiles genome.fasta --sjdbGTFfile genome.gtf --sjdbOverhang 149

--runThreadN：線程數(shù)。

--runMode genomeGenerate：構(gòu)建基因組索引。

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

--genomeFastaFiles：基因組文件。

--sjdbGTFfile：基因組注釋文件。

--sjdbOverhang：reads長(zhǎng)度減1。

索引構(gòu)建完成后，就可以看到index_dir中生成了以下文件：

有了索引后，我們就可以進(jìn)行reads比對(duì)了。

## 進(jìn)行 reads 比對(duì)

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模式進(jìn)行reads比對(duì)。簡(jiǎn)單來(lái)說(shuō)就是先按索引進(jìn)行第一次比對(duì)，而后把第一次比對(duì)發(fā)現(xiàn)的新剪切位點(diǎn)信息加入到索引中進(jìn)行第二次比對(duì)。

--quantMode TranscriptomeSAM GeneCounts：將reads比對(duì)至轉(zhuǎn)錄本序列。

--runThreadN：線程數(shù)。

--genomeDir：索引目錄。

--alignIntronMin：最短的內(nèi)含子長(zhǎng)度。（根據(jù)GTF文件計(jì)算）

--alignIntronMax：最長(zhǎng)的內(nèi)含子長(zhǎng)度。（根據(jù)GTF文件計(jì)算）

--outSAMtype BAM SortedByCoordinate：輸出BAM文件并進(jìn)行排序。

--sjdbOverhang：reads長(zhǎng)度減1。

--outSAMattrRGline：ID代表樣本ID，SM代表樣本名稱，PL為測(cè)序平臺(tái)。在使用GATK進(jìn)行SNP Calling時(shí)同一SM的樣本可以合并在一起。

--outFilterMismatchNmax：比對(duì)時(shí)允許的最大錯(cuò)配數(shù)。

--outSJfilterReads Unique：對(duì)于跨越剪切位點(diǎn)的reads（junction reads)，只考慮跨越唯一剪切位點(diǎn)的reads。

--outSAMmultNmax：每條reads輸出比對(duì)結(jié)果的數(shù)量。

--outFileNamePrefix：輸出文件前綴。

--outSAMmapqUnique 60：將uniquely mapping reads的MAPQ值調(diào)整為60，滿足下游使用GATK進(jìn)行分析的需要。

--readFilesCommand：對(duì)FASTQ文件進(jìn)行操作。

--readFilesIn：輸入FASTQ文件的路徑。

比對(duì)完成后，我們可以看到輸出目錄下有以下文件：

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

##?查看 BAM 文件

samtools?view?CK-1_Aligned.sortedByCoord.out.bam?|head?-n 5

可以看到，以"Aligned.sortedByCoord.out.bam"為后綴的BAM文件中，reads比對(duì)到的位置是基因組位置。

以"Aligned.toTranscriptome.out.bam"為后綴的BAM文件中，reads比對(duì)到的位置是轉(zhuǎn)錄本位置。

"Log.final.out"里記錄了許多比對(duì)情況的統(tǒng)計(jì)信息。

STAR的參數(shù)非常多，大家在實(shí)際應(yīng)用過(guò)程中可以參考它的Manual。

雖然STAR的uniquely mapping?reads比例比較高，但運(yùn)行時(shí)所需的內(nèi)存非常大，大家在使用時(shí)一定要注意提供足夠大的內(nèi)存。

參考資料：

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

歡迎關(guān)注微信公眾號(hào)“生信小王子”！