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格式是目前用來存放大量核酸比對結(jié)果信息的通用格式，也是人類能夠“直接”閱讀的格式類型当编，而BAM和CRAM是為了方便傳輸届慈，降低存儲壓力將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. 索引
大于號：將一條命令執(zhí)行結(jié)果（標(biāo)準(zhǔn)輸出臊泌，或者錯誤輸出，本來都要打印到屏幕上面的）重定向其它輸出設(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
//雖然我不知道第三步有什么用礁凡。高氮。。
如果你要比較同一個樣本(within-sample)不同基因之間的表達(dá)情況顷牌，你就需要考慮到轉(zhuǎn)錄本長度剪芍，因?yàn)檗D(zhuǎn)錄本越長，那么檢測的片段也會更多窟蓝，直接比較等于讓小孩和大人進(jìn)行賽跑罪裹。如果你是比較不同樣本（across sample）同一個基因的表達(dá)情況，雖然不必在意轉(zhuǎn)錄本長度运挫，但是你要考慮到測序深度（sequence depth)状共，畢竟測序深度越高，檢測到的概率越大谁帕。除了這兩個因素外峡继，你還需要考慮GC%所導(dǎo)致的偏差，以及測序儀器的系統(tǒng)偏差匈挖。目前對read count標(biāo)準(zhǔn)化的算法有RPKM（SE）, FPKM（PE）碾牌，TPM, TMM等，不同算法之間的差異與換算方法已經(jīng)有文章進(jìn)行整理和吐槽了儡循。
在轉(zhuǎn)錄本水平上舶吗，一般常用工具為Cufflinks和它的繼任者StringTie， eXpress择膝。這些軟件要處理的難題就時轉(zhuǎn)錄本亞型（isoforms）之間通常是有重疊的誓琼，當(dāng)二代測序讀長低于轉(zhuǎn)錄本長度時，如何進(jìn)行區(qū)分调榄？這些工具大多采用的都是expectation maximization（EM）踊赠。好在我們有三代測序。上述軟件都是alignment-based每庆，目前許多alignment-free軟件筐带，如kallisto, silfish, salmon，能夠省去比對這一步缤灵，直接得到read count伦籍，在運(yùn)行效率上更高蓝晒。不過最近一篇文獻(xiàn)[1]指出這類方法在估計豐度時存在樣本特異性和讀長偏差。
-f bam/sam： 指定輸入文件格式帖鸦，默認(rèn)SAM
-r name/pos: 你需要利用samtool sort對數(shù)據(jù)根據(jù)read name或者位置進(jìn)行排序芝薇，默認(rèn)是name
-s yes/no/reverse: 數(shù)據(jù)是否來自于strand-specific assay。DNA是雙鏈的作儿，所以需要判斷到底來自于哪條鏈洛二。如果選擇了no， 那么每一條read都會跟正義鏈和反義鏈進(jìn)行比較攻锰。默認(rèn)的yes對于雙端測序表示第一個read都在同一個鏈上晾嘶，第二個read則在另一條鏈上。
-a 最低質(zhì)量娶吞， 剔除低于閾值的read
-m 模式 union（默認(rèn)）, intersection-strict and intersection-nonempty垒迂。一般而言就用默認(rèn)的，作者也是這樣認(rèn)為的妒蛇。
-i id attribute: 在GTF文件的最后一欄里机断，會有這個基因的多個命名方式（如下）， RNA-Seq數(shù)據(jù)分析常用的是gene_id绣夺， 當(dāng)然你可以寫一個腳本替換成其他命名方式吏奸。
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
