歡迎關(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ì)。
## 下載 STARwget -c https://github.com/alexdobin/STAR/archive/2.7.3a.tar.gz## 解壓 STARtar -xvzf 2.7.3a.tar.gz## 運(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)“生信小王子”乡革！
