仍然是兩年前的筆記1. prepare-reference
如果用rsem對比對後的bam進行轉錄本定量，則在比對過程中要確保比對用到的索引是由rsem-prepare-reference產生的。
~/software/rsem/rsem-prepare-reference \
--transcript-to-gene-map ~/project/rna-seq/ref_cds/gene_transcript.txt \ #作用是在後面的定量結果檔案中，新增gene名稱, 轉錄本名稱兩列，該檔案每一行都是gene_id\ttranscript_id的形式，eg: cluster_11236 cluster_11236.1
--bowtie2 \ #rsem可呼叫bowtie, bowtie2, star三種比對工具；這裡選用bowtie2
可以看到，單純用bowtie2建的索引和rsem呼叫bowtie2建的索引是不一樣的。
2. calculate-expression
用法分為兩類，分別是從fa/fq得到表達矩陣，和從sam/bam得到表達矩陣（仍然要求是比對到rsem-prepare-reference生成的索引）。以單端的fq資料為例。
rsem-calculate-expression [options] upstream_read_file(s) reference_name sample_name
rsem-calculate-expression [options] --paired-end upstream_read_file(s) downstream_read_file(s) reference_name sample_name
rsem-calculate-expression [options] --sam/--bam [--paired-end] input reference_name sample_name
cat ~/project/rna-seq/dir.txt | while read id
do~/software/rsem/rsem-calculate-expression -p 8 --bowtie2 \
~/project/data/rna-seq/$.fastq.gz \
~/project/rna-seq/ref_cds/cds.byrsem \
--samtools-sort-mem 2g --fragment-length-mean 50 \ # 單端資料建議使用--fragment-length-mean和--fragment-length-sd
完成之後得到這些檔案，其中，rsem.genes.results和rsem.isoforms.results分別表示gene水平和轉錄本水平的定量結果。每一列含義：
less rsem.genes.results|head -n 1
gene_id transcript_id(s) length effective_length expected_count tpm fpkm
less rsem.isoforms.results|head -n 1
transcript_id gene_id length effective_length expected_count tpm fpkm isopct
後面用ebseq檢驗差異基因/轉錄本時，會使用到這兩個檔案。
3. differential expression analysis using ebseq
下面以gene水平差異分析為例。
3.1 generate-data-matrix
這一步提取上一步得到的每個樣本定量結果檔案中的expected_count列，組成資料矩陣。
呼叫ebseq進行檢驗
~/software/rsem/rsem-run-ebseq \
genemat.txt 2,2,2,2 genemat.results #2,2,2,2表示4個condition, 每個condition有兩個重複；順序要和3.1中輸入檔案表示的condition的順序一致
#會得到三個檔案
genemat.results.condmeans genemat.results genemat.results.pattern
#genemat.results.pattern
"c1" "c2" "c3" "c4"
"pattern1" 1 1 1 1
"pattern2" 1 1 1 2
"pattern3" 1 1 2 1
"pattern4" 1 1 2 2
"pattern5" 1 2 1 1
"pattern6" 1 2 1 2
"pattern7" 1 2 2 1
"pattern8" 1 2 2 2
"pattern9" 1 1 2 3
"pattern10" 1 2 1 3
"pattern11" 1 2 2 3
"pattern12" 1 2 3 1
"pattern13" 1 2 3 2
"pattern14" 1 2 3 3
"pattern15" 1 2 3 4
#以pattern14為例，1 2 3 3表示某基因表達：c1與c2不同，c3與c4相同
#四種condition如果有基因表達存在差異，就這些情況了
#genemat.results
#第一列是各個基因名稱，接著15列是該基因符合該種parttern的概率
#"map"為該基因最可能的模式；"ppde"：posterior probability of being differentially expressed，越大越好
"pattern1" "pattern2" "pattern3" "pattern4" "pattern5" "pattern6" "pattern7" "pattern8" "pattern9" "pattern10" "pattern11" "pattern12" "pattern13" "pattern14" "pattern15" "map" "ppde"
#genemat.results.condmeans
#為每個樣本合併重複之後的定量結果，如下圖，這個結果可以用來控制fold change
控制fdr(錯誤發現率)來挑選差異基因
~/software/rsem/rsem-control-fdr \
genemat.results 0.05 genemat.de.txt
將genemat.results檔案中，ppde大於0.95的記錄提取出來
因水平有限，有錯誤的地方，歡迎批評指正！
如何使用GMAP GSNAP進行轉錄組序列比對
gmap最早用於講est cdna序列比對到參考基因組上，可以用於基因組結構注釋。後來高通量測序時代，又開發了gsnap支援高通量資料比對，這篇文章主要介紹gmap，畢竟高通量轉錄組資料比對大家更喜歡用star,hista2等軟體。下面是我原始碼安裝的 wget tar xf gmap gsnap ...
如何使用GMAP GSNAP進行轉錄組序列比對
gmap最早用於講est cdna序列比對到參考基因組上，可以用於基因組結構注釋。後來高通量測序時代，又開發了gsnap支援高通量資料比對，這篇文章主要介紹gmap，畢竟高通量轉錄組資料比對大家更喜歡用star,hista2等軟體。下面是我原始碼安裝的 wget tar xf gmap gsnap ...
使用Mikado挑選最好的轉錄本進行注釋
mikado是基於python3寫的基因組結構注釋工具，它主要做的是從多個轉錄組組裝工具得到的轉錄本裡挑選出最好的結果作為基因組的結構注釋。此外，它還會基於同源蛋白比對結果對轉錄本打分。換句話說這個軟體主要是根據轉錄組資料進行注釋，沒有 ab inito 軟體安裝比較方法，我們可以使用biocond...