hisat2 -t -x ~/rna_seq_analysis/reference/index/hg19/genome -1 ~/rna_seq_analysis/fastq/SRR5894154_1.fastq.gz -2 ~/rna_seq_analysis/fastq/SRR5894154_2.fastq.gz -S ~/rna_seq_analysis/aligned/SRR5894154.sam
之前出問題好像是因?yàn)樯倏樟丝崭瘛?/p>
http://blog.sciencenet.cn/blog-3334560-1078097.html
得到的sam格式文件40多G。。。。
我哭了。。。
下一步把sam文件轉(zhuǎn)化為bam文件，用samtools
SAM格式是目前用來存放大量核酸比對(duì)結(jié)果信息的通用格式，也是人類能夠“直接”閱讀的格式類型，而BAM和CRAM是為了方便傳輸，降低存儲(chǔ)壓力將SAM進(jìn)行壓縮得到的格式形式。?注，BAM格式必須要建立索引才能快速讀取指定位置的信息。
# 1.3版本前
samtools view -bS bwa.sam > bwa.bam
samtools sort bwa.bam > bwa_sorted.bam
samtools index bwa_sorted.bam
# 1.3版本后
samtools sort bwa.sam > bwa_sorted.bam
samtools index bwa_sorted.bam
1. 格式轉(zhuǎn)換
2. 排序
3. 索引
大于號(hào)：將一條命令執(zhí)行結(jié)果（標(biāo)準(zhǔn)輸出，或者錯(cuò)誤輸出，本來都要打印到屏幕上面的）重定向其它輸出設(shè)備（文件，打開文件操作符，或打印機(jī)等等）
1. samtools view -S SRR5894154.sam -b > SRR5894154.bam? ?
bam文件不到8g，于是趕緊把sam文件刪了～
2.?samtools sort SRR5894154.bam > SRR5894154_sorted.bam
3. samtools index SRR5894154_sorted.bam
//雖然我不知道第三步有什么用。。。
如果你要比較同一個(gè)樣本(within-sample)不同基因之間的表達(dá)情況，你就需要考慮到轉(zhuǎn)錄本長(zhǎng)度，因?yàn)檗D(zhuǎn)錄本越長(zhǎng)，那么檢測(cè)的片段也會(huì)更多，直接比較等于讓小孩和大人進(jìn)行賽跑。如果你是比較不同樣本（across sample）同一個(gè)基因的表達(dá)情況，雖然不必在意轉(zhuǎn)錄本長(zhǎng)度，但是你要考慮到測(cè)序深度（sequence depth)，畢竟測(cè)序深度越高，檢測(cè)到的概率越大。除了這兩個(gè)因素外，你還需要考慮GC%所導(dǎo)致的偏差，以及測(cè)序儀器的系統(tǒng)偏差。目前對(duì)read count標(biāo)準(zhǔn)化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之間的差異與換算方法已經(jīng)有文章進(jìn)行整理和吐槽了。
在轉(zhuǎn)錄本水平上，一般常用工具為Cufflinks和它的繼任者StringTie， eXpress。這些軟件要處理的難題就時(shí)轉(zhuǎn)錄本亞型（isoforms）之間通常是有重疊的，當(dāng)二代測(cè)序讀長(zhǎng)低于轉(zhuǎn)錄本長(zhǎng)度時(shí)，如何進(jìn)行區(qū)分？這些工具大多采用的都是expectation maximization（EM）。好在我們有三代測(cè)序。上述軟件都是alignment-based，目前許多alignment-free軟件，如kallisto, silfish, salmon，能夠省去比對(duì)這一步，直接得到read count，在運(yùn)行效率上更高。不過最近一篇文獻(xiàn)[1]指出這類方法在估計(jì)豐度時(shí)存在樣本特異性和讀長(zhǎng)偏差。
-f bam/sam： 指定輸入文件格式，默認(rèn)SAM
-r name/pos: 你需要利用samtool sort對(duì)數(shù)據(jù)根據(jù)read name或者位置進(jìn)行排序，默認(rèn)是name
-s yes/no/reverse: 數(shù)據(jù)是否來自于strand-specific assay。DNA是雙鏈的，所以需要判斷到底來自于哪條鏈。如果選擇了no， 那么每一條read都會(huì)跟正義鏈和反義鏈進(jìn)行比較。默認(rèn)的yes對(duì)于雙端測(cè)序表示第一個(gè)read都在同一個(gè)鏈上，第二個(gè)read則在另一條鏈上。
-a 最低質(zhì)量， 剔除低于閾值的read
-m 模式 union（默認(rèn)）, intersection-strict and intersection-nonempty。一般而言就用默認(rèn)的，作者也是這樣認(rèn)為的。
-i id attribute: 在GTF文件的最后一欄里，會(huì)有這個(gè)基因的多個(gè)命名方式（如下）， RNA-Seq數(shù)據(jù)分析常用的是gene_id， 當(dāng)然你可以寫一個(gè)腳本替換成其他命名方式。
gene_id "ENSG00000223972.5_2"; transcript_id "ENST00000456328.2_1"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-002"; exon_number 2; exon_id "ENSE00003582793.1_1"; level 2;.
htseq-count -s no?-r pos? -f bam ~/rna_seq_analysis/aligned/SRR5894154_sorted.bam ~/rna_seq_analysis/human_genome/gencode.v28lift37.annotation.sorted.gff3 > ~/rna_seq_analysis/aligned/SRR5894154.count
