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 -munion
eg_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-bed
file
.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...