hisat2比對 生信筆記 轉錄組分析HISAT2

2021-10-12 05:13:22 字數 4898 閱讀 9553

hisat2是一款快速、敏感的序列比對軟體。使用改進的bwt演算法,相比bowtie/tophat2具有更高的敏感性和更快的運算速度。

安裝hisat2

我採用conda安裝,也是最簡單的方法

conda install -c bioconda hisat2

conda install -c bioconda/label/cf201901 hisat2

解壓:unzip hisat2-2.1.0-beta-linux_x86_64.zip

新增環境變數:vim ~/.bashrc

export path=/user/your_path:$path

儲存退出:source ~/.bashrc

建立索引

建立基因組索引:hisat2-build –p 4 genome.fa genome

建立基因組+轉錄組+snp索引:

hisat2提供了兩個python指令碼將gtf檔案轉換成hisat2-build能使用的檔案

extract_exons.py genome.gtf>genome.exon

extract_splice_sites.py genome.gtf>genome.ss

另外,分析snp可以將其資訊加入至索引中。

extract_snps.py snp.txt>genome.snp

最後,使用將基因組,轉錄組,snp建立索引。

hisat2-build -p 8 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tran_index

不新增轉錄本和snp也可以建立

hisat2-build -p 8 genome.fa genome_index

官方操作手冊簡要

3.1 用法:

hisat2 [options]* -x [-s ]

e. g., 單端測序的比對:

hisat2 -f -x genome_index -u reads_1.fas -s eg1.sam
e.g., 雙端測序的比對:
hisat -f -x genome_index -1 reads_1.fas -2 reads_2.fas -s eg2.sam
3.3 輸入選項:

-q輸入檔案為fastq格式。fastq格式為預設引數。

-qseq

輸入檔案為qseq格式。

-f輸入檔案為fasta格式。

-r輸入檔案中,每一行代表一條序列,沒有序列名和測序質量等。選擇此項時,–ignore-quals引數也會被選擇。

-c此引數後是直接比對的序列,而不是包含序列的檔名。序列間用逗號隔開。選擇此項時,–ignore-quals引數也會被選擇。

-s/–skip

跳過輸入檔案中前條序列進行比對。

-u/–qupto

只使用輸入檔案中前條序列進行比對,預設是沒有限制。

-5/–trim5

比對前去除每條序列5』端個鹼基

-3/–trim3

比對前去除每條序列3』端個鹼基

–phred33

輸入的fastq檔案鹼基質量值編碼標準為phred33,phred33為預設引數。

–phred64

輸入的fastq檔案鹼基質量值編碼標準為phred64。

–solexa-quals

將solexa的鹼基質量值編碼標準轉換為phred。

–int-quals

將sam檔案轉為bam檔案並排序

samtools 安裝: conda install samtools

samtools view -s .eg.sam -b > eg.bam

samtools sort -n  eg.bam eg_sort.bam 因為雙端測序,必須對bam檔案進行排序

counts計數

推薦三個軟體——heseq,featurecount和bedtools

4.1 heseq

4.1.1 heseq安裝

conda install -c bcbio htseq

4.1.2 heseq使用

htseq-count -f bam -t exon -i gene_id -munioneg_sort.bamgenome.gtf > counts.txt

####命令引數

-f | --format default: sam 設定輸入檔案的格式,該引數的值可以是sam或bam。

-r | --order default: name 設定sam或bam檔案的排序方式,該引數的值可以是name或pos。前者表示按read名進行排序,後者表示按比對的參考基因組位置進行排序。若測序資料是雙末端測序,當輸入sam/bam檔案是按pos方式排序的時候,兩端reads的比對結果在sam/bam檔案中一般不是緊鄰的兩行,程式會將reads對的第乙個比對結果放入記憶體,直到讀取到另一端read的比對結果。因此,選擇pos可能會導致程式使用較多的記憶體,它也適合於未排序的sam/bam檔案。而pos排序則表示程式認為雙末端測序的reads比對結果在緊鄰的兩行上,也適合於單端測序的比對結果。很多其它表達量分析軟體要求輸入的sam/bam檔案是按pos排序的,但htseq推薦使用name排序,且一般比對軟體的預設輸出結果也是按name進行排序的。

-s | --stranded default: yes 設定是否是鏈特異性測序。該引數的值可以是yes,no或reverse。no表示非鏈特異性測序;若是單端測序,yes表示read比對到了基因的正義鏈上;若是雙末端測序,yes表示read1比對到了基因正義鏈上,read2比對到基因負義鏈上;reverse表示雙末端測序情況下與yes值相反的結果。根據說明檔案的理解,一般情況下雙末端鏈特異性測序,該引數的值應該選擇reverse(本人暫時沒有測試該引數)。

-a | --a default: 10 忽略比對質量低於此值的比對結果。在0.5.4版本以前該引數預設值是0。

-t | --type default: exon 程式會對該指定的feature(gtf/gff檔案第三列)進行表達量計算,而gtf/gff檔案中其它的feature都會被忽略。

-i | --idattr default: gene_id 設定feature id是由gtf/gff檔案第9列那個標籤決定的;若gtf/gff檔案多行具有相同的feature id,則它們來自同乙個feature,程式會計算這些features的表達量之和賦給相應的feature id。

-m | --mode default: union 設定表達量計算模式。該引數的值可以有union, intersection-strict and intersection-nonempty。這三種模式的選擇請見上面對這3種模式的示意圖。從圖中可知,對於原核生物,推薦使用intersection-strict模式;對於真核生物,推薦使用union模式。

-o | --samout 輸出乙個sam檔案,該sam檔案的比對結果中多了乙個xf標籤,表示該read比對到了某個feature上。

-q | --quiet 不輸出程式執行的狀態資訊和警告資訊。

-h | --help 輸出幫助資訊。

4.2 featurecount <推薦使用>

featurecount通過subread安裝 conda install -y subread

4.2.1 用法

featurecounts -t 6 -p -t exon

-g gene_id -a gencode.vm13.annotation.gtf
-oeg_sort.bam
###引數
4.2.2 資料提取
則只需要提取出第1和第7列的資訊。
cut -f 1,7 count.txt |grep -v '^#' >feacturcounts.txt
4.2.2 資料矯正計算公式

4.3 bedtools用法:bedtools multicov -bams 1.bam 2.bam 3.bam 4.bam-bedfile.bed >read.count.txt

比對軟體hisat2的使用

官方手冊 基因組比對軟體常用bwa,轉錄組比對軟體常用bowtie2 hisat2等,其中有參考基因組的常用hisat2,沒有參考基因組的常用bowtie2。下面我們來介紹一下hisat2的使用方法 一 建立索引 建立基因組索引 hisat2 build p 4 genome.fa genome 建...

轉 A2W W2A T2A T2W 等巨集

如果你覺得使用 widechartomultibyte,multibytetowidechar 等函式比較麻煩 眾多的引數,緩衝區的分配與銷毀等。那麼可以使用 a2w w2a t2a t2w 等巨集來代替,它們對上面兩個函式進行了封裝。在使用這些巨集之前,應該包含標頭檔案 atlconv.h 並在呼...

2 基礎控制項2

transform 的預設值為 1,0,0,1,0,0 nslog nsstringfrom 可以列印其他型別的 比如類 affine transform a ffine transform 的初始化為 make make 只能改變一次 不可重複改變 如果想要重複改變 則使用去掉 make 的方法t...