發(fā)現(xiàn)服務(wù)器上沒有安裝STAR (Spliced Transcripts Alignment to a Reference)，這個(gè)轉(zhuǎn)錄組最常用的比對(duì)工具之一，也是我之前一直的用的轉(zhuǎn)錄組比對(duì)工具，今天安裝一下并重新學(xué)習(xí)，好好理解之前設(shè)置的參數(shù)是否正確。

STAR是ENCODE計(jì)劃(ENCyclopedia Of DNA Elements,人類基因組DNA元件百科全書計(jì)劃)的御用pipeline工具，在轉(zhuǎn)錄組的文章中出鏡率極高，別人說其準(zhǔn)確率高，映射速度快，但需要占用大量內(nèi)存，對(duì)計(jì)算資源有較高的要求。在之前Hisat2安裝使用過程中，提到了2017年的一篇NC比較轉(zhuǎn)錄組比對(duì)工具的文章，又查了一下，這樣總結(jié)的：STAT相比較TopHat和Hisat2，有較高的唯一比對(duì)率；STAR會(huì)將沒有paired mapping上的reads都剔除，避免single reads比對(duì)到基因組上；并且STAR對(duì)lower-quality（包括more soft-clipped和錯(cuò)配堿基）比對(duì)有較高的容忍度，這對(duì)一些雜合率較高的基因組優(yōu)勢(shì)比較明顯；這次注意到，在用GATK對(duì)RNA-Seq進(jìn)行 Call Variants時(shí)，采用STAR的STAR 2-pass模式，估計(jì)以后也會(huì)用到。

下載安裝軟件

https://github.com/alexdobin/STAR

選擇其中一個(gè)版本下載后， tar -zxvf 進(jìn)行解壓：

?tar -zxvf STAR-2.7.9a.tar.gz

?cd STAR/source

?make STAR

然后這次我注意到在bin目錄下有兩個(gè)帶有linux目錄及source目錄下都有STAR命令，都可以運(yùn)行，我翻看之前的命令行，用的第二個(gè)里面的STAR命令，初步判斷三個(gè)均可以，這次還是選用2中的STAR命令：? ? ? ?

二、構(gòu)建基因組索引Index

和Hisat2一樣，需要先構(gòu)建基因組索引，索引文件作用現(xiàn)在還只記得是在比對(duì)過程中，我們并不是把幾十萬條reads直接比對(duì)到基因組上去，而是和Index進(jìn)行比較，使比對(duì)過程變地可行，希望等課題結(jié)束后，再回過頭來好好學(xué)習(xí)一下索引文件作用的原理，先上腳本：

參數(shù)解釋：

--runThreadN：線程數(shù)為10

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

--genomeDir：指定索引生成目錄；

--genomeFastaFiles：指定參考基因組；

--sjdbGTFfile：指定參考基因組的注釋文件；

--sjdbOverhang：這個(gè)是reads長度的最大值減1，默認(rèn)是100，我不是很理解很多人分析的學(xué)習(xí)方法中都設(shè)置100，二代測序都是150bp的序列長度，我設(shè)置了149 (有時(shí)間時(shí)改一下數(shù)值比較一下對(duì)結(jié)果是否有影響);

發(fā)現(xiàn)有三個(gè)反斜杠“\”異常成了黃色，暫時(shí)不清楚原因，結(jié)果報(bào)錯(cuò)了：

其實(shí)我也不知道為啥，將運(yùn)行命令行的反斜杠去掉，再試一下：

剛才的問題解決了，又報(bào)了其它錯(cuò)誤信息：

居然是gtf文件的錯(cuò)誤，第一次遇見這個(gè)問題，然后找原因：

我們看一下gtf的開頭是CM023448.1，如下圖：

我的參考基因組開頭是>GWHAMMI00000001，如下圖：

原來是染色體的命名方式不一樣，一個(gè)是CM開頭，另一個(gè)是GWHAMMI開頭，我回到NCBI去下載序列文件又看了一下,居然是我之前下錯(cuò)文件了(從另一個(gè)數(shù)據(jù)庫下載的參考基因組，兩個(gè)數(shù)據(jù)庫同一物種染色體編號(hào)規(guī)則不同)，之前做的工作又浪費(fèi)了，重新下載，指定序列文件，30min后，成功建立索引，索引目錄如下：

reads比對(duì)：

相比于Hisat2，STAR太多的參數(shù)設(shè)置了，對(duì)于模式生物還好，很多默認(rèn)參數(shù)就可以，但對(duì)于我的課題研究，就得仔細(xì)看看這些參數(shù)了，著實(shí)用去了我不少時(shí)間，先上我的腳本，如下圖：

我的參數(shù)設(shè)置：

未用的其它參數(shù)：

--outFilterMismatchNmax：比對(duì)時(shí)允許的最大錯(cuò)配數(shù)(可根據(jù)結(jié)果修改);

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

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

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

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

--alignIntronMin：最短的內(nèi)含子長度設(shè)定了20，(根據(jù)GTF文件計(jì)算);

--alignIntronMax：最長的內(nèi)含子長度設(shè)定了50000，(根據(jù)GTF文件計(jì)算);

--bamRemoveDuplicatesType?? 輸出BAM文件時(shí)，STAR還可以對(duì)BAM進(jìn)行一些預(yù)處理，用于去重。

四：結(jié)果如下圖，

1、使用samtools查看生成的BAM文件。

samtoolsview sample_Aligned.sortedByCoord.out.bam |head -n 5

2、結(jié)果內(nèi)容：

Aligned.sortedByCoord.out.bam:reads比對(duì)到基因組的位置;

Aligned.toTranscriptome.out.bam:reads比對(duì)到轉(zhuǎn)錄本的位置；

Log.final.out：統(tǒng)計(jì)了比對(duì)情況的信息,是非常重要的結(jié)果;

SJ.out.tab:splice junctions的一些信息，其中需要注意的是：對(duì)于junction的位置信息，STAR則是按照intron的起始和終止位置來定，而其他的一些軟件則是按照exon的位置來決定的

?

附：我比較了一下star和Hisat2的結(jié)果差異，在運(yùn)行時(shí)間和比對(duì)率上，star并沒有表現(xiàn)出明顯的優(yōu)越性上。

參考：

https://blog.csdn.net/weixin_28913137/article/details/112281831

本文使用 文章同步助手 同步