新媒體經(jīng)理

在所有是非的場景中，我最喜歡你。

大寫字母

這次給大家?guī)淼氖枪俜降膶?duì)比軟件STAR。 該項(xiàng)目是法國國家人類基因組研究所（NHGRI）于2003年9月發(fā)起的公共聯(lián)合研究項(xiàng)目，致力于找出人類基因組中的所有功能。 成分[。 這是人類基因組計(jì)劃完成后，國家人類基因組研究所發(fā)起的最重要的項(xiàng)目之一。 本項(xiàng)目形成的所有數(shù)據(jù)將及時(shí)發(fā)布在公共數(shù)據(jù)庫中。

在我之前的RNA-seq數(shù)據(jù)分析——方法論文章的實(shí)踐文章中，STAR在對(duì)比軟件的對(duì)比中也表現(xiàn)出了不錯(cuò)的表現(xiàn)。 所以在進(jìn)行比較時(shí)，我也考慮過將其與STAR一起使用來檢查它們的性能并選擇合適的比較工具。

明星安裝

cd biosoft && mkdir STAR && cd STAR
wget https://github.com/alexdobin/STAR/archive/2.5.3a.tar.gz
tar -xzf 2.5.3a.tar.gz
cd STAR-2.5.3a
# for easy use, add bin/ to your PATH

下載需要參考基因組和索引創(chuàng)建

# downloading dna index fasta file
nohup wget -r -np -nH -nd -R index.html -L ftp://ftp.ensembl.org/pub/release-90/fasta/homo_sapiens/dna_index/ &
# download gft annotation file
nohup wget ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens/Homo_sapiens.GRCh38.90.chr_patch_hapl_scaff.gtf.gz &
mkdir STAR_index && cd STAR_index
STAR --runMode genomeGenerate --genomeDir ~/reference/STAR_index/ --genomeFastaFiles ~/reference/genome/hg38/Homo_sapiens.GRCh38.dna.toplevel.fa --sjdbGTFfile ~/reference/genome/hg38/Homo_sapiens.GRCh38.90.chr_patch_hapl_scaff.gtf --sjdbOverhang 199
# --sjdbOverhang 數(shù)值為reads長度-1
# Mode 為generate
# --genomeFastaFiles --sjdbGTFfile 分別對(duì)應(yīng)fasta文件和GTF文件

STAR的使用

# STAR的manual里面給了最基本的比對(duì)參數(shù)示例
STAR
--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2 ]
# 基本示例，
針對(duì)fastq.gz文件增加--readFilesCommand gunzip -c 參數(shù)/--readFilesCommand zcat參數(shù)，針對(duì)bzip2文件使用--readFilesCommand bunzip2 -c參數(shù)
STAR --runThreadN 20 --genomeDir ~/reference/STAR_index/ --readFilesCommand zcat --readFilesIn ~/RNA-seq/LiuPing_data/RNA-seq/SC_w2q20m35_N_1.fq.gz ~/RNA-seq/LiuPing_data/RNA-seq/SC_w2q20m35_N_2.fq.gz
# 輸出unsorted or sorted bam file
--outSAMtype BAM Unsorted 實(shí)際上就是-name 的sort，下游可以直接接HTSeq
--outSAMtype BAM SortedByCoordinate
--outSAMtype BAM Unsorted SortedByCoordinate 兩者都輸出

附加參數(shù)說明

# 單獨(dú)指定注釋文件，而不用在構(gòu)建的時(shí)候使用
--sjdbGTFfile /path/to/ann.gtf
--sjdbFileChrStartEnd /path/to/sj.tab
# ENCODE參數(shù)
# 減少偽junction的幾率
--outFilterType BySJout
# 最多允許一個(gè)reads被匹配到多少個(gè)地方
--outFilterMultimapNmax 20
# 在未有注釋的junction區(qū)域，最低允許突出多少個(gè)bp的單鏈序列
--alignSJoverhangMin 8
# 在有注釋的junction區(qū)域，最低允許突出多少個(gè)bp的單鏈序列
--alignSJDBoverhangMin 1
# 過濾掉每個(gè)paired read mismatch數(shù)目超過N的數(shù)據(jù)，999代表著忽略這個(gè)過濾
--outFilterMismatchNmax 999
# 相對(duì)paired read長度可以允許的mismatch數(shù)目，如果read長度為100，數(shù)值設(shè)定為0.04，則會(huì)過濾掉100*2*0.04=8個(gè)以上的數(shù)據(jù)
--outFilterMismatchNoverReadLmax 0.04
# 最小的intro長度
--alignIntronMin 20
# 最大的intro長度
--alignIntronMax 1000000
# maximum genomic distance between mates，翻譯不出來，自行理解
--alignMatesGapMax 1000000

星輸出

STAR可以根據(jù)您的參數(shù)設(shè)置輸出多個(gè)結(jié)果文件，包括各類信息。 下面是帶有默認(rèn)參數(shù)的輸出文件的詳細(xì)顯示。 對(duì)于一些比較難的翻譯dnastar序列比對(duì)，我選擇使用原文

E00516:168:H37WKCCXY:8:1101:6400:59130 ? ?99 ? ?1 ? ?92836373 ? ?255 ? ?20M1063N129M ? ?= ? ?92837548 ? ?4244 ? ?GGCTTGTCTATCCCTCACAGTACCAAACGATTCCCTGGTTATGATTCTGAAAGCAAGGAATTTAATGCAGAAGTACATCGGAAGCACATCATGGGCCAGAATGTTGCAGATTACATGCGCTACTTAATGGAAGAAGATGAAGATGCTTA ? ?AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ ? ?NH:i:1 ? ?HI:i:1 ? ?AS:i:289 ? ?nM:i:0
# 我截取了一條比對(duì)信息
我們來看一下最后面的 NH:i:1 ?HI:i:1 ?AS:i:289 ? ?nM:i:0
NH:i:后面的數(shù)值代表著此條read比對(duì)到幾個(gè)loci，1代表著unique map，數(shù)值大于1代表著multi-mappers
HI:i:后面的數(shù)值attrbiutes enumerates multiple
alignments of a read starting with 1，下游分析接cufflinks or stringtie的時(shí)候需要使用參數(shù)--outSAMattrIHstart 0
AS:i:的數(shù)值代表著local alignment score (paired for paired-edn reads)
nM:i:的數(shù)值代表著the number of mismatches per (paired) alignment, not to be confused with NM, which is the number of mismatches in each mate
關(guān)于下游處理工具的兼容性還需要使用者自己仔細(xì)參考manual

tail Log.out
Joined thread # 12
Completed: thread #13
Joined thread # 13
Joined thread # 14
Joined thread # 15
Joined thread # 16
Joined thread # 17
Joined thread # 18
Joined thread # 19
ALL DONE!

tail Log.progress.out
Sep 08 17:57:52 ? ? 33.1 ? ?23115987 ? ? ?285 ? ?94.1% ? ?284.0 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 17:58:53 ? ? 34.0 ? ?24349711 ? ? ?285 ? ?94.1% ? ?284.0 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:00:23 ? ? 33.5 ? ?24789186 ? ? ?285 ? ?94.1% ? ?284.1 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:01:51 ? ? 33.3 ? ?25493588 ? ? ?285 ? ?94.1% ? ?284.0 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:02:58 ? ? 33.5 ? ?26284824 ? ? ?285 ? ?94.1% ? ?284.1 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:04:23 ? ? 33.7 ? ?27163519 ? ? ?285 ? ?94.1% ? ?284.1 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:05:36 ? ? 33.1 ? ?27428080 ? ? ?285 ? ?94.1% ? ?284.1 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:06:54 ? ? 33.8 ? ?28659661 ? ? ?285 ? ?94.1% ? ?284.1 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:08:00 ? ? 34.3 ? ?29741743 ? ? ?285 ? ?94.1% ? ?283.9 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
ALL DONE!

head Log.progress.out 
 ? ? ? ? ? Time ? ?Speed ? ? ? ?Read ? ? Read ? Mapped ? Mapped ? Mapped ? Mapped Unmapped Unmapped Unmapped Unmapped
 ? ? ? ? ? ? ? ? ? ?M/hr ? ? ?number ? length ? unique ? length ? MMrate ? ?multi ? multi+ ? ? ? MM ? ?short ? ?other
Sep 08 17:17:47 ? ? ?2.9 ? ? ? 88583 ? ? ?288 ? ?94.2% ? ?287.4 ? ? 0.1% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.7% ? ? 0.0%
Sep 08 17:18:53 ? ? 14.5 ? ? ?711158 ? ? ?282 ? ?94.1% ? ?281.9 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
Sep 08 18:08:00 ? ? 34.3 ? ?29741743 ? ? ?285 ? ?94.1% ? ?283.9 ? ? 0.2% ? ? 4.0% ? ? 0.1% ? ? 0.0% ? ? 1.8% ? ? 0.0%
ALL DONE!

head SJ.out.tab 
1 ? ?14830 ? ?14969 ? ?2 ? ?2 ? ?0 ? ?1 ? ?9 ? ?69
1 ? ?14844 ? ?14969 ? ?2 ? ?2 ? ?0 ? ?0 ? ?2 ? ?30
1 ? ?15039 ? ?15795 ? ?2 ? ?2 ? ?1 ? ?2 ? ?7 ? ?53
1 ? ?15948 ? ?16606 ? ?2 ? ?2 ? ?1 ? ?1 ? ?1 ? ?41
1 ? ?16028 ? ?16606 ? ?2 ? ?2 ? ?0 ? ?0 ? ?1 ? ?57
1 ? ?16311 ? ?16606 ? ?2 ? ?2 ? ?0 ? ?2 ? ?0 ? ?67
1 ? ?16766 ? ?16853 ? ?2 ? ?2 ? ?0 ? ?2 ? ?0 ? ?43
1 ? ?16766 ? ?16857 ? ?2 ? ?2 ? ?1 ? ?17 ? ?108 ? ?73
1 ? ?16766 ? ?16875 ? ?2 ? ?2 ? ?0 ? ?0 ? ?1 ? ?61
1 ? ?16789 ? ?16875 ? ?2 ? ?2 ? ?0 ? ?0 ? ?1 ? ?53
# 參數(shù)釋義
column 1: chromosome
column 2: first base of the intron (1-based)
column 3: last base of the intron (1-based)
column 4: strand (0: undened, 1: +, 2: -)
column 5: intron motif: 0: non-canonical; 1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5:AT/AC, 6: GT/AT
column 6: 0: unannotated, 1: annotated (only if splice junctions database is used)
column 7: number of uniquely mapping reads crossing the junction
column 8: number of multi-mapping reads crossing the junction
column 9: maximum spliced alignment overhang

寫在最后

也許我探索STAR的最終目的是借助STAR的and來實(shí)現(xiàn)的。 我自己處理的數(shù)據(jù)上有一個(gè)-，其余的對(duì)比軟件還沒有發(fā)現(xiàn)這個(gè)功能。

使用-參數(shù)時(shí)，STAR可以將reads拆分成兩部分，分別進(jìn)行比較

STAR-是一個(gè)，可以進(jìn)行STARdnastar序列比對(duì)，點(diǎn)我看代碼

其實(shí)STAR也可以做2-pass，比較新穎

使用-參數(shù)也可以達(dá)到HTSeq的療效。 它可以幫助您生成計(jì)數(shù)，從而節(jié)省您 HTSeq 的精力。 有時(shí)間再回去做一下比較，看看HTSeq和HTSeq的效率如何。

參考：

日常鮑勃聯(lián)排別墅

如有侵權(quán)請(qǐng)聯(lián)系刪除！