要求
實(shí)現(xiàn)這個(gè)功能的軟件也很多影暴，還是煩請(qǐng)大家先自己搜索幾個(gè)教程型宙，入門請(qǐng)統(tǒng)一用htseq-count伦吠，對(duì)每個(gè)樣本都會(huì)輸出一個(gè)表達(dá)量文件毛仪。
需要用腳本合并所有的樣本為表達(dá)矩陣。參考：生信編程直播第四題：多個(gè)同樣的行列式文件合并起來
對(duì)這個(gè)表達(dá)矩陣可以自己簡單在excel或者R里面摸索腺逛，求平均值刨晴，方差狈癞。
看看一些生物學(xué)意義特殊的基因表現(xiàn)如何茂契，比如GAPDH,β-ACTIN等等掉冶。
理論基礎(chǔ)
在上篇的比對(duì)中，我們需要糾結(jié)是否真的需要比對(duì)恢共，如果你只需要知道已知基因的表達(dá)情況璧亚，那么可以選擇alignment-free工具（例如salmon, sailfish)癣蟋，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）濒生。到了這一篇的基因（轉(zhuǎn)錄本）定量罪治，需要考慮的因素就更加多了，以至于我不知道如何說清才能理清邏輯恒序。
定量分為三個(gè)水平
- 基因水平(gene-level)
- 轉(zhuǎn)錄本水平(transcript-level)
- 外顯子使用水平(exon-usage-level)谁撼。
在基因水平上厉碟，常用的軟件為HTSeq-count，featureCounts崭参，BEDTools, Qualimap, Rsubread, GenomicRanges等何暮。以常用的HTSeq-count為例铐殃，這些工具要解決的問題就是根據(jù)read和基因位置的overlap判斷這個(gè)read到底是誰家的孩子富腊。值得注意的是不同工具對(duì)multimapping reads處理方式也是不同的，例如HTSeq-count就直接當(dāng)它們不存在是整。而Qualimpa則是一人一份浮入，平均分配羊异。
對(duì)每個(gè)基因計(jì)數(shù)之后得到的count matrix再后續(xù)的分析中球化，要注意標(biāo)準(zhǔn)化的問題筒愚。如果你要比較同一個(gè)樣本(within-sample)不同基因之間的表達(dá)情況，你就需要考慮到轉(zhuǎn)錄本長度句伶，因?yàn)檗D(zhuǎn)錄本越長考余，那么檢測(cè)的片段也會(huì)更多，直接比較等于讓小孩和大人進(jìn)行賽跑疫蔓。如果你是比較不同樣本（across sample）同一個(gè)基因的表達(dá)情況身冬，雖然不必在意轉(zhuǎn)錄本長度酥筝，但是你要考慮到測(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)行整理和吐槽了。但是看铆，有一些下游分析的軟件會(huì)要求是輸入的count matrix是原始數(shù)據(jù)弹惦，未經(jīng)標(biāo)準(zhǔn)化，比如說DESeq2石抡，這個(gè)時(shí)候你需要注意你上一步所用軟件會(huì)不會(huì)進(jìn)行標(biāo)準(zhǔn)化啰扛。
在轉(zhuǎn)錄本水平上嗡贺，一般常用工具為Cufflinks和它的繼任者StringTie暑刃， eXpress。這些軟件要處理的難題就時(shí)轉(zhuǎn)錄本亞型（isoforms）之間通常是有重疊的溜嗜，當(dāng)二代測(cè)序讀長低于轉(zhuǎn)錄本長度時(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í)存在樣本特異性和讀長偏差惰匙。
在外顯子使用水平上，其實(shí)和基因水平的統(tǒng)計(jì)類似铃将。但是值得注意的是為了更好的計(jì)數(shù)，我們需要提供無重疊的外顯子區(qū)域的gtf文件[2]劲阎。用于分析差異外顯子使用的DEXSeq提供了一個(gè)Python腳本（dexseq_prepare_annotation.py）執(zhí)行這個(gè)任務(wù)绘盟。
小結(jié)
計(jì)數(shù)分為三個(gè)水平： gene-level, transcript-level, exon-usage-level
標(biāo)準(zhǔn)化方法： FPKM RPKM TMM TPM
輸出表達(dá)矩陣
在RNA-Seq分析中，每一個(gè)基因就是一個(gè)feature（特征？）奥此，而基因被認(rèn)為是它的所有外顯子的和集弧哎。在可變剪切分析中，可以單獨(dú)把每個(gè)外顯子當(dāng)作一個(gè)feature稚虎。而在ChIP-Seq分析中撤嫩，feature則是預(yù)先定義的結(jié)合域。但是確定一個(gè)read到底屬于哪一個(gè)feature有時(shí)會(huì)非常棘手蠢终。因此HTSeq提供了三種模式序攘，示意圖見前一幅圖
- the union of all the sets S(i) for mode union. This mode is recommended for most use cases.
- the intersection of all the sets S(i) for mode intersection-strict.
- the intersection of all non-empty sets S(i) for mode intersection-nonempty.
基本用法非常的簡單：
# 安裝
conda install htseq
# 使用
# htseq-count [options] <alignment_file> <gtf_file>
htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count
用一個(gè)循環(huán)處理多個(gè)BAM文件(在/mnt/f/Data目錄下)
mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log
done
運(yùn)行的時(shí)間會(huì)比較久，所以可以去了解不同參數(shù)的用法了寻拂，其中比較常用的為：
- -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;...
Jimmy的伏筆
我們這次分析是人類mRNA-Seq測(cè)序的結(jié)果宗侦，但是我們其實(shí)只下載了3個(gè)sra文件愚臀。一般而言RNA-Seq數(shù)據(jù)分析都至少要有2個(gè)重復(fù)，所以必須要有4個(gè)sra文件才行矾利。我在仔細(xì)讀完文章的方法這一段以后，發(fā)現(xiàn)他們有一批數(shù)據(jù)用的是其他課題組的： For 293 cells, the mRNA-seq results of the control samples include (1) those from the doxycycline-treated parental Flp-In T-REx 293 cells by us and (2) those from the doxycycline-treated control Flp-In T-REx 293 cells performed by another group unrelated to us (sample GSM1095127 in GSE44976)馋袜。 然后和Jimmy交流之后男旗，他也承認(rèn)自己只分析了小鼠的數(shù)據(jù)，而沒有分析人類的數(shù)據(jù)欣鳖。所以我們需要根據(jù)文章提供的線索下載另外一份數(shù)據(jù)察皇，才能進(jìn)行下一步的分析。
這個(gè)時(shí)候就有一個(gè)經(jīng)常被問到的問題：不同來源的RNA-Seq數(shù)據(jù)能夠直接比較嗎？甚至說如果不同來源的RNA-seq數(shù)據(jù)的構(gòu)建文庫都不一樣該如何比較?不同來源的RNA-Seq結(jié)果之間比較需要考慮 批次效應(yīng)（batch effect) 的影響什荣。
處理批次效應(yīng)矾缓，根據(jù)我搜索的結(jié)果，是不能使用FPKM/RPKM稻爬，關(guān)于這個(gè)標(biāo)準(zhǔn)化的吐槽嗜闻，我在biostars上找到了如下觀點(diǎn)：
- FPKM/RPKM 不是標(biāo)準(zhǔn)化的方法，它會(huì)引入文庫特異的協(xié)變量
- FPKM/RPKM has never been peer-reviewed, it has been introduced as an ad-hoc measure in a supplementary 沒有同行評(píng)審
- One of the authors of this paper states, that it should not be used because of faulty arithmetic 作者說算法有問題
- All reviews so far have shown it to be an inferior scale for DE analysis of genes Length normalization is mostly dispensable imo in DE analysis because gene length is constant
有人建議使用一個(gè)Bioconductor包http://www.bioconductor.org/packages/devel/bioc/html/sva.html 我沒有具體了解桅锄，有生之年去了解補(bǔ)充琉雳。
還有人引用了一篇文獻(xiàn) IVT-seq reveals extreme bias in RNA-sequencing 證明不同文庫的RNA-Seq結(jié)果會(huì)存在很大差異。
結(jié)論： 可以問下原作者他們是如何處理數(shù)據(jù)的友瘤，居然有一個(gè)居然沒有重復(fù)的分析也能過審翠肘。改用小鼠數(shù)據(jù)進(jìn)行分析”柩恚或者使用無重復(fù)的分析方法束倍，或者模擬一份數(shù)據(jù)出來，先把流程走完盟戏。
合并表達(dá)矩陣
HTSeq-count輸出結(jié)果是一個(gè)個(gè)獨(dú)立的文件绪妹，后續(xù)分析需要把多個(gè)文件合并成一個(gè)行為基因名，列為樣本名抓半，中間為count的行列式文件喂急。肯定是不會(huì)用Excel手動(dòng)處理（雖然可以寫一個(gè)VBA腳本笛求，但是數(shù)據(jù)量過大不好處理了）廊移，這里使用的Python寫一個(gè)腳本。
基本邏輯：
- 讀取文件
- 建立一個(gè)字典探入，如果key不在字典中狡孔，新增key和value，如果key在字典中蜂嗽，追加value苗膝。
- 輸出
#!/usr/bin/python3
import sys
mydict = {}
for file in sys.argv[1:]:
for line in open(file, 'r'):
key,value = line.strip().split('\t')
if key in mydict:
mydict[key] = mydict[key] + '\t' + value
else:
mydict[key] = value
for key,value in mydict.items():
print(key + '\t' + value +'\n')
這幾行代碼寫了2個(gè)番茄鐘，但是debug花了我一個(gè)番茄鐘植旧。問題出在str和list兩種數(shù)據(jù)格式的混亂使用辱揭。還有一個(gè)bug： 由于詞典是無序的，所以原本代表樣本來源的第一行病附，會(huì)跑到其他行问窃。
在論壇上找到一個(gè)更加簡潔的代碼（要求基因名順序排列），用到paste,awk,printf這三個(gè)shell命令完沪。
paste *.txt | awk '{printf $1 "\t";for(i=2;i<=NF;i+=2) printf $i"\t";printf $i}'
保存為countCombiner.py域庇，輸入文件為count, 輸出為標(biāo)準(zhǔn)輸出嵌戈，需要重定向。
簡單分析
這一步需要用到R語言或者是Excel讀取數(shù)據(jù)听皿。
1.導(dǎo)入數(shù)據(jù)
options(stringsAsFactors = FALSE)
# import data if sample are small
control <- read.table("F:/Data/RNA-Seq/matrix/SRR3589956.count",
sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589957.count",
sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589958.count",
sep="\t",col.names = c("gene_id","rep2"))
2.數(shù)據(jù)整合熟呛。gencode的注釋文件中的gene_id（如ENSG00000105298.13_3）在EBI是不能搜索到的，所以我就只保留ENSG00000105298這部分尉姨。
# merge data and delete the unuseful row
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL
3.總體情況, 大部分基因都為0庵朝，所以可以刪掉節(jié)省體積。
summary(raw_count_filt)
control rep1 rep2
Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
Median : 0.0 Median : 0.0 Median : 0.0
Mean : 356.1 Mean : 370.3 Mean : 316.6
3rd Qu.: 15.0 3rd Qu.: 15.0 3rd Qu.: 10.0
Max. :161867.0 Max. :121902.0 Max. :105565.0
4.看幾個(gè)具體基因
在EBI上搜索GAPDH找到ID為ENSG00000111640啊送。
GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]
gene_id control rep1 rep2
ENSG00000111640 ENSG00000111640.14_2 41857 53902 55302
文章研究的AKAP95（ENSG00000105127）的表達(dá)量在KD中都降低了
> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> AKAP95
gene_id control rep1 rep2
ENSG00000105127 ENSG00000105127.8_2 1168 539 506
下面的差異基因表達(dá)偿短，讓我想想，該如何收拾Jimmy挖的坑馋没。
參考文獻(xiàn)
[1] Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
[2] Detecting differential usage of exons from RNA-seq data
[3] RNA-seq Data Analysis-A Practical Approach(2015)
